disp(' ');
disp(' mdof_modal_arbit_force_rk4.m  ver 1.3  December 21, 2011 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the response of an MDOF ');
disp(' system to arbitrary force excitation via the RK4 method.');
disp('  ');
disp(' The system is decoupled using normal modes as an  ');
disp(' intermediate step. ');
disp('  ');
disp(' It is currently set-up for a two-degree-of-freedom ');
disp(' system, but can be readily modified for additional degrees. ');
disp(' ');
%
clear damp
clear DI;
clear VI;
clear w;
clear t;
clear M;
clear C;
clear K;
clear accel;
clear length;
clear tt;
clear ff;
clear FI;
clear tint;
clear fint;
clear a;
clear x;
clear y;
%
disp(' ');
disp(' Enter the units system ');
disp(' 1=English  2=metric ');
iu=input(' ');
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
if(iu==1)
     disp(' Select input mass unit ');
     disp('  1=lbm  2=lbf sec^2/in  ');
     imu=input(' ');

else
    disp(' mass unit = kg ');
end
%
disp(' ');
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
%
if(iu==1 && imu==1)
  M=M/386;
end
%
sz=size(M);
ndof=sz(1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select damping input method ');
disp('   1=uniform damping ratio ');
disp('   2=damping ratio vector ');
%
idm=input(' ');
if(idm==1)
    disp(' Enter damping ratio ');
    udamp=input(' ');
    for(i=1:ndof)
        damp(i)=udamp;
    end    
else
%
    if(file_choice==1)
        damp = input(' Enter the damping vector name:  ');
    end
    if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        damp = xlsread(xfile);
%         
    end
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
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
%
disp(' ');
disp(' Enter the initial velocity vector ');
VI=input(' ');
disp(' ');
disp(' Enter the initial displacement vector ');
DI=input(' ');
disp(' ');
% 
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(K,M,1);
%
Tmin=1/(max(fn));
Tmax=1/(min(fn));
%
disp(' ');
disp(' Enter duration(sec)');
dur=input(' ');
%
srr=20*max(fn);
disp(' ');
disp(' Enter sample rate (samples/sec)');
out1=sprintf(' (recommend %8.4g)',srr);
disp(out1);
sr=input(' ');
%
dt=1/sr;
%
NT=round(dur/dt);
%
t=linspace(0,dur,NT);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[FFI,force_dof,nff]=ODE_force_input(iu,ndof,NT,t,file_choice);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nx=zeros(NT,ndof);
ny=zeros(NT,ndof);
%
X1=zeros(ndof,1);
X2=zeros(ndof,1);
X3=zeros(ndof,1);
X4=zeros(ndof,1);
%
Y1=zeros(ndof,1);
Y2=zeros(ndof,1);
Y3=zeros(ndof,1);
Y4=zeros(ndof,1);
%
nx(1,:)=(MST*(M*DI))';
ny(1,:)=(MST*(M*VI))';
%
for(j=1:ndof)
    cc(j)=2*damp(j)*omegan(j);
    kk(j)=(omegan(j))^2;
end    
%
for(i=1:(NT-1))
% 
	    h=t(i+1)-t(i);
%
	    hd2=h/2;
%
        FA=zeros(ndof,1);
        FB=zeros(ndof,1);
%
        for(j=1:ndof)
%
            j_index=force_dof(j);
%
            if(j_index~=-999)
                FA(j)=FFI(i,j_index);
                FB(j)=FFI(i+1,j_index);
            end
%
        end
%
        nFA=MST*FA;
        nFB=MST*FB;   
%
        FF1= nFA;
        FF2= (nFB + nFA)/2.;
        FF3= FF2;
        FF4= nFB;
%
        for(j=1:ndof)
            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     
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear x;
clear y;
clear a;
%
x=zeros(NT,ndof);
y=zeros(NT,ndof);
a=zeros(NT,ndof);
%
for(i=1:NT)
    x(i,:)=(ModeShapes*(nx(i,:))')';
    y(i,:)=(ModeShapes*(ny(i,:))')';
    a(i,:)=(ModeShapes*(na(i,:))')';    
end
%
num=ndof;
%
mdof_plot(t,x,y,a,num,iu);