disp(' ');
disp(' rod_modal_transient_FEA.m  ver 1.0   June 25, 2012 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the longitudinal vibration, ');
disp(' modal transient response to of a fixed-free rod an  '); 
disp(' arbitrary base input via the finite element method. ');
%
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;
%
clear damp;
clear ea;
clear f;
clear ngw;
%
clear ttend;
clear ttstart;
clear length;
%
clear inv_KH;
clear inv_KH_diag;
%
clear dT;
clear vT;
clear accT;
%
clear ZdT;
clear ZvT;
clear ZaccT;
%
clear dvelox;
clear ddisp;
%
clear forig;
%
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(' ');
%
ibc=1;
%
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
%
%  S2
%  M2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% disp(' ');
%% disp('  mdof_modal_enforced_acceleration_ri.m  ver 1.1  February 9, 2012 ');
%% disp('  by Tom Irvine ');
%
%% disp(' ');
%% disp('  Response of a multi-degree-of-freedom system to enforced ');
%% disp('  acceleration at selected point mass nodes via a ramp invariant ');
%% disp('  digital recursive filtering relationship. ');
%% disp(' ');
%
close all hidden;
%
clear Mdd;  % driven
clear Mdf;
clear Mfd;
clear Mff;  % free
%
clear Mdw;
clear Mwd;
clear Mww;
%
clear Kdw;
clear Kwd;
clear Kww;
%
clear MP;
clear MT;
%
clear Kdd;  
clear Kdf;
clear Kfd;
clear Kff;
%
clear invKff;
%
%
clear KP;
clear KT;
%
clear n;
clear na;
clear nv;
clear nd;
clear nx;
%
clear m;
clear k;
%
clear TT;
clear T1;
clear T2;
%
clear I;
clear Idd;
clear Iff;
%
clear damp;
clear ea;
clear f;
clear ngw;
%
clear ttend;
clear ttstart;
clear length;
%
clear inv_KH;
clear inv_KH_diag;
%
clear dT;
clear vT;
clear accT;
%
clear ZdT;
clear ZvT;
clear ZaccT;
%
clear dvelox;
clear ddisp;
%
clear forig;
%
%% disp(' ');
%% disp(' Enter the units system ');
%% disp(' 1=English  2=metric ');
%% iu=input(' ');
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
mass_scale=1;
%
if(iu==1)
     disp(' Select input mass unit ');
     disp('  1=lbm  2=lbf sec^2/in  ');
%%      imu=input(' ');
     imu=2;
     if(imu==1)
         mass_scale=386;
     end
else
    disp(' mass unit = kg ');
end
%
if(iu==1)
    disp(' stiffness unit = lbf/in ');
else
    disp(' stiffness unit = N/m ');
end
%
disp(' ');
%% disp(' Select file input method ');
%% disp('   1=file preloaded into Matlab ');
%% disp('   2=Excel file ');
%% file_choice = input('');
%
file_choice=1;
%
disp(' ');
disp(' Mass Matrix ');
%
%% if(file_choice==1)
%%         m = input(' Enter the matrix name:  ');
%% end
%% if(file_choice==2)
%%         [filename, pathname] = uigetfile('*.*');
%%         xfile = fullfile(pathname, filename);
%
%%         m = xlsread(xfile);
%
%% end
%
m=M2;
%
m=m/mass_scale;
%
mass=m;
num=max(size(m));
%
%% disp(' ');
%% disp(' Stiffness Matrix ');
%
%% if(file_choice==1)
%%         k = input(' Enter the matrix name:  ');
%% end
%% if(file_choice==2)
%%         [filename, pathname] = uigetfile('*.*');
%%         xfile = fullfile(pathname, filename);
%
%%         k = xlsread(xfile);
%
%% end
%% stiff=k;
%
stiff=S2;
%
disp(' Input Matrices ');
mass
stiff
%
[fn,omega,ModeShapes,MST]=Generalized_Eigen(stiff,mass,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' Select modal damping input method ');
disp('  1=uniform damping for all modes ');
disp('  2=damping vector ');
idm=input(' ');
%
disp(' ');
if(idm==1)
   disp(' Enter damping ratio ');
   dampr=input(' ');
   damp(1:num)=dampr;
else
   disp(' Modal Damping Vector ');
%
   if(file_choice==1)
        damp = input(' Enter the vector:  ');
   end
   if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        damp = xlsread(xfile);
%         
   end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sz=size(mass);
ndof=(sz(1));
num=ndof;
%
out1=sprintf('\n number of dofs =%d \n',num);
disp(out1);
%
disp(' ');
disp(' Enter the duration (sec) ');
dur=input(' ');
disp(' ');
sug=20*(max(omega)/(2*pi));
%
disp(' Enter the sample rate (samples/sec)');
out1=sprintf(' (suggest > %8.4g ) ',sug);
disp(out1);
sr=input(' ');
dt=1/sr;
nt=1+round(dur/dt);
%
t =  linspace(0,dur,nt);
nt=length(t);
t=t';
%
[FFI,accel_dof,nff,nem,nfree,enforced_dof,free_dof,ngw]=...
    ODE_acceleration_input_rod(iu,num,nt,t,file_choice);
