% Frank Fisher
% ME345
% Eigenvalue problem... solution via Matlab
% This is for the two spring and two mass system

%           --------                ------------
%    k1     |       |       k2      |           |
%--/\/\/\---|   m1  |------/\/\/\---|    m2     |
%           |       |               |           |
%           --------                -------------

clear
clc
close all

% INPUT system parameters (in class, used k1=3, k2=2, m1=1, m2=1)
k1=3;
k2=2;
m1=1;
m2=1;

% Input initial conditions (assumed initial velocity and accel of masses)
vm1_t0=0;           % init velocity of m1
vm2_t0=0;           % init velocity of m2
am1_t0=3.45;        % init accel of m1
am2_t0=.775;        % init accel of m2

% If want to set specific values for the constants (to get normal mode
% behavior) set this

normal_mode=0;      % if this is set equal to 1, use THESE constants in final solution

%%%%%%%% DO NOT CHANGE ANYTHING BELOW THIS LINE %%%%%%%% 

A = [-k1/m1 - k2/m1, k2/m1;
    k2/m2, -k2/m2];             %this is the A matrix

[e_vectors, e_values]=eig(A);   %Matlab built-in command for eigenvalues and eigenvectors
                                %e_vectors - columns give the eigenvectors
                                %e_values - diagonals given the eigenvalues
                                
evect1=e_vectors(:,1)          %this is one of the eigenvectors
evect2=e_vectors(:,2)          %this is the other eigenvector
eval1=-e_values(1,1)           %one of the eigenvalues (note: negative here)
eval2=-e_values(2,2)           %one of the eigenvalues (note: negative here)
                                
%Following the notation of class, once we have these values we can write
%the velocities of the masses as a function of time (and unknown constants
%based on the initial conditions).

%To find the accelerations, also need to take the derivatives of these
%velocities (See class notes)

%To solve for unknowns a1, a2, b1, b2, use boundary conditions. Will solve
%in matrix form. Note: assuming that all BCs are at t = 0!
temp_matrix1=[evect1(1),evect2(1);evect1(2),evect2(2)];
a1_a2=inv(temp_matrix1)*[vm1_t0;vm2_t0];
a1=a1_a2(1);
a2=a1_a2(2);

temp_matrix2=[evect1(1)*sqrt(eval1),evect2(1)*sqrt(eval2);evect1(2)*sqrt(eval1),evect2(2)*sqrt(eval2)];
b1_b2=inv(temp_matrix2)*[am1_t0;am2_t0];
b1=b1_b2(1);
b2=b1_b2(2);


%%% BELOW THIS LINE MIGHT WANT TO CHANGE!!!!

if normal_mode==1
    %a1=0; b1=0; a2=1; b2=1;     %first mode
    a1=1; b1=1; a2=0; b2=0;     %second mode
end

% Now plug these into our equations and can plot versus time (t)
t=0:.1:10;
vm1=evect1(1)*(a1*cos(sqrt(eval1)*t)+b1*sin(sqrt(eval1)*t)) + evect2(1)*(a2*cos(sqrt(eval2)*t)+b2*sin(sqrt(eval2)*t));
vm2=evect1(2)*(a1*cos(sqrt(eval1)*t)+b1*sin(sqrt(eval1)*t)) + evect2(2)*(a2*cos(sqrt(eval2)*t)+b2*sin(sqrt(eval2)*t));

%This was the analytical solution FOR PARTICULAR SYSTEM PARAMETERS
vm1_class=sin(t)+sin(sqrt(6)*t);
vm2_class=2*sin(t)-0.5*sin(sqrt(6)*t);

% PLOT TO COMPARE THE SOLUTIONS ('Class solution' for given BCs
plot(t,vm1,t,vm2,t,vm1_class,'x',t,vm2_class,'o')
legend('vm1 - Matlab', 'vm2 - Matlab', 'vm1 - Class', 'vm2 - class');
