disp(' ');
disp(' mdof_rk4_free.m  ver 1.5  December 21, 2011 ');
disp(' by Tom Irvine ');
disp(' ');
disp(' This program calculates the response of an MDOF ');
disp(' system to initial conditions using the RK4 method.');
disp(' ');
%
clear DI;
clear VI;
clear w;
clear t;
clear M;
clear C;
clear K;
clear aaa;
clear x;
clear y;
clear accel;
clear length;
clear dampV;
clear ccc;
%
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(' damping unit = lbf sec/in ');
else
    disp(' damping unit = N sec/m ');
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
ndof=max(size(M));
%
if(iu==1 && imu==1)
  M=M/386;
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select damping input method ')
disp('  1=damping coefficient matrix ')
disp('  2=uniform modal damping ')
disp('  3=modal damping vector ')
%
idamp=input(' ');
%
if(idamp==1)
    disp(' ');
    disp(' Damping Matrix ');
%
    if(file_choice==1)
        C = input(' Enter the matrix name:  ');
    end
    if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        C = xlsread(xfile);
%         
    end
end
if(idamp==2)
   disp(' Enter uniform damping ratio ');
   damp=input(' ');  
end    
if(idamp==3)
   disp(' Enter damping ratio vector ');
%
    if(file_choice==1)
        dampV = input('  Enter the damping vector name:  ');
    end
    if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        dampV = 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
% 
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(K,M,1);
%
if(idamp==2 || idamp==3)
    ccc=zeros(ndof,ndof);
%
    for i=1:ndof
        if(idamp==2)
            ccc(i,i)=2*damp*(omegan(i));
        else
            ccc(i,i)=2*dampV(i)*(omegan(i));        
        end    
    end
%   
    s1=(MST*M);
    s2=(ccc*s1);
    s3=(ModeShapes*s2);
    C=(M*s3);
    C
end
%
disp(' ');
disp(' Enter the initial velocity vector ');
VI=input(' ');
disp(' ');
disp(' Enter the initial displacement vector ');
DI=input(' ');
disp(' ');
%
Tmin=1/(max(fn));
Tmax=1/(min(fn));
%
dur=20*Tmax;
dt=Tmin/50.;
%
NT=round(dur/dt);
%
t=linspace(0,dur,NT);
%
ndof=max(size(M));
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
x=zeros(NT,ndof);
y=zeros(NT,ndof);
accel=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);
%
x(1,:)=DI;
y(1,:)=VI;
%
for(i=1:(NT-1))
% 
	    h=t(i+1)-t(i);
%
	    hd2=h/2;
%
        xx=x(i,:)';
        yy=y(i,:)';
%
        X1= xx;
        Y1= yy;
        F1= M\(-C*Y1-K*X1);
%
        X2= X1+Y1*hd2;
        Y2= Y1+F1*hd2;
        F2= M\(-C*Y2-K*X2);
%
        X3= X1+Y2*hd2;
        Y3= Y1+F2*hd2;
        F3= M\(-C*Y3-K*X3);
%
        X4= X1+Y3*h;
        Y4= Y1+F3*h;
        F4= M\(-C*Y4-K*X4);
%
        x(i+1,:)= (X1+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 ))';   % displacement
        y(i+1,:)= (Y1+(h/6.)*( F1 + 2*F2 + 2*F3 + F4 ))';   % velocity
        accel(i+1,:)= M\(-C*y(i+1,:)'-K*x(i+1,:)');         % acceleration
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
num=ndof;
%
mdof_plot(t,x,y,accel,num,iu);