disp(' ');
disp('  mdof_modal_enforced_acceleration_rk4.m  ver 1.0  December 27, 2011 ');
disp('  by Tom Irvine ');
%
disp(' ');
disp('  Response of a multi-degree-of-freedom system to enforced ');
disp('  acceleration at selected point mass nodes via the Runge-Kutta ');
disp('  fourth order method. ');
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 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(' ');
     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('');
%
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=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;
%
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(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;
%
d=zeros(nt,nfree);
v=zeros(nt,nfree);
acc=zeros(nt,nfree);
%
nx=zeros(nt,nfree);
ny=zeros(nt,nfree);
na=zeros(nt,nfree);
%%
for i=1:(nt-1)
%    
  	    h=t(i+1)-t(i);
%
	    hd2=h/2;
%
        nFA=-QT*Mwd*FFI(i,:);
        nFB=-QT*Mwd*FFI(i+1,:);   
%
        FF1= nFA;
        FF2= (nFB + nFA)/2.;
        FF3= FF2;
        FF4= nFB;
%
        for(j=1:nfree)
            xx=nx(i,j)';
            yy=ny(i,j)';
%
            X1= xx;
            Y1= yy;
            F1= (FF1(j)-cc(j)*Y1-kk(j)*X1);
%
            X2= X1+Y1*hd2;
            Y2= Y1+F1*hd2;
            F2= (FF2(j)-cc(j)*Y2-kk(j)*X2);
%
            X3= X1+Y2*hd2;
            Y3= Y1+F2*hd2;
            F3= (FF3(j)-cc(j)*Y3-kk(j)*X3);
%
            X4= X1+Y3*h;
            Y4= Y1+F3*h;
            F4= (FF4(j)-cc(j)*Y4-kk(j)*X4);
%
            nx(i+1,j)= (X1+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 ))';    % displacement
            ny(i+1,j)= (Y1+(h/6.)*( F1 + 2*F2 + 2*F3 + F4 ))';    % velocity
            na(i+1,j)= FF4(j) -cc(j)*ny(i+1,j) -kk(j)*nx(i+1,j); % acceleration  
        end 
%
        d(i+1,:)=(Q*nx(i+1,:));
        v(i+1,:)=(Q*ny(i+1,:));
        acc(i+1,:)=(Q*na(i+1,:));   
%
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(t,d,v,acc,num,iu)
%