clear all
close all

%
% first, we declare the constants for the problem
%
% better display values for plotting
fontsize=20;
linewidth=2;
markersize=10;
speed=30;
logswitch=1;    %if 1, space betas logarithmically, if 0, space linearly

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

% global variables
global beta omega_o omega f_o
N = 300;
beta_factor = linspace(2/N,2,N)';
plt = @plot;
if logswitch
    beta_factor = logspace(-1,1,100)';
    plt = @semilogx;
end
%check out what it looks like when omega loops in logspace.  This way, you
%can clearly see the solution settle as omega -> infinity

omega_o = 1;
omega = 1;  %this is on resonance, but you could have chosen something else
beta = beta_factor*omega_o;
om_1 = sqrt(beta.^2 - omega_o^2);
decay(1:dsearchn(beta,omega_o)) = beta(1:dsearchn(beta,omega_o));
decay(dsearchn(beta,omega_o)+1:length(beta)) = beta(dsearchn(beta,omega_o)+1:length(beta)) - om_1(dsearchn(beta,omega_o)+1:length(beta));
delta = atan2(omega^2-omega_o^2,2*beta*omega);

f_o = 1;
Amax = f_o./(2*beta*omega_o);

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

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

time = linspace(t_start,t_finish,1000);
x_h = 0.5*x_o*real(exp((-beta+om_1)*time)+exp((-beta-om_1)*time));
x_f = f_o*cos(omega*time);
subplot(2,2,3)
plot(tspan,zeros(length(tspan),1),'k--','LineWidth',linewidth/2)
hold all
plot(time,x_f,'r:','LineWidth',linewidth);
h_4 = 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)
plt(beta,decay,'b-','LineWidth',linewidth)
hold all
h_5 = plt(beta,decay,'b-','LineWidth',linewidth);
xlabel('\beta [rad/s]','FontSize',fontsize)
ylabel('decay par. [rad/s]','FontSize',fontsize)
xlim([beta_factor(1) beta_factor(end)])
ylim([0 1.2*x_o])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

subplot(2,2,4)
plt(beta,Amax,'b-','LineWidth',linewidth)
hold all
h_6 = plt(beta,Amax,'b-','LineWidth',linewidth);
xlabel('\beta [rad/s]','FontSize',fontsize)
ylabel('A [m]','FontSize',fontsize)
xlim([beta_factor(1) beta_factor(end)])
ylim([0 1.2*max(Amax)])
set(gca,'FontSize',fontsize,'LineWidth',linewidth)

for i=1:length(beta_factor)
    beta = beta_factor(i)*omega_o;
    om_1 = sqrt(beta^2 - omega_o^2);
    %we're redefining delta here; first, it was delta(omega) for all omega, 
    %just for the sake of subplot(2,2,4), but now it will be delta(omega)
    %for the specific omega we're looping through
    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);
    
    subplot(2,2,1)
    delete(h_1)
    delete(h_2)
    delete(h_3)
    h_1 = plot(tspan,Amax(i)*ones(length(tspan),1),'r:','LineWidth',linewidth);
    h_2 = plot(tspan,-Amax(i)*ones(length(tspan),1),'r:','LineWidth',linewidth);
    h_3 = plot(t,x(:,1),'b-','LineWidth',linewidth);
    if Amax(i) > x_o
        ylim([-1.2*Amax(i) 1.2*Amax(i)])
    else
        ylim([-1.2*Amax(i) x_o])
    end

    x_h = 0.5*x_o*real(exp((-beta+om_1)*t)+exp((-beta-om_1)*t));   
    %instead of shifting the homogenous solution by -delta, we'll just
    %shift the forcing function by +delta as we loop through omega
    subplot(2,2,3)
    delete(h_4)
    h_4 = plot(t,x_h,'b-','LineWidth',linewidth);
    
    subplot(2,2,2)
    delete(h_5)
    h_5 = plt([beta,beta],[0,1.2*omega_o],'r:','LineWidth',linewidth);

    subplot(2,2,4)
    delete(h_6)
    h_6 = plt([beta,beta],[0,1.2*Amax(1)],'r:','LineWidth',linewidth);
    
    pause(speed*0.005)    %this is just to make the "movie" run slower
    
end