%% pendul - Program to compute the motion of a simple pendulum
% using the Euler or Verlet method
clear all;  help pendul      % Clear the memory and print header

%% * Select the numerical method to use: Euler or Verlet
NumericalMethod = menu('Choose a numerical method:', ...
                       'Euler','Verlet');
					   
%% * Set initial position and velocity of pendulum
theta0 = input('Enter initial angle (in degrees): ');
theta = theta0*pi/180;   % Convert angle to radians
omega = 0;               % Set the initial velocity

%% * Set the physical constants and other variables
g_over_L = 1;            % The constant g/L
time = 0;                % Initial time
irev = 0;                % Used to count number of reversals
period = [];             % Used to record period estimates
tau = input('Enter time step: ');

%% * Take one backward step to start Verlet
accel = -g_over_L*sin(theta);    % Gravitational acceleration
theta_old = theta - omega*tau + 0.5*tau^2*accel;    

%% * Loop over desired number of steps with given time step
%    and numerical method
nstep = input('Enter number of time steps: ');
for istep=1:nstep  

  %* Record angle and time for plotting
  t_plot(istep) = time;            
  th_plot(istep) = theta*180/pi;   % Convert angle to degrees
  time = time + tau;
  
  %* Compute new position and velocity using 
  %    Euler or Verlet method
  accel = -g_over_L*sin(theta);    % Gravitational acceleration
  if( NumericalMethod == 1 )
    theta_old = theta;               % Save previous angle
    theta = theta + tau*omega;       % Euler method
    omega = omega + tau*accel; 
  else  
    theta_new = 2*theta - theta_old + tau^2*accel;
    theta_old = theta;			   % Verlet method
    theta = theta_new;  
  end
  
  %* Test if the pendulum has passed through theta = 0;
  %    if yes, use time to estimate period
  if( theta*theta_old < 0 )  % Test position for sign change
    fprintf('Turning point at time t= %f \n',time);
    if( irev == 0 )          % If this is the first change,
      time_old = time;       % just record the time
    else
      period(irev) = 2*(time - time_old);
      time_old = time;
    end
    irev = irev + 1;       % Increment the number of reversals
  end
end

%% * Estimate period of oscillation, including error bar
AvePeriod = mean(period);
ErrorBar = std(period)/sqrt(irev);
fprintf('Average period = %g +/- %g\n', AvePeriod,ErrorBar);

%% * Graph the oscillations as theta versus time
figure(1); clf;    % Clear figure window #1 and bring it forward
plot(t_plot,th_plot,'+');
xlabel('Time');  ylabel('\theta (degrees)');