disp(' ');
disp(' ss_beam_mass.m   ver 1.0  Feb 15, 2011 ');
disp(' ');
disp(' by Tom Irvine  Email: tomirvine@aol.com ');
disp(' ');
disp(' This program calculates natural frequencies of ');
disp(' a simply-supported beam with a discrete mass via');
disp(' Rayleigh-Ritz Method. ');
disp(' ');
%
clear m;
clear k;
clear fn;
%
pm=0.;
xpm=0.;
disp(' ');
L =input(' Enter the total length (in) ');
disp(' ');
L1=input(' Enter the length to the discrete mass (in) ');
disp(' ');
icross=input(' Enter the cross-section: \n 1=rectangle  2=solid cylinder  3=other ');
disp(' ');
if(icross==1)
			width=input(' Enter the width (inch) ');
			thick=input(' Enter the thickness (inch) ');
			A=width*thick;
			I=(1./12.)*width*(thick^3.);
end
if(icross==2)
			diam=input(' Enter the diameter (inch) ');
			A=(pi/4.)*(diam^2.);
			I=(pi/64.)*(diam^4.);
end
if(icross ~=1 & icross ~=2)
            disp(' ');
            I = input(' Enter the area moment of inertia (in^4) ');
            %
            disp(' ');
            A = input(' Enter the cross-section area (in^2) ');
end
%
disp(' ');
imat=input(' Enter material: \n 1=aluminum  2=steel  3=other ');
%
if(imat==1)
    rho=0.1;
	E=1.0e+07;
end
if(imat==2)
    rho=0.29;
    E=3.0e+007;
end
if(imat ~=1 & imat ~=2)
    disp(' ');   
    E = input(' Enter the elastic modulus (lbf/in^2) ');
    disp(' ');
    rho = input(' Enter the mass density (lbm/in^3) ');           
 end
%
tm=rho*A*L;
%
rho=rho*A/386.;
%
disp(' ');
M=input(' Enter the discrete mass (lbm) ');
%
disp(' ');
out1=sprintf('  Beam Mass = %8.4g lbm ',tm);
out2=sprintf(' Total Mass = %8.4g lbm ',tm+M);
%
disp(out1);
disp(out2);
%
M=M/386.;
%
PL=pi/L;
PL2=PL^2;
PL3=PL^3;
PL4=PL^4;
PI=pi;
TwoP=2*PI;
ThreeP=3*PI;
FourP=4*PI;
FiveP=5*PI;
SixP=6*PI;
%
LR=L1/L;
PL1R=pi*LR;
%
r2=rho/2;
%
c11=  L-(L/TwoP)*sin(TwoP);
c12=  0;
c13=  0;
c22=  L-(L/FourP)*sin(FourP);
c23=  0;
c33=  L-(L/SixP)*sin(SixP);
%
m(1,1)= M*(sin(PL1R))^2               + r2*c11;
m(1,2)= M*(sin(PL1R))*(sin(2*PL1R))   + r2*c12;
m(1,3)= M*(sin(PL1R))*(sin(3*PL1R))   + r2*c13;
m(2,2)= M*(sin(2*PL1R))^2             + r2*c22;
m(2,3)= M*(sin(2*PL1R))*(sin(3*PL1R)) + r2*c23;
m(3,3)= M*(sin(3*PL1R))^2             + r2*c33;
%
k=zeros(3,3);
%
k(1,1)=   c11;
k(2,2)=16*c22;
k(3,3)=81*c33;
%
k=0.5*PL4*E*I*k;
%
for(i=1:3)
    for(j=i+1:3)
        m(j,i)=m(i,j);
        k(j,i)=k(i,j);    
    end
end
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen(k,m,2);
%
disp(' ');
disp('   n     fn(Hz) ');
for i=1:3
    out1 = sprintf('   %d  %8.4g ',i,fn(i));
    disp(out1);
end