disp(' ');
disp(' rod_FEA.m  ver 1.3   March 10, 2012 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the longitudinal vibration frequencies');
disp(' of a rod.  It also calculates the free vibration response to  '); 
disp(' initial velocity. ');
%
clear PF;
clear EMM;
clear A;
clear E;
clear h;
clear rho;
clear T;
clear S2;
clear M2;
clear Kj;
clear Mj;
clear S;
clear M;
%
clear acc;
clear dis;
clear vel;
clear k;
clear m;
clear dam;
clear damp;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear omegan;
clear fn;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
clear t;
clear f1;
clear f2;
clear x;
clear d;
clear v;
clear na;
clear nv;
clear nd;
clear Q;
clear mass;
clear stiffness;
clear max;
%
clear A;
clear B;
clear n;
clear nv;
clear x;
clear v;
clear t;
clear tt;
clear omegad;
clear acc;
clear AF;
clear VF;
clear XF;
%   
tpi=2.*pi;
%
iu=1;  % English units
%
disp(' ');
total_length = input(' Enter the length (inch) ');
diam=input(' Enter the diameter (inch) ');
A=(pi/4.)*(diam^2.);
Area=A;
%    
disp(' ');
[E,rho,mu]=materials(iu);
%
m=rho*A;
%
ml=m*total_length;
%
disp(' ');
out1=sprintf(' mass=%8.4g lbm ',ml*386);
out2=sprintf('     =%8.4g lbf sec^2/in ',ml);
disp(out1);
disp(out2);
%
disp(' ');
ne = input(' Enter the number of elements ');
nn=ne+1;
%
disp(' ');
disp(' Enter boundary conditions ');
disp('  1=free-free ');
disp('  2=fixed-free ');
disp('  3=fixed-fixed ');
ibc=input(' ');
%
h=total_length/ne;
%
Kj=(E*A/h)*[1 -1; -1 1];
Mj=(h*m/6)*[2 1; 1 2];
%
S=zeros(nn,nn);
M=zeros(nn,nn);
%
for(i=1:ne)
    T=zeros(2,nn);
    T(1,i)=1;
    T(2,i+1)=1;
    S=S+T'*Kj*T;
    M=M+T'*Mj*T;
end
%
if(ibc==1) % free-free
    M2=M;
    S2=S;
end
if(ibc==2) % fixed-free
    M2=M(2:nn,2:nn);
    S2=S(2:nn,2:nn);
end
if(ibc==3) % fixed-fixed
    nm=nn-1;
    M2=M(2:nm,2:nm);
    S2=S(2:nm,2:nm);
end
%
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(S2,M2,3);
Q=ModeShapes;
%
nnn=max(size(fn));
%
if ibc==2 || ibc==3
    mass_hat=MST*M2*ModeShapes;
%
    r=ones(nnn,1);
%
    Lbar=MST*M2*r;
    for i=1:nnn
        PF(i)=Lbar(i)/mass_hat(i,i);
        EMM(i)=(Lbar(i))^2/mass_hat(i,i);
    end
%
    disp(' ');
    disp('  PF = Participation Factor ');
    disp(' EMM = Effective Modal Mass ');
    disp(' ');
    disp('  n    fn(Hz)      PF         EMM ');
%
    for i=1:nnn
        out1=sprintf(' %d\t  %7.5g\t  %7.4g\t  %7.4g',i,fn(i),PF(i),EMM(i));
        disp(out1);
    end 
    disp(' ');
    TE=sum(EMM);
    out1=sprintf(' Total EMM = %8.4g lbf sec^2/in',TE);
    disp(out1);
%
end
%
disp(' ');
disp(' Perform free vibration analysis for initial velocity? ');
disp(' 1=yes  2=no ');
fv=input(' ');
%
if(fv==1)
   disp(' ');
   dam=input(' Enter common, viscous damping ratio ');
   disp(' ');
   v=input(' Enter common, initial velocity (in/sec) ');
   disp(' ');
   T=input(' Enter duration (sec) ');
   disp(' ');      
   mf=max(fn);
   qqq=sprintf('Suggested sample rate > %8.4g ',mf);
   disp(qqq);
   sr=input(' Enter sample rate (samples/sec) ');
   dt=1/sr;
