% Hayden Molloy
% AERO3220 HW4
% Due 18 Feb 2026

clear; close all; clc;


%% Initial Parameters

    k1 = 100;      % N/m
    k2 = 50;       % N/m
    m1 = 500;      % kg
    m2 = 250;      % kg

% Time Parameters
    t0 = 0;
    tf = 200;
    dt = 0.01;
    tspan = t0:dt:tf;

% Initial conditions [z1 z1dot z2 z2dot]
    x0 = [0; 0; 0; 0];

% Initialize Force Functions
    F_1  = @(t) 100;                 % N
    F_2  = @(t) 100*sin(0.25*t);     % N
    F_3  = @(t) 100*sin(0.50*t);     % N


%% Define Cases and ODEs

% CASE 1:
    b = 100;
    F = F_1;
[t1,x1] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);

% Case 2:
    b = 100;
    F = F_2;
[t2,x2] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);

% CASE 3:
    b = 0;
    F = F_2;
[t3,x3] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);

% Case 4:
    b = 100;
    F = F_3;
[t4,x4] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);

% Case 5:
    b = 0;
    F = F_3;
[t5,x5] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);


%% Report Results

% Create Legend
lgnd = { ...
    'Case 1: F=100, b=100', ...
    'Case 2: F=100sin(0.25t), b=100', ...
    'Case 3: F=100sin(0.25t), b=0', ...
    'Case 4: F=100sin(0.50t), b=100', ...
    'Case 5: F=100sin(0.50t), b=0'};


% Plot 1: Top mass displacement vs. Time
figure(1);

plot(t1,x1(:,1),'LineWidth',1.3); 
hold on;
plot(t2,x2(:,1),'LineWidth',1.3);
plot(t3,x3(:,1),'LineWidth',1.3);
plot(t4,x4(:,1),'LineWidth',1.3);
plot(t5,x5(:,1),'LineWidth',1.3);
grid on;
xlabel('Time (s)');
ylabel('z_1 (m)');
title('Top Mass Displacement z_1(t)');
legend(lgnd,'Location','southoutside');


% Plot 2: Bottom mass displacement vs. Time
figure(2);

plot(t1,x1(:,3),'LineWidth',1.3); 
hold on;
plot(t2,x2(:,3),'LineWidth',1.3);

plot(t3,x3(:,3),'LineWidth',1.3);

plot(t4,x4(:,3),'LineWidth',1.3);

plot(t5,x5(:,3),'LineWidth',1.3);

grid on;
xlabel('Time (s)');
ylabel('z_2 (m)');
title('Bottom Mass Displacement z_2(t)');
legend(lgnd,'Location','southoutside');


% Plot 3: Top mass velocity vs. Time
figure(3);

plot(t1,x1(:,2),'LineWidth',1.3); 
hold on;
plot(t2,x2(:,2),'LineWidth',1.3);

plot(t3,x3(:,2),'LineWidth',1.3);

plot(t4,x4(:,2),'LineWidth',1.3);

plot(t5,x5(:,2),'LineWidth',1.3);

grid on;
xlabel('Time (s)');
ylabel('dz_1/dt (m/s)');
title('Top Mass Velocity z_1''');
legend(lgnd,'Location','southoutside');


% PLOT 4: Bottom mass velocity vs. Time
figure(4);

plot(t1,x1(:,4),'LineWidth',1.3); 
hold on;
plot(t2,x2(:,4),'LineWidth',1.3);

plot(t3,x3(:,4),'LineWidth',1.3);

plot(t4,x4(:,4),'LineWidth',1.3);

plot(t5,x5(:,4),'LineWidth',1.3);

grid on;
xlabel('Time (s)');
ylabel('dz_2/dt (m/s)');
title('Bottom Mass Velocity z_2''');
legend(lgnd,'Location','southoutside');


% PLOT 5: z1(t) and z2(t)
% (solid = z1, dashed = z2)

c1 = [0 114 178]/255;
c2 = [213 94 0]/255;
c3 = [230 159 0]/255;
c4 = [117 112 179]/255;
c5 = [0 158 115]/255;

figure(5);

plot(t1,x1(:,1),'-','LineWidth',1.2, 'Color', c1); 
hold on;
plot(t1,x1(:,3),'--','LineWidth',1.2, 'Color', c1);

plot(t2,x2(:,1),'-','LineWidth',1.2, 'Color', c2);
plot(t2,x2(:,3),'--','LineWidth',1.2, 'Color', c2);

plot(t3,x3(:,1),'-','LineWidth',1.2, 'Color', c3);
plot(t3,x3(:,3),'--','LineWidth',1.2, 'Color', c3);

plot(t4,x4(:,1),'-','LineWidth',1.2, 'Color', c4);
plot(t4,x4(:,3),'--','LineWidth',1.2, 'Color', c4);

plot(t5,x5(:,1),'-','LineWidth',1.2, 'Color', c5);
plot(t5,x5(:,3),'--','LineWidth',1.2, 'Color', c5);

grid on;
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Top and Bottom Mass Displacements');

legend({ ...
    'Case 1 z_1', ...
    'Case 1 z_2', ...
    'Case 2 z_1', ...
    'Case 2 z_2', ...
    'Case 3 z_1', ...
    'Case 3 z_2', ...
    'Case 4 z_1', ...
    'Case 4 z_2', ...
    'Case 5 z_1', ...
    'Case 5 z_2'}, ...
    'Location','eastoutside');


%% Local ODE function
function dx = eom_two_mass(t, x, m1, m2, k1, k2, b, Ffun)
    z1    = x(1);
    z1dot = x(2);
    z2    = x(3);
    z2dot = x(4);

    u = Ffun(t);

    % m1*z1ddot + k1*z1 - k1*z2 = 0
    z1ddot = (-k1/m1)*z1 + (k1/m1)*z2;

    % m2*z2ddot + b*z2dot + (k1+k2)*z2 - k1*z1 = u
    z2ddot = (k1/m2)*z1 - ((k1+k2)/m2)*z2 - (b/m2)*z2dot + (1/m2)*u;

    dx = [z1dot; z1ddot; z2dot; z2ddot];
end