%% ME - 221 Dynamical systems, Spring 2025
% Solution to MATLAB PROBLEM SET 1
% Student Name: 
% SCIPER:
% Submission date:

%% Initiate the script
clc; clear; close all

%% Define the parameters
l_0 = ;    % m - COMPLETE
M = ;    % kg - COMPLETE
k = ;    % N/m - COMPLETE
c = ;    % Ns/m - COMPLETE
g = ;     % m/s^2 - COMPLETE

%% Q2 Simulation of dynamic response
% Harmonic force parameters
A = ;            % Force amplitude [N] - COMPLETE
w1 = ; %COMPLETE
w2 = ; %COMPLETE
w3 = ; %COMPLETE
w = [w1 w2 w3];	% rad/s
% Initial conditions
y0 = ; %COMPLETE
v0 = ; %COMPLETE
IC = [y0; v0];
tt = ;  % define the time vectors
opts = odeset('RelTol',1e-4,'AbsTol',1e-6);

% Simulate the three cases
[~,z_1] = ode45(@(t,z)  A1_P1_ode_2025(t,z,M,k,c,A,w(1),l_0,g), tt, IC, opts);
[~,z_2] = ode45(@(t,z)  A1_P1_ode_2025(t,z,M,k,c,A,w(2),l_0,g), tt, IC, opts);
[~,z_3] = ode45(@(t,z)  A1_P1_ode_2025(t,z,M,k,c,A,w(3),l_0,g), tt, IC, opts);

% Data visualization 
figure(1);
yyaxis left
    plot(tt,z_1(:,1),'linewidth',2);
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    set(gca,'FontSize',11)
yyaxis right
    plot(tt,z_1(:,2),'linewidth',2);
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',11)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 15 10])
exportgraphics(gcf,'Figure_Q2_1.png','Resolution',300);    % Save a png image of the first graph with a resolution of 300 dpi

figure(2);
%COMPLETE
figure(3);
%COMPLETE
%% Q4 Plot the rail model
L0 = ; %COMPLETE
l_01 = ;%COMPLETE
l_02 = ;%COMPLETE
l_03 = ;%COMPLETE
l_0V = [l_01 l_02 l_03];  
alpha = ; %COMPLETE
beta = ; %COMPLETE

xx = ; %COMPLETE
yy = ; %COMPLETE

for i = 1:3
    l0 = l_0V(i);
    V = (1/2).*k.*((-1).*l0+L0+alpha.*((-1)+cos(beta.*xx))).^2+alpha.*g.*M.*(( ...
          -1)+cos(beta.*xx));
    dVf = @(x) alpha.*beta.*((-1).*g.*M+k.*(alpha+l0+(-1).*L0+(-1).*alpha.*cos( ...
          beta.*x))).*sin(beta.*x);
    ddVf = @(x) alpha.*beta.^2.*((-1).*g.*M.*cos(beta.*x)+k.*cos(beta.*x).*( ...
          alpha+l0+(-1).*L0+(-1).*alpha.*cos(beta.*x))+alpha.*k.* ...
          sin(beta.*x).^2);

    %Finding the stable equilibrium points
    S = csapi(xx,???(xx));      %REPLACE ??? by appropriate anonymous function to get the non-discrete function which we want to find the zeros of.
    z = fnzeros(S);             % Find zeros of the function. To understand structure read DOC
    temp = ~(z(1,:)-z(2,:));    % Find duplicates.
    roots = z(1,temp).';        % Remove duplicates.
    roots(???(roots)<=0) = [];  %REPLACE ??? appropriate anonymous function to keep zeros which are minima
    
  
    figure;clf
    yyaxis left
    plot(,'linewidth',1.5);hold all %COMPLETE
    for jj = 1:length(roots)
        xline(roots(jj),'--k')
    end
        xlabel('$x$(m)','interpreter','Latex')
        ylabel('$y$(m)','interpreter','Latex')
        set(gca,'fontsize',8)
        xlim([xx(1) xx(end)])
    yyaxis right
    plot(,'linewidth',1.5);hold all %COMPLETE
        xlabel('$x$(m)','interpreter','Latex')
        ylabel('$V$(J)','interpreter','Latex')
        set(gca,'fontsize',8)
        set(gcf,'color','w','units','centimeters','outerposition',[5 5 12 8])
        xlim([xx(1) xx(end)])
    
    exportgraphics(gcf,strcat('Figure_Q4_',num2str(i),'.png'),'Resolution',300);  

end

%% Q5 Simulation
% Harmonic force parameters
l0 = ; %COMPLETE

V = (1/2).*k.*((-1).*l0+L0+alpha.*((-1)+cos(beta.*xx))).^2+alpha.*g.*M.*(( ...
      -1)+cos(beta.*xx));
dVf = @(x) alpha.*beta.*((-1).*g.*M+k.*(alpha+l0+(-1).*L0+(-1).*alpha.*cos( ...
      beta.*x))).*sin(beta.*x);
