disp(' ')
disp(' composite_beam.m   ver 1.0   November 15, 2010')
disp(' by Tom Irvine   Email: tomirvine@aol.com')
disp(' ')
%
disp(' This program calculates the natural frequency of multi-layer beam. ');
disp(' ');
%
clear width;
clear thick;
clear inertia;
clear area;
clear AE;
clear EI;
clear Yw;
%
disp(' Select system ');
disp(' 1=English  2=metric ');
iunits=input(' ');
disp(' ');
n=input(' Enter the number of layers ');
%
disp(' ');
disp(' Start with the top layer in the following data entry .');
disp(' ');
%
for(i=1:n)
    out1=sprintf('\n ** Layer %d **\n',i);
    disp(out1);
%
    imat=input(' Enter material: \n 1=aluminum  2=steel  3=other ');
%
    if(iunits==1)
        if(imat==1)
            rho(i)=0.1;
            E(i)=1.0e+07;
        end
        if(imat==2)
            rho(i)=0.29;
            E(i)=3.0e+007;
        end
    else
        if(imat==1)
            rho(i)=2768;
            E(i)=6.891e+010;
        end
        if(imat==2)
            rho(i)=7887;
            E(i)=2.067e+011;
        end        
    end
%
    if(iunits==1)
       if(imat==3)
           rho(i)=input(' Enter mass density (lbm/in^3)   ');
             E(i)=input(' Enter the elastic modulus (psi) ');
       end
       disp(' ');
       width(i)=input(' Enter width (inch)     ');
       thick(i)=input(' Enter thickness (inch) ');
    else
       if(imat==3)
           rho(i)=input(' Enter mass density (kg/m^3)  ');
           E(i)=input(' Enter the elastic modulus (Pa) ');
       end
       disp(' ');
       width(i)=input(' Enter width (mm)     ');
       thick(i)=input(' Enter thickness (mm) ');
%       
       width(i)=width(i)/1000;
       thick(i)=thick(i)/1000; 
    end
    area(i)=width(i)*thick(i);
    inertia(i)=(1/12)*width(i)*(thick(i))^3;
    AE(i)=area(i)*E(i);
    EI(i)=E(i)*inertia(i);
end
%
num=0.;
den=0.;
for(i=1:n)
    ut=0;
    for(j=i+1:n)
        ut=ut+thick(j);
    end    
    Yw(i)=ut+thick(i)/2;
end
for(i=1:n)
    num=num+AE(i)*Yw(i);
    den=den+AE(i);
end
Ybar=num/den;
%
if(iunits==1)
    out1=sprintf('\n Ybar = %8.4g inch  (AE centroid) \n',Ybar);
else
    out1=sprintf('\n Ybar = %8.4g mm (AE centroid) \n',Ybar/1000);
end
disp(out1);
%
for(i=1:n)
    c(i)=abs(Yw(i)-Ybar);
    c2(i)=(c(i))^2;
    AEc2(i)=AE(i)*c2(i);
end
%
AEc2_sum=sum(AEc2);
EI_sum=sum(EI);
%
if(iunits==1)
    out1=sprintf('\n  AEc2 = %8.4g lbf in^2 ',AEc2_sum);
    disp(out1);
    out1=sprintf('\n    EI = %8.4g lbf in^2 \n',EI_sum);
    disp(out1);
else
    out1=sprintf('\n  AEc2 = %8.4g N m^2 ',AEc2_sum);
    disp(out1);
    out1=sprintf('\n    EI = %8.4g N m^2 \n',EI_sum);
    disp(out1);   
end
%
composite_EI=AEc2_sum+EI_sum;
%
if(iunits==1)
    out1=sprintf('\n composite EI = %8.4g lbf in^2 \n',composite_EI);
else
    out1=sprintf('\n composite EI = %8.4g N m^2 \n',composite_EI);
end
disp(out1);
%
rho_sum=sum(area.*rho);
%
if(iunits==1)
    out1=sprintf('\n composite rho = %8.4g lbm/in \n',rho_sum);
    rho_sum=rho_sum/386;
else
    out1=sprintf('\n composite rho = %8.4g kg/m   \n',rho_sum);
end
disp(out1);
%
EI_d_rho=composite_EI/rho_sum;
%
disp(' ');
if(iunits==1)
    L=input(' Enter length (inch)' );
else
    L=input(' Enter length (m)' );    
end
%
disp(' ');
disp(' Enter boundary conditions ');
disp('  1=free-free ');
disp('  2=fixed-free ');
disp('  3=fixed-fixed ');
disp('  4=pinned-pinned ');
ibc=input(' ');
%
R=sqrt(EI_d_rho)/L^2;
%
if(ibc==1 || ibc==3 )  % free-free or fixed-fixed
    fn = (22.373/(2*pi))*R; 
end
%
if(ibc==2)  % fixed-free
    fn = 3.5156/(2*pi)*R;     
end
%
if(ibc==4)  % pinned-pinned
    fn = (pi/2)*R;     
end
%
out1=sprintf('\n fn=%8.4g Hz \n',fn);
disp(out1);