%   two_dof_force_sine.m   ver 1.0  June 1, 2010
%
%   by Tom Irvine   Email: tomirvine@aol.com
%
disp(' ');
disp(' This program finds the eigenvalues and eigenvectors for a ');
disp(' two-degree-of-freedom system.');
disp(' The equation of motion is:   M (d^2x/dt^2) + K x = 0 ');
%
clear k;
clear m;
clear Eigenvalues;
clear ModeShapes;
clear omega;
clear fn;
clear MST;
clear ModalMass;
clear part;
clear QTMQ;
clear r;
%
tpi=2.*pi;
%
n=2;
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
disp(' ');
disp(' Enter m11 ');
m(1,1)=input(' ');
%
disp(' Enter m12 ');
m(1,2)=input(' ');
%	
disp(' Enter m22 ');	
m(2,2)=input(' ');
%
m(2,1)=m(1,2);
%
disp(' ');
disp(' Enter k11 ');
k(1,1)=input(' ');
%
disp(' Enter k12 ');
k(1,2)=input(' ');
%	
disp(' Enter k22 ');	
k(2,2)=input(' ');
%
k(2,1)=k(1,2);
%
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);
%
MST=ModeShapes';
%
r(1)=1;
r(2)=1;
%
part = MST*m*r';
%
disp(' ')
disp(' Particpation Factors = ');
disp(' ');
out1=sprintf('  %8.4g   ', part(1));
out2=sprintf('  %8.4g   ', part(2));
disp(out1);
disp(out2);
%
ModalMass = (part.*part)/1.;
%
disp(' ');
disp(' Effective Modal Mass =');
disp(' ');
out1=sprintf('  %8.4g   ', ModalMass(1));
out2=sprintf('  %8.4g   ', ModalMass(2));
disp(out1);
disp(out2);
%
tm=ModalMass(1)+ModalMass(2);
%
out1=sprintf('\n Total Modal Mass = %12.4f ',tm);
disp(out1)
%
disp(' ');
disp(' Enter the damping ratio for mode 1 ');
damp(1)=input(' ');
%
disp(' ');
disp(' Enter the damping ratio for mode 2 ');
damp(2)=input(' ');
%
disp(' ');
disp(' Enter the amplitude (lbf) for force 1 ');
B(1)=input(' ');
%
disp(' ');
disp(' Enter the frequency (Hz) for force 1 ');
f1=input(' ');
alpha=tpi*f1;
%
disp(' ');
disp(' Enter the amplitude (lbf) for force 2 ');
B(2)=input(' ');
%
disp(' ');
disp(' Enter the frequency (Hz) for force 2 ');
f2=input(' ');
beta*tpi*f2;
%
disp(' ');
disp(' Enter the duration (sec) ');
dur=input(' ');
%
disp(' ');
disp(' Enter the samppe rate (samples/sec) ');
sr=input(' ');
%
dt=1/sr;
nt=round(dur/dt);
%
for(j=1:2)
   den1(j)=omegan(j)^2-alpha^2;
   den2(j)=omegan(j)^2-beta^2;
end
%
for(i=0:nt)
   t=dt*i;
   for(j=1:2)
       n(j)=(ModeShapes(xx)*(B(j)*alpha)*( )/den1;
   end
end    
%