%
   DI= zeros(nnn,1);
   VI=v*ones(nnn,1);
%
   sz=size(DI);
   if(sz(2)>sz(1))
       DI=DI';
   end
%
   sz=size(VI);
   if(sz(2)>sz(1))
       VI=VI';
   end
%
   nDI=MST*(M2*DI);
   nVI=MST*(M2*VI);
%   
   damp=dam*ones(nn,1);
%
   iu=1;
%
   nt=round(T/dt);
%
   clear domegan;
   clear omd2;
   for(i=1:nnn)
      omegad(i)=omegan(i)*sqrt(1-(damp(i))^2);
      domegan(i)=damp(i)*omegan(i);
      omd2(i)=(omegad(i))^2;
   end
%
   clear x;
   clear v;
   clear acc;
%
   for(i=1:nt)
       t=dt*(i-1);
       tt(i)=t;
%
       for(j=1:nnn)
           arg=omegad(j)*t;
           eee=exp(-domegan(j)*t);
           dee=-domegan(j)*eee;
%
           cos_arg=cos(arg); 
           sin_arg=sin(arg);        
%           
           A=nDI(j);
           B=(domegan(j)*nDI(j)+nVI(j))/omegad(j);  
%
            n(j)= eee*(A*cos_arg+B*sin_arg);
           nv(j)= dee*(A*cos_arg+B*sin_arg)+omegad(j)*eee*(-A*sin_arg+B*cos_arg);
           na(j)= -(omegan(j))^2*n(j) - 2*damp(j)*omegan(j)*nv(j);  
       end    
       x(i,:)=Q*n';
       v(i,:)=Q*nv';
     acc(i,:)=Q*na'; 
%
   end
%
   disp(' ');
   npl=input(' Enter node for plotting ');
%
   clear tx;
   clear tv;
   clear tacc;
%
   tx=x;
   tv=v;
   tacc=acc;
%
     x=zeros(nt,nn);
     v=zeros(nt,nn);
   acc=zeros(nt,nn);
%
   if(ibc==1) % free-free
       x=tx;
       v=tv;
       acc=tacc;
   end
   if(ibc==2) % fixed-free
         x(:,2:nn)=tx;
         v(:,2:nn)=tv;
       acc(:,2:nn)=tacc;
   end
   if(ibc==3) % fixed-fixed
         x(:,2:nm)=tx;
         v(:,2:nm)=tv;
       acc(:,2:nm)=tacc;     
   end
%
   figure(1);
   plot(tt,x(:,npl));
   grid on;
   xlabel('Time(sec)');
   title('Displacement');
   if(iu==1)
       ylabel('Disp(in)');
   else
       ylabel('Disp(m)');    
   end
%
   figure(2);
   plot(tt,v(:,npl));
   grid on;
   xlabel('Time(sec)');
   title('Velocity');
   if(iu==1)
       ylabel('Vel(in/sec)');
   else
       ylabel('Vel(m/sec)');    
   end
%
   figure(3);
   plot(tt,acc(:,npl));
   grid on;
   xlabel('Time(sec)');
   title('Acceleration');
   if(iu==1)
       ylabel('Accel(in/sec^2)');
   else
       ylabel('Accel(m/sec^2)');    
   end
%
   XF=[tt' x];
   VF=[tt' v];
   AF=[tt' acc]; 
%
   clear FF;
   clear ff;
%
   ff(:,1)=x(:,2)-x(:,1);
   for(i=2:ne)
       ff(:,i)=(x(:,i+1)-x(:,i-1))/2;
   end
   ff=E*Area*ff/h;
%
   FF=[tt' ff];
%
max(ff(:,1))
end