function doall=pendulum_self_sustained
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos.

% A modified pendulum equation that allows for self-sustained
% oscillations even in the absence of forcing.  This allows us to see
% the full quasi-periodicity route to chaos, including the stage of
% the torus in phase space, unlike the normal pendulum which cannot
% display the torus, see homework 8.  Change the forcing amplitude
% below to see all stages in the route to chaos.  Solve numerically
% using the "modified Euler scheme" (a second order, constant step
% scheme), see Strogatz page 33.

% The equation:
% d^2\theta/dt^2=(g/L)*\sin(\theta) + gamma*([d\theta/dt]-[d\theta/dt]^3) + f*cos(omega_f*t)

%  First, translate to dynamical system of first order ODEs:
%  dx/dt = y
%  dy/dt = -(g/L)*sin(x) - gamma*y + f*cos(omega_f*z)
%  dz/dt = 1

% define parameters:
% ------------------
% pendulum parameters: gravity and length:
g=9.8;
L=9.8;
% Linear frequency:
omega_L=sqrt(g/L);
% Strength of friction (gamma1=1.0; gamma3=0.1; for chaos):
gamma1=1.0;
gamma3=0.1;
% f is forcing amplitude (3.0 gives chaos; 2.0 mode locked; 1.0 torus; 0.0 limit cycle):
f=3.0;
% Forcing frequency (1.5*sqrt(g/L) for chaos):
omega_f=1.5*sqrt(g/L);
% time step: needs to be small enough for a given problem:
delta_t=0.01*2*pi/omega_f;  % 0.01*2*pi/(max(omega_f,omega_L));
% Print (and save) every nprint steps:
nprint=1;
% Integration time:
T=1500.0;

% Number of steps:
N=floor(T/delta_t);
NP=floor(N/nprint);

% Initial conditions:
x=zeros(3,NP);
n=1;
x(:,n)=[pi/8.0,0,0]';

% Integrate:
for n=1:1:NP
    xp=x(:,n);
    for j=1:1:nprint
         x_tilde=xp+delta_t*rhs(xp,g,L,gamma1,gamma3,omega_f,f);
         xp=xp+0.5*(rhs(xp,g,L,gamma1,gamma3,omega_f,f)  ...
      +                rhs(x_tilde,g,L,gamma1,gamma3,omega_f,f))*delta_t;
    end
    x(:,n+1)=xp;
end

% Plot simple time series:
plot(x(3,:),x(2,:))
xlabel('time')
ylabel('d theta/dt (rad/time)')

% Plot reconstructed phase space using delay coordinates:
delay=floor(1.0/(delta_t*sqrt(g/L)));
b=floor(NP/2);
e=NP-3*delay;
figure
plot3(x(2,b:e),x(2,b+delay:e+delay),x(2,b+2*delay:e+2*delay))



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ff=rhs(x,g,L,gamma1,gamma3,omega_f,f)
% return the rhs of the dynamical system:

ff=[x(2),-(g/L)*sin(x(1)) + gamma1*x(2)-gamma3*x(2)^3 + f*cos(omega_f*x(3)),1.0]';