% Example 1
% Two-Degree of Freedom
clear all;
clc;

k=1;
m=1;

m1=m;
m2=2*m;
k1=k;
k2=k1;
k3=2*k;

% Using Classnotes Formulation
fprintf('Using Class-Notes Formulation...\n')
term1=((k1+k2)*m2+(k2+k3)*m1)/m1/m2;
term2=term1;
term3=((k1+k2)*(k2+k3)-k2^2)/m1/m2;

% Eigen-values
r1=1/2*term1-1/2*sqrt(term2^2-4*term3);
r2=1/2*term1+1/2*sqrt(term2^2-4*term3);
fprintf('\n Eigen-Value 1 = %g \n', r1)
fprintf('\n Eigen-Value 2 = %g \n', r2)

% Eigen-Vectors
X2_X1_mode1=(-m1*r1+(k1+k2))/k2;
X2_X1_mode2=(-m1*r2+(k1+k2))/k2;

% Normal-Modes or Modal Vectors
fprintf('\n Eigen-Vectors \n')
% Modal Vector for Mode 1
%Assume X1=1
X1=1;
mode1=[X1 ; X2_X1_mode1*X1]
% Modal Vector for Mode 2
%Assume X1=1
X1=1;
mode2=[X1 ; X2_X1_mode2*X1]

fprintf('================================================================ \n')
% Using Matlab and Matrix Notation
fprintf('Using Matlab and Matrix Formulation...\n')
M=[m1 0;0 m2];
K=[k1+k2 -k2;-k2 k2+k3];
% function eig in Matlab to eigen-Values/Vectors
[Eigen_vectors,Eigen_values]=eig(K,M);

Eg_values=diag(Eigen_values);
fprintf('\n Eigen-Value 1 = %g \n', Eg_values(1))
fprintf('\n Eigen-Value 2 = %g \n', Eg_values(2))

fprintf('\n Eigen-Vectors \n')
Eigen_vectors(:,1)
Eigen_vectors(:,2)


%Normalize to 1
fprintf('\n Normalization of Eigen-Vectors to 1... \n')
for i=1:2
    work1=max(Eigen_vectors(:,i));
    work2=min(Eigen_vectors(:,i));
    if abs(work1)>=abs(work2);
    New_Eig_Vector(:,i)=Eigen_vectors(:,i)/work1;
    else abs(work1)<abs(work2);
    New_Eig_Vector(:,i)=Eigen_vectors(:,i)/work2;  
    end
end

fprintf('\n Normalized Eigen-Vectors \n')
New_Eig_Vector(:,1)
New_Eig_Vector(:,2)


fprintf('================================================================ \n')

fprintf('\n Check Orthogonality of Mode Vectors \n')
fprintf('Using Class-Notes Formulation...\n')
Phi=[mode1 mode2];
Phi'*M*Phi
Phi'*K*Phi

fprintf('Using Matlab and Matrix Formulation...\n')
Phi=New_Eig_Vector;
Phi'*M*Phi
Phi'*K*Phi

fprintf('\n Note:')
fprintf('\n The orthogonal property of the mode shapes is very important since it')
fprintf('\n decouples the equation of motion.  Decoupling means that the first')
fprintf('\n equation of motion will be in function of x1 only and the second equation of motion')
fprintf('\n in function of x2 only.  In other words the 2 degree  of freedom systen is')
fprintf('\n reduced to two 1 single degree of freedom systems')

