function [t,phi] = pendulum_leapfrog(T,dt,phi0,gamma)

% function [t,phi] = pendulum_leapfrog(T,dt,phi0,gamma)
%
% Leapfrog method (with Asselin smoothing) applied to nonlinear pendulum eqn:
%   dphi(1)/dt = phi(2)                 
%   dphi(2)/dt = -sin(phi(1))
% Arguments:
%   T     final integration time 
%   dt    timestep  
%   phi0  initial phi at time t=0
%   gamma Asselin smoothing coefficient (0 for pure leapfrog)
% Outputs:
%   t     Vector of times 0:dt:T of length nt
%   phi   Matrix phi(nt,2) of solutions at the times t. 
  
F = inline('[phi(2) -sin(phi(1))]','phi');
nt = round(T/dt);
phi = zeros(nt,2);
phi(1,:) = phi0;

% Use forward Euler starting step (midpoint RK2 would be better)

phi(2,:) = phi(1,:) + dt*F(phi(2,:)); 

for n = 2:nt
  phi(n+1,:) = phi(n-1,:) + 2*dt*F(phi(n,:)); % Leapfrog step
  phi(n,:) = phi(n,:) ...
             + gamma*(phi(n+1,:)-2*phi(n,:)+phi(n-1,:)); % Asselin smooth
end

t = (0:nt)*dt;