clear all
close all
%
% 4 POINTS FOR DECLARING ALL NECESSARY VARIABLES
% first, we declare the constants for the problem
%
% better display values for plotting
fontsize=20;
linewidth=2;
markersize=10;
speed=10;

% default values for plotting
% fontsize=12;
% linewidth=0.5;
% markersize=6;

% global variables
global beta omega_o omega f_o
omega_factor = linspace(0,2,300)';
%omega_factor = logspace(-1,1,100)';
beta = 0.1;
omega_o = 1;
om = omega_factor*omega_o;
om_1 = sqrt(omega_o^2 - beta^2);

f_o = 1;
Amax = f_o/(2*beta*omega_o);
A = f_o./sqrt((om.^2-omega_o^2).^2+(2*beta*om).^2);
delta = atan2(om.^2-omega_o^2,2*beta*om);

% initial parameters
x_o = 1;
v_o = 0;
t_start = 0;
t_finish = 100;
tspan=[t_start,t_finish];
init=[x_o,v_o];     %[x_o, (dx/dt)_o]

subplot(2,2,1)
plot(tspan,Amax*ones(length(tspan),1),'r:','LineWidth',linewidth)
hold all
plot(tspan,-Amax*ones(length(tspan),1),'r:','LineWidth',linewidth)
plot(tspan,zeros(length(tspan),1),'k--','LineWidth',linewidth/2)
h_1 = plot(tspan,-Amax*ones(length(tspan),1),'r:','LineWidth',linewidth);
xlabel('t [s]','FontSize',fontsize)
ylabel('x [m]','FontSize',fontsize)
ylim([-1.2*Amax 1.2*Amax])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

time = linspace(t_start,t_finish,100);
x_h = x_o*exp(-beta*time).*cos(om_1*time);
subplot(2,2,3)
plot(time,x_h,'b-','LineWidth',linewidth)
hold all
plot(tspan,zeros(length(tspan),1),'k--','LineWidth',linewidth/2)
h_2 = plot(time,x_h,'b-','LineWidth',linewidth);
xlabel('t [s]','FontSize',fontsize)
ylabel('x_h and x_f [m]','FontSize',fontsize)
ylim([-1.2*x_o 1.2*x_o])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(2,2,2)
plot(om,A,'b-','LineWidth',linewidth)
hold all
h_3 = plot(om,A,'b-','LineWidth',linewidth);
xlabel('\omega [rad/s]','FontSize',fontsize)
ylabel('A [m]','FontSize',fontsize)
xlim([0 max(om)])
ylim([0 1.2*Amax])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(2,2,4)
plot(om,delta,'b-','LineWidth',linewidth)
hold all
plot(tspan,zeros(length(tspan),1),'k--','LineWidth',linewidth/2)
h_4 = plot(om,delta,'b-','LineWidth',linewidth);
xlabel('\omega [rad/s]','FontSize',fontsize)
ylabel('\delta [rad]','FontSize',fontsize)
xlim([0 max(om)])
ylim([-pi/2 pi/2])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)
    
for i=1:length(omega_factor)
    omega = omega_factor(i)*sqrt(omega_o^2-beta^2);
    % 10 POINTS FOR CORRECT FUNCTION CALL AND CORRECT FUNCTION INTERNALS
    % options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6]);   %Matlab default values
    options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6]);
    [t,x]=ode45(@harmonic_oscillator,tspan,init,options);
    %[t,x]=ode45(@harmonic_oscillator_quad,tspan,init,options);

    % 4 POINTS FOR ACTUALLY PLOTTING WHAT I ASKED FOR
    
    subplot(2,2,1)
    delete(h_1)
    h_1 = plot(t,x(:,1),'b-','LineWidth',linewidth);

    x_f = f_o*cos(omega*t);
    subplot(2,2,3)
    delete(h_2)
    h_2 = plot(t,x_f,'r:','LineWidth',linewidth);

    subplot(2,2,2)
    delete(h_3)
    h_3 = plot([omega,omega],[0,1.2*Amax],'r:','LineWidth',linewidth);

    subplot(2,2,4)
    delete(h_4)
    h_4 = plot([omega,omega],[-pi,pi],'r:','LineWidth',linewidth);
    
    pause(speed*0.005)
    
end