%% clear all; close all;
%% Eli Tziperman: example for course in nonlinear dynamics and chaos.
%% improvements by Lauren Kuntz

%% Solve simple periodically driven pendulum problem numerically
%% using ode45, plot time series, phase space, reconstructed phase
%% space. 

%% The equation:
%% d^2\theta/dt^2 + gamma*(d\theta/dt) + (g/L)*\sin(\theta) = a*cos(omega_f*t)

%%  First, transform to a system of 1st order ODEs, with x=theta:
%%  dx/dt = y
%%  dy/dt = -(g/L)*sin(x) - gamma*y + a*cos(omega_f*z)
%%  dz/dt = 1

%% define parameters:
% -------------------
%% pendulum parameters: 
%% calculate gravity & length from linear frequency using omega=sqrt(g/l):
g=9.8;
L=9.8;
%% linear period:
omega_L=sqrt(g/L);
period_L=2*pi/omega_L;
%% specify forcing period:
omega_f=pi/5;
period_f=2*pi/omega_f;
%% Strength of friction:
gamma=1/2;
%% forcing amplitude:
% a=0.2; MarkerSize_poincare=12; % periodic
 a=1.6; MarkerSize_poincare=1; % chaos
%% length of integration:
T_end=5000*pi;
%% time to start plotting, allow plotting only part after transient
%% decay when relevant:
T_start_plot=4950*pi;
T_start_plot_poincare=500*pi;
%% increment of output times from ode45:
dt=0.01;
%% Start plotting after transients dissipate
i_beg = round(T_start_plot/dt);
i_beg_poincare = round(T_start_plot_poincare/dt);
%% use forcing period as sub-sampling time for Poincare sections:
T_poincare_subsample=period_f;
i_poincare_subsample=round(T_poincare_subsample/dt);
%% specify delay for phase space reconstructions:
reconstruction_delay_time=0.1*period_f;
i_delay=round(reconstruction_delay_time/dt);

fprintf(1,'solving driven pendulum with following parameters:\n')
fprintf(1,['linear period=%g, forcing period=%g, friction time scale=%g,\n' ...
           'start plot time=%g, reconstruction delay time=%g\n'] ...
        ,period_L,period_f,1/gamma ...
        ,T_start_plot,reconstruction_delay_time);

%% the rhs of the ODEs to integrate:
F=@(t,x) [x(2); -(g/L)*sin(x(1))-gamma*x(2)+a*cos(omega_f*x(3)); 1];

%% Specify error tolerance:
opts=odeset('reltol',1.e-6);

%% Solve using ode45
initial_conditions=[pi-0.01; 0; 0];
output_times=[0:dt:T_end];
[t,x]=ode45(F,output_times,initial_conditions,opts);
x(:,1)=mod(x(:,1)+pi,2*pi)-pi;

%% Plot time series:
time=t;
set(0,'defaulttextfontsize',14); set(0,'defaultaxesfontsize',14);
figure(1); clf;

subplot(2,3,1);
plot(time(i_beg:end)-T_start_plot,x(i_beg:end,2));
xlim([time(i_beg)-T_start_plot,time(end)-T_start_plot]);
xlabel('time')
ylabel('d\theta/dt')
title('time series')

subplot(2,3,4);
plot(time(i_beg:end)-T_start_plot,x(i_beg:end,1),'.','MarkerSize',1);
xlim([time(i_beg)-T_start_plot,time(end)-T_start_plot]);
xlabel('time')
ylabel('\theta')
title('time series')

%% Plot phase space:
subplot(2,3,2);
plot(x(i_beg:end,1),x(i_beg:end,2),'.','MarkerSize',1);
xlim([min(x(i_beg:end,1)),max(x(i_beg:end,1))]);
xlabel('\theta');
ylabel('d\theta/dt')
title('phase space')

subplot(2,3,5);
plot(x(i_beg:end,2) ...
      ,x(i_beg-i_delay:end-i_delay,2));
%      ,x(i_beg-2*i_delay:end-2*i_delay,2));
xlabel('$\dot{\theta}(t)$','interpreter','LaTeX');
ylabel('$\dot{\theta}(t-\tau)$','interpreter','LaTeX')
%zlabel('$\dot{\theta}(t-2\tau)$','interpreter','LaTeX')
title('reconstructed phase space')

%% Plot Poincare sections in phase space:
subplot(2,3,3);
plot(x(i_beg_poincare:i_poincare_subsample:end,1),x(i_beg_poincare:i_poincare_subsample:end,2),'.','MarkerSize',MarkerSize_poincare)
xlim([min(x(i_beg_poincare:end,1)),max(x(i_beg_poincare:end,1))]);
ylim([min(x(i_beg_poincare:end,2)),max(x(i_beg_poincare:end,2))]);
xlabel('\theta');
ylabel('d\theta/dt');
title('poincare section');

subplot(2,3,6);
plot(x(i_beg_poincare:i_poincare_subsample:end,2) ...
     ,x(i_beg_poincare-i_delay:i_poincare_subsample:end-i_delay,2) ...
     ,'.','MarkerSize',MarkerSize_poincare);
%     ,x(i_beg_poincare-2*i_delay:i_poincare_subsample:end-2*i_delay,2) ...
xlim([min(x(i_beg_poincare:end,2)),max(x(i_beg_poincare:end,2))]);
ylim([min(x(i_beg_poincare:end,2)),max(x(i_beg_poincare:end,2))]);
%zlim([min(x(i_beg_poincare:end,2)),max(x(i_beg_poincare:end,2))]);
xlabel('$\dot{\theta}(t)$','interpreter','LaTeX'); 
ylabel('$\dot{\theta}(t-\tau)$','interpreter','LaTeX'); 
%zlabel('$\dot{\theta}(t-2\tau)$','interpreter','LaTeX'); 
title('reconstructed poincare section');

%% save as pdf:
set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [11 8]);
set(gcf, 'PaperPosition', [0 0 11 8]); % [left, bottom, width, height];
saveas(gcf,sprintf('Figure-driven-pendulum.pdf'));
