%
%   rocket_bending.m,  ver 1.1, July 9, 2004
%
%   by Tom Irvine
%   Email: tomirvine@aol.com
%
%   This program calculates the natural frequencies of a rocket modeled as
%   a beam
%
disp(' ref 1 ');
%
clear ibc;
clear ne;
clear total_length;
clear xv;
clear dx;
clear dof;
clear ipm;
clear pm;
clear xpm;
clear klocal;
clear mlocal;
clear dof;
clear mode;
clear eigenvalue;
clear eigv;
clear fn;
clear emm;
clear mmm;
clear pff;
clear pf;
clear EI;
clear rho;
clear stiff;
clear mass;
clear v;
%
disp(' ');
disp(' The EI file must have two columns:  x(inch) & EI (lbf in^2) ');
EIF = input(' Enter the EI filename:  ');
%
disp(' ');
disp(' The mass file must have two columns:  x(inch) & mass/length (lbm/in) ');
MF = input(' Enter the mass filename:  ');
%
pm=0.;
%
i=1;
[EI,rho,ibc,ne,total_length,xv,dx,dof,ipm,pm,xpm]=rocket_bending_data_entry(EIF,MF);
%
[klocal] = rocket_local_stiffness(EI,ne,dx);
[mlocal] = rocket_local_mass(rho,ne,dx);
%
[stiff,mass,dof]=assembly(ne,klocal,mlocal,ibc,dof,ipm,pm,xpm,dx);
%
[mode,eigenvalue]=eig(stiff,mass);
%
disp(' ');
disp(' Eigenvalues ');
%
for i=1:dof
    if(abs(eigenvalue(i,i))<1.0e-22)
        eigenvalue(i,i)=0.;
    end
    eigv(i)=eigenvalue(i,i);
end
%
%  The eigenvalues should already be sorted.
%
disp(' ');
for i=1:dof
  fn(i)=sqrt(eigv(i))/(2.*pi);
end
for i=1:rank(mass)
  v(i)=1.;
end      
%
%  The eigenvectors should already be normalized with respect to the mass
%  matrix.  Verify that modes'*mass*modes = identity matrix.
%
LM=mode'*mass*v';
disp('        Natural    Participation    Effective  ');
disp('Mode   Frequency      Factor        Modal Mass ');
%
pf=LM;
sum=0;
%
%  Only half of the values are printed since accuracy decreases with mode
%  number.
%
mmm=mode'*mass*mode;
for i=1:dof/2
    pff(i)=pf(i)/mmm(i,i);
    emm(i)=pf(i)^2/mmm(i,i);
    if(emm(i) < rho*total_length/100000.)
        emm(i)=0.;
        pff(i)=0.;
    end
    out1 = sprintf('%d  %10.4g Hz   %10.4g   %10.4g',i,fn(i),pff(i),emm(i) );
    disp(out1)
    sum=sum+pf(i)^2;
end
%
disp(' ')
out1 = sprintf(' Total Effective Modal Mass =   = %10.4g lbf sec^2/in',sum );
disp(out1)
%
total_mass=0.;
for i=1:ne
    total_mass=total_mass+rho(i)*dx(i);
end
%
disp(' ')
out1 = sprintf(' Total Mass = %10.4g lbf sec^2/in',total_mass );
disp(out1)
out1 = sprintf('            = %10.4g lbm',total_mass*386. );
disp(out1)
%
[pmode,xvv]=plot_modes(ibc,mode,dof,dx);