ddVf = @(x) alpha.*beta.^2.*((-1).*g.*M.*cos(beta.*x)+k.*cos(beta.*x).*( ...
      alpha+l0+(-1).*L0+(-1).*alpha.*cos(beta.*x))+alpha.*k.* ...
      sin(beta.*x).^2);

%Finding the stable equilibrium points
S = csapi(xx,???(xx));      %REPLACE ??? by appropriate anonymous function to get the non-discrete function which we want to find the zeros of.
z = fnzeros(S);             % Find zeros of the function. To understand structure read DOC
temp = ~(z(1,:)-z(2,:));    % Find duplicates.
roots = z(1,temp).';        % Remove duplicates.
roots(???(roots)<=0) = [];  %REPLACE ??? appropriate anonymous function to keep zeros which are minima

    
ICM = []; %COMPLETE

%% Initial conditions (1)
i=1;
IC = ICM(i,:)';
tt = ;  %COMPLETE
opts = ; %COMPLETE

% Simulate
[~,z_1] = ode45(@(t,z)  ); %COMPLETE

% Plot position and velocity of particle over time.
figure(7); clf
subplot(1,2,1)
yyaxis left
    plot(,'linewidth',1.5); %COMPLETE
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    for jj = 1:length(roots)
        yline(roots(jj),'--k')
    end
    set(gca,'FontSize',8)
    ylim([0.9*min(z_1(:,1)) 1.1*max(z_1(:,1))])
yyaxis right
    plot(,'linewidth',1.5); %COMPLETE
    ylim([1.1*min(z_1(:,2)) 1.1*max(z_1(:,2))])
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',8)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])

%Plot particle position on rail and potential energy in space.
figure(7);
subplot(1,2,2)
yyaxis left
plot(xx,yy,'linewidth',1.5);hold all
COL = cool(size(z_1,1));
for jj = 1:length(roots)
    xline(roots(jj),'--k')
end
for kk = 1:size(z_1,1)/50:size(z_1,1)
    plot(z_1(kk,1),alpha*(cos(beta*z_1(kk,1))-1),'o','color',COL(kk,:),'MarkerFaceColor',COL(kk,:));
end
    xlim([xx(1) xx(end)])
    xlabel('${x}$(m)','interpreter','Latex')
    ylabel('$y$(m)','interpreter','Latex')
    set(gca,'fontsize',8)
    set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])
yyaxis right
    plot(,'linewidth',2);hold all %COMPLETE
    xlabel('$x$(m)','interpreter','Latex')
    ylabel('$V$(J)','interpreter','Latex')
    set(gca,'fontsize',8)
    set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])
    xlim([xx(1) xx(end)]) 

exportgraphics(gcf,strcat('Figure_Q5_1.png'),'Resolution',300); 
        
%% Initial conditions (2)
i=2;
IC = ICM(i,:)';
opts = ; %COMPLETE

% Simulate
[~,z_1] = ode45(@(t,z)  ); %COMPLETE

% Plot position and velocity of particle over time.
figure(8); clf
subplot(1,2,1)
yyaxis left
    plot(,'linewidth',1.5); %COMPLETE
    xlabel('$t$(s)','interpreter','Latex')
    ylabel('$x$(m)','interpreter','Latex')
    for jj = 1:length(roots)
        yline(roots(jj),'--k')
    end
    set(gca,'FontSize',8)
    ylim([0.9*min(z_1(:,1)) 1.1*max(z_1(:,1))])
    set(gca,'FontSize',8)
yyaxis right
    plot(,'linewidth',1.5); %COMPLETE
    ylim([1.1*min(z_1(:,2)) 1.1*max(z_1(:,2))])
    ylabel('$\dot{x}$(m/s)','interpreter','Latex')
    set(gca,'FontSize',8)
set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])

%Plot particle position on rail and potential energy in space.
figure(8);
subplot(1,2,2)

yyaxis left
plot(,'linewidth',1.5);hold all %COMPLETE
COL = cool(size(z_1,1));
for jj = 1:length(roots)
    xline(roots(jj),'--k')
end
for kk = 1:size(z_1,1)/50:size(z_1,1)
    plot(z_1(kk,1),alpha*(cos(beta*z_1(kk,1))-1),'o','color',COL(kk,:),'MarkerFaceColor',COL(kk,:));
end
    xlim([xx(1) xx(end)])
    xlabel('${x}$(m)','interpreter','Latex')
    ylabel('$y$(m)','interpreter','Latex')
    set(gca,'fontsize',8)
    set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])

yyaxis right
    plot(,'linewidth',2);hold all %COMPLETE
    xlabel('$x$(m)','interpreter','Latex')
    ylabel('$V$(J)','interpreter','Latex')
    set(gca,'fontsize',8)
    set(gcf,'color','w','units','centimeters','outerposition',[5 5 12*2 8])
    xlim([xx(1) xx(end)]) 
exportgraphics(gcf,strcat('Figure_Q5_2.png'),'Resolution',300);