%  Script to plot leapfrog method with dt = pi/10, with and without
%  Asselin filtering with gamma = 0.05, applied to
%  the nonlinear pendulum equation: 
%    dphi(1)/dt = phi(2)
%    dphi(2)/dt = -sin(phi(1))
%  with IC of pendulum starting at rest 1 radian from the vertical:
%    phi(1) = 1, phi(2) = 0
%  Plots show phi(1) over roughly the 1-2 and the 49-50 pendulum
%  periods for each dt.

T = 100*pi;
dt = pi/10;
nt = round(T/dt);
phi0 = [1 0];

gamma = 0;  % No Asselin smoothing
[t,phi] = pendulum_leapfrog(T,dt,phi0,gamma);
subplot(2,2,1)
irange1 = 1 + round(trange1(1)/dt):round(trange1(2)/dt);
plot(t(irange1)/(2*pi), phi(irange1,1))
title('Leapfrog, \Delta t = \pi/10')
xlabel('t/2\pi')
ylabel('\phi_1')
subplot(2,2,2)
irange2 = 1 + round(trange2(1)/dt):round(trange2(2)/dt);
plot(t(irange2)/(2*pi), phi(irange2,1))
xlabel('t/2\pi')
ylabel('\phi_1')

gamma = 0.05; %  Add Asselin smoothing
[t,phi] = pendulum_leapfrog(T,dt,phi0,gamma);
subplot(2,2,3)
plot(t(irange1)/(2*pi), phi(irange1,1))
title('Leapfrog-Asselin, \Delta t = \pi/10, \gamma=0.05')
xlabel('t/2\pi')
ylabel('\phi_1')
subplot(2,2,4)
plot(t(irange2)/(2*pi), phi(irange2,1))
xlabel('t/2\pi')
ylabel('\phi_1')