%
forig=zeros(nt,nem);
%
k=1;
for i=1:num
    nfi=accel_dof(i);
    if(nfi>0)
        forig(:,k)=FFI(:,nfi);
        k=k+1;
    end    
end
%
clear FFI;
FFI=forig;
%
for(i=1:nem)
    A=FFI(:,i);
    dv=cumtrapz(A);
    dvelox(:,i)=dv*dt;
end
%
for(i=1:nem)
    A=dvelox(:,i);
    dd=cumtrapz(A);
    ddisp(:,i)=dd*dt;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Partition Mass & Stiffness Matrices
%
[fn,omega,ModeShapes,MST,Mwd,Kwd,Mww,Kww,TT,T1,T2,KT,MT]=...
    partition_matrices(mass,stiff,num,nfree,nem,free_dof,enforced_dof,'accel');
%
%
omegan=omega;
%
clear r;
clear part;
clear ModalMass;
%
r=ones(nfree,1);
%
part = MST*Mww*r;
%
disp(' Participation Factors  ');
part
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Modal Transient
%
clear omegad;
for i=1:nfree
    omegad(i)=omegan(i)*sqrt(1-damp(i)^2);
end
%
clear Q;
Q=ModeShapes;
QT=Q';
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  displacement
%
%  d
%  v
%  acc
%
%
%%
%%
kk=zeros(nfree,1);
cc=zeros(nfree,1);
%
for i=1:nfree
    kk(i)=(omegan(i))^2.;
    cc(i)=2*damp(i)*omegan(i);
end
%
t=linspace(0,dur,nt);
%
clear d;
clear v;
clear acc;
clear nFA;
%
d=zeros(nt,nfree);
v=zeros(nt,nfree);
acc=zeros(nt,nfree);
%
nx=zeros(nt,nfree);
ny=zeros(nt,nfree);
na=zeros(nt,nfree);
%
nFA=zeros(nt,nfree);
%
for i=1:nt
%    
        nFA(i,:)=-QT*Mwd*FFI(i,:);  
%
end    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Calculate Filter Coefficients
%
disp(' ');
disp(' Calculate Filter Coefficients ');
% 
[a1,a2,df1,df2,df3,vf1,vf2,vf3,af1,af2,af3]=ramp_invariant_filter_coefficients(nfree,omegan,damp,dt);
%
%  Numerical Engine
%
disp(' ')
disp(' Calculating response...');
%
nx=zeros(nt,nfree);
nv=zeros(nt,nfree);
na=zeros(nt,nfree);
%
for j=1:nfree
%
%  displacement
%
    d_forward=[   df1(j),  df2(j), df3(j) ];
    d_back   =[     1, -a1(j), -a2(j) ];
    d_resp=filter(d_forward,d_back,nFA(:,j));
%    
%  velocity
%
    v_forward=[   vf1(j),  vf2(j), vf3(j) ];
    v_back   =[     1, -a1(j), -a2(j) ];
    v_resp=filter(v_forward,v_back,nFA(:,j));
%    
%  acceleration
%   
    a_forward=[   af1(j),  af2(j), af3(j) ];
    a_back   =[     1, -a1(j), -a2(j) ]; 
    a_resp=filter(a_forward,a_back,nFA(:,j));
%
    nx(:,j)=d_resp;  % displacement
    nv(:,j)=v_resp;  % velocity
    na(:,j)=a_resp;  % acceleration  
%
end
%
%
for i=1:nt
    d(i,:)=(Q*(nx(i,:))')';
    v(i,:)=(Q*(nv(i,:))')';
    acc(i,:)=(Q*(na(i,:))')';    
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear d_save;
clear v_save;
clear a_save;
d_save=d;
v_save=v;
a_save=acc;
%
%    Transformation back to xd xf
%
%%    dT=TT*d';
%%    vT=TT*v';
%%    accT=TT*acc';
%
dT=(T1*ddisp' + T2*d')';
%
vT=(T1*dvelox' + T2*v')';
%
accT=(T1*FFI' + T2*acc')';
%
%  Put in original order
%
ZdT=[ ddisp  dT ];
ZvT=[ dvelox  vT ];
ZaccT=[ forig  accT ];
%
clear d;   
clear v;
clear acc;
d=zeros(nt,num);
v=zeros(nt,num);
acc=zeros(nt,num);
%
for i=1:num
   for j=1:nt
        d(j,ngw(i))=  (ZdT(j,i));        
        v(j,ngw(i))=  (ZvT(j,i));     
      acc(j,ngw(i))=(ZaccT(j,i));
   end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
mdof_plot_rod(t,d,v,acc,num,iu)
%