%   three_dof.m   ver 1.0  March 22, 2005
%
%   by Tom Irvine   Email: tomirvine@aol.com
%
disp(' ');
disp(' This program finds the eigenvalues and eigenvectors for a ');
disp(' three-degree-of-freedom system.');
disp(' ');
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=3;
%
disp(' Assume symmetric mass and stiffness matrices. ');
%
disp(' ');
disp(' Enter m11 ');
m(1,1)=input(' ');
%
disp(' Enter m12 ');
m(1,2)=input(' ');
%	
disp(' Enter m13 ');	
m(1,3)=input(' ');
%
disp(' Enter m22 ');	
m(2,2)=input(' ');
%
disp(' Enter m23 ');	
m(2,3)=input(' ');
%
disp(' Enter m33 ');	
m(3,3)=input(' ');
%
disp(' ');
disp(' Enter k11 ');
k(1,1)=input(' ');
%
disp(' Enter k12 ');
k(1,2)=input(' ');
%	
disp(' Enter k13 ');	
k(1,3)=input(' ');
%
disp(' Enter k22 ');	
k(2,2)=input(' ');
%
disp(' Enter k23 ');	
k(2,3)=input(' ');
%
disp(' Enter k33 ');	
k(3,3)=input(' ');
%
% Symmetry
%
m(2,1)=m(1,2);
m(3,1)=m(1,3);
m(3,2)=m(2,3);
k(2,1)=k(1,2);
k(3,1)=k(1,3);
k(3,2)=k(2,3);
%
disp(' ');
disp(' The mass matrix is');
m
disp(' ');
disp(' The stiffness matrix is');
k
%
%  Calculate eigenvalues and eigenvectors
%
[ModeShapes,Eigenvalues]=eig(k,m);
%
omega = sqrt(Eigenvalues);
fn=omega/tpi;
%
disp(' ');
disp('  Natural Frequencies = ');
out1=sprintf('\n  %8.4g Hz',fn(1,1));
out2=sprintf('  %8.4g Hz',fn(2,2));
out3=sprintf('  %8.4g Hz',fn(3,3));
disp(out1);
disp(out2);
disp(out3);
%
MST=ModeShapes';
temp=zeros(3,3);
temp=m*ModeShapes;
QTMQ=MST*temp;
%   
scale1=1./sqrt(QTMQ(1,1));
scale2=1./sqrt(QTMQ(2,2));
scale3=1./sqrt(QTMQ(3,3));
%
ModeShapes(:,1) = ModeShapes(:,1)*scale1;  
ModeShapes(:,2) = ModeShapes(:,2)*scale2; 
ModeShapes(:,3) = ModeShapes(:,3)*scale3; 
%
MST=ModeShapes';
%
disp(' ');
disp('  Modes Shapes (column format) =');
disp(' ');
out1=sprintf(' %12.4g   %12.4g   %12.4g ',ModeShapes(1,1),ModeShapes(1,2),ModeShapes(1,3));
out2=sprintf(' %12.4g   %12.4g   %12.4g ',ModeShapes(2,1),ModeShapes(2,2),ModeShapes(2,3));
out3=sprintf(' %12.4g   %12.4g   %12.4g ',ModeShapes(3,1),ModeShapes(3,2),ModeShapes(3,3));
disp(out1);
disp(out2);
disp(out3);