disp(' ');
disp(' mdof.m   ver 1.2  February 10, 2010 ');
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 ');
%
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 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,omega,ModeShapes,MST]=Generalized_Eigen(k,m,1);
%
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
%
out1=sprintf('\n Total Modal Mass = %12.4f ',tm);
disp(out1)
%