disp(' ');
disp(' mdof_free.m   ver 1.3   December 21, 2011 ');
disp(' ');
disp(' by Tom Irvine   Email: tomirvine@aol.com ');
disp(' ');
disp(' This program solves the following equation of motion: ');
disp('    M (d^2x/dt^2) + K x = 0 ');
disp('   ');
disp(' where the system is subjected to initial conditions. ');
%
clear acc;
clear dis;
clear vel;
clear k;
clear m;
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 a;
clear x;
clear d;
clear v;
clear na;
clear nv;
clear nd;
clear Q;
clear mass;
clear stiffness;
clear max;
%
tpi=2.*pi;
%
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
%
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
%
m=m/mass_scale;
%
mass=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
stiffness=k;
%
size(m)
size(k);
%
num=max(size(m));
%
disp(' ');
disp(' The mass matrix is');
m
disp(' ');
disp(' The stiffness matrix is');
k
%
%  Calculate eigenvalues and eigenvectors
%
[fn,omegan,ModeShapes,MST]=Generalized_Eigen(k,m,1);
clear Q;
Q=ModeShapes;
%
clear r;
clear part;
clear ModalMass;
%
r=ones(num,1);
%
part = MST*m*r;
%
disp(' Participation Factors  ');
part
%
ModalMass = (part.*part)/1.;
%
if(iu==1)
    ModalMass=ModalMass*386;
end
disp(' Effective Modal Mass ');
ModalMass
%
tm=0;
for(i=1:num)
   tm=tm+ModalMass(i,1);
end
nor=num;
%
out1=sprintf('\n Total Modal Mass = %12.4f ',tm);
disp(out1)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear VI;
clear DI;
clear T;
clear nt;
clear nDI;
clear nVI;
clear damp;
%
disp(' ');
disp(' Select damping input method ');
disp('  1=uniform modal damping ');
disp('  2=modal damping vector ');
idm=input(' ');
%
if(idm==1)
% 
    disp(' ');
    disp(' Enter modal damping ratio ');
    udamp=input(' ');
    damp=ones(num,1);
    damp=damp*udamp;
% 
else
    if(file_choice==1)
        disp(' ');
        damp = input(' Enter the modal damping vector  ');
    end
    if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        damp = xlsread(xfile);
%         
    end
end
%
if(file_choice==1)
        disp(' ');    
        VI = input(' Enter the initial velocity vector:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        VI = xlsread(xfile);
%         
end
%
if(file_choice==1)
        disp(' ');    
        DI = input(' Enter the initial displacement vector:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        DI = xlsread(xfile);
%         
end
%
sz=size(DI);
if(sz(2)>sz(1))
    DI=DI';
end
%
sz=size(VI);
if(sz(2)>sz(1))
    VI=VI';
end
%
nDI=MST*(m*DI);
nVI=MST*(m*VI);
%
nic=num;
%
disp(' ');
disp(' Enter the duration (sec) ');
T=input(' ');
%
disp(' ');
minSR=10.*(max(fn(nic)));
%
out1=sprintf(' The suggest sample rate > %8.4g (samples/sec)',minSR);
disp(out1);
%
disp(' ');
disp(' Enter the sample rate (samples/sec) ');
sr=input(' ');
%
dt=1/sr;
%
nt=T/dt;
%
clear A;
clear B;
clear n;
clear nv;
clear na;
clear x;
clear v;
clear a;
clear t;
clear tt;
clear omegad;
%
for(i=1:num)
   omegad(i)=omegan(i)*sqrt(1-(damp(i))^2);
end
%
sz=size(omegan);
if(sz(2)>sz(1))
    omegan=omegan';
end
%
sz=size(omegad);
if(sz(2)>sz(1))
    omegad=omegad';
end
%
sz=size(damp);
if(sz(2)>sz(1))
    damp=damp';
end
%
progressbar % Create figure and set starting time
%
clear A;
clear B;
clear om2;
clear b1;
clear b2;
clear b3;
%
om2=omegan.*omegan;
%
A=nDI; 
b1=omegan.*A;
b2=damp.*b1;
b3=b2+nVI;
%
for(j=1:num)
   B(j)=b3(j)/omegad(j);
end
%
sz=size(A);
if(sz(2)>sz(1))
    A=A';
end
%
sz=size(B);
if(sz(2)>sz(1))
    B=B';
end
%
for(i=1:nt)
    progressbar(i/nt) % Update figure        
    t=dt*(i-1);
    tt(i)=t;
%
    arg=omegad*t;
%
    sz=size(arg);
    if(sz(2)>sz(1))
       arg=arg';
    end
%
    cos_arg=cos(arg);
    sin_arg=sin(arg);
    eee=exp(-damp.*omegan*t); 
%
    sz=size(cos_arg);
    if(sz(2)>sz(1))
       cos_arg=cos_arg';
    end    
%
    sz=size(sin_arg);
    if(sz(2)>sz(1))
       sin_arg=sin_arg';
    end  
%
    sz=size(eee);
    if(sz(2)>sz(1))
       eee=eee';
    end      
%
    A_cos_arg=A.*cos_arg;
    A_sin_arg=A.*sin_arg;
    B_cos_arg=B.*cos_arg;    
    B_sin_arg=B.*sin_arg;   
%
    term1=A_cos_arg+B_sin_arg;
    n(1:num)=eee.*term1;
%
    sz=size(n);
    if(sz(2)>sz(1))
       n=n';
    end      
%
    term2=-A_sin_arg+B_cos_arg;
%
    nv=-damp.*(omegan.*n) + omegad.*( eee.*term2 );
%
    na=-2*damp.*(omegan.*nv) - om2.*n;    
%
    x(i,:)=Q(:,1:nic)*n(1:nic);
    v(i,:)=Q(:,1:nic)*nv(1:nic);
    a(i,:)=Q(:,1:nic)*na(1:nic);  
%
end
%
%
[colororder,color_rows]=line_colors();
%
figure(1);
hold('all');
%
j=1;
for i = 1:num
  if(j>color_rows)
      j=j-color_rows;
  end
  plot(tt',x(:,i),'Color', colororder(j,:));
  j=j+1;
end
hold off;
grid on;
xlabel('Time(sec)');
title('Displacement');
if(iu==1)
    ylabel('Disp(in)');
else
    ylabel('Disp(m)');    
end
plot_legend(num);
%
figure(2);
hold('all');
%
j=1;
for i = 1:num
  if(j>color_rows)
      j=j-color_rows;
  end
  plot(tt',v(:,i),'Color', colororder(j,:));
  j=j+1;
end
hold off;
grid on;
xlabel('Time(sec)');
title('Velocity');
if(iu==1)
    ylabel('Vel(in/sec)');
else
    ylabel('Vel(m/sec)');    
end
plot_legend(num);
%
figure(3);
%%
if(iu==1)
    a=a/386;
end
%%
hold('all');
%
j=1;
for i = 1:num
  if(j>color_rows)
      j=j-color_rows;
  end
  plot(tt',a(:,i),'Color', colororder(j,:));
  j=j+1;
end
hold off;
%%
xlabel('Time(sec)');
title('Acceleration');
if(iu==1)
    ylabel('Accel(G)');
else
    ylabel('Accel(m/sec^2)');    
end
grid on;
plot_legend(num);
%
XF=[tt' x];
VF=[tt' v];
AF=[tt' a];
%