%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Driver script for 'pendulum' and 'lpendulum' which solve the nonlinear and 
% linear pendulum equations using second-order finite difference techniques
% as discussed in class.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Dummy assignment statement so that octave knows that this is a script
% file, not a function file.
dummy = 0;

% Turn off the pager to ensure that all terminal output appears immediately
% and in the proper order.
more off;

infotp('Greetings student! Welcome to the pendulum demonstration script!');
infotp('');
infotp('Note that this script executes with ''more off'', so that from');
infotp('time to time you may have to scroll back through the terminal output');
infotp('to see material after it flashes by.');
infotp('');

% Examine definition of 'pendulum' which solves the nonlinear 
% pendulum equation using O(h^2) finite difference techniques 
% as discussed in class.  Recall that 'type <function name>'
% will display the definition of a function, as well as the 
% location of the .m file.
infotp('');
prompt('display definition of ''pendulum''');
type pendulum

% There is a similar function, lpendulum, which solves the 
% *linear* pendulum equation (i.e. the one that you study in 
% elementary physics courses.
infotp('');
prompt('display definition of ''lpendulum''');
type lpendulum

infotp('');
infotp('We will now perform some computations using ''pendulum'' and');
infotp('''lpendulum''.');
infotp('');
infotp('Note that in all cases below, we will set theta0 = 0');
infotp('(i.e. so that the pendulum is hanging vertically), and will treat the');
infotp('initial angular velocity, omega0 as our control parameter');
infotp('');

prompt('set values for tmax, theta0, level and omega0');

tmax = 40.0, theta0 = 0, level = 8, omega0 = 0.01,

infotp('');
infotp('Note that omega0 has been set to a very small value.');
infotp('We should therefore expect the dynamics of the nonlinear and');
infotp('linear pendula to be very similar.');
infotp('');

prompt('invoke ''pendulum'' and ''lpendulum''');

[t theta omega]   =  pendulum(tmax, level, theta0, omega0);
[t ltheta lomega] = lpendulum(tmax, level, theta0, omega0);

prompt('plot theta(t): nonlinear (red), linear (green)');
clf; hold on; plot(t, theta, 'r'); plot(t, ltheta, 'og');
prompt('continue');

infotp('');
infotp('We will now perform a basic convergence test of the solutions');
infotp('produced by ''pendulum'' by computing at discretization levels');
infotp('6, 7 and 8, using the same input parameters defined previously');
infotp('');
prompt('perform the 3-level convergence test');

infotp('Computing level 6 solution')
[t6 theta6 omega6]   =  pendulum(tmax, 6, theta0, omega0);

infotp('Computing level 7 solution')
[t7 theta7 omega7]   =  pendulum(tmax, 7, theta0, omega0);

infotp('Computing level 8 solution')
[t8 theta8 omega8]   =  pendulum(tmax, 8, theta0, omega0);

infotp('');
infotp('We will now plot theta(t) from all 3 calculations on a');
infotp('single plot level 6 (red)  level 7 (green)  level 8 (blue)');
infotp('');

prompt('make the plot');
clf; hold on; plot(t6, theta6, 'r-.o'); 
plot(t7, theta7, 'g-.+'); plot(t8, theta8, 'b-.*');
prompt('continue');

infotp('');
infotp('We will now filter the values of theta(t) from the level 7 and 8');
infotp('arrays so that they have the same length as the level 6 arrays');
infotp('(with values defined at the times in t6).');
infotp('We will then compute');
infotp('');
infotp('(dtheta67 = theta6 - theta7)  and (dtheta78 = theta7 - theta8)');
infotp('');
infotp('and make a plot of the resulting values:');
infotp('dtheta67 (red) dtheta78 (green)');
infotp('');

theta7 = theta7(1:2:length(theta7));
theta8 = theta8(1:4:length(theta8));

dtheta67 = theta6 - theta7;
dtheta78 = theta7 - theta8;

prompt('make the plot');
clf; hold on; plot(t6, dtheta67, 'r-.o'); plot(t6, dtheta78, 'g-.+');
prompt('continue');

infotp('');
infotp('We now rescale dtheta78 by 4 and replot with dtheta67.');
infotp('If the convergence is second order, the curves should be');
infotp('nearly coincident.');
infotp('');
dtheta78 = 4 * dtheta78;

prompt('make the plot');
clf; hold on; plot(t6, dtheta67, 'r-.o'); plot(t6, dtheta78, 'g-.+');

infotp('');
infotp('Note that the curves ARE nearly conincident.');
infotp('Also note how the magnitude of the error increases with time.');
infotp('');
prompt('continue');

infotp('');
infotp('We will now run a series of calculations at level 8 using');
infotp('increasingly larger values of omega0.');
infotp('');
infotp('Note that this means that we will be giving the pendula ever');
infotp('increasing initial angular velocities, so the resulting amplitudes');
infotp('of the oscillations will also increase.');
infotp('');

infotp('A key result that you can see immediately is, whereas the period');
infotp('(frequency) of the linear pendulum (green) stays constant');
infotp('(i.e. is NOT dependent on the amplitude of the oscillation),');
infotp('the period of the nonlinear pendulum (red) DOES depend on the');
infotp('amplitude (or, equivalently, on omega0).');
infotp('');

infotp('Resetting fixed input arguments for ''pendulum''');
tmax = 40, level = 8, theta0 = 0

prompt('continue');
infotp('');
infotp('We will now loop over the following values for omega0');
infotp('');
infotp('omega0 = [0.1 0.3 1.0 1.5 1.9]');
infotp('');
infotp('After each calculation we will plot theta(t) for the nonlinear (red)');
infotp('and linear (green) pendula, as well as a phase space plots,');
infotp('omega(t) vs theta(t) with the same colour scheme.');
infotp('');
for omega0 = [0.1 0.3 1.0 1.5 1.9]
   infotp('');
   fprintf('+++ omega0 = %g\n',omega0);
   infotp('');
    
   prompt('perform computation');
   [t theta omega]   =  pendulum(tmax, level, theta0, omega0);
   [t ltheta lomega] = lpendulum(tmax, level, theta0, omega0);

   infotp('Making theta(t) plots');
   clf; hold on; ptitle = sprintf('omega_0 = %g',omega0);
   plot(t, theta, 'r'); plot(t, ltheta, '-.og');
   title(ptitle);
   prompt('continue');

   infotp('Making phase-space plots');
   clf; hold on; ptitle = sprintf('Phase space plot: omega_0 = %g',omega0);
   plot(theta, omega, 'r'); plot(ltheta, lomega, '-.og');
   title(ptitle);
   prompt('continue');
end

infotp('');
infotp('Now, let''s give the nonlinear pendulum an even larger');
infotp('initial kick, by setting omega0 to 2.05');
infotp('');

omega0 = 2.05

infotp('');
prompt('perform the computation');

[t theta omega]   =  pendulum(tmax, level, theta0, omega0);

infotp('Plotting theta(t)');
clf; plot(t,theta,'r');
infotp('');
infotp('This curve is clearly different from the others we have seen.');
infotp('Can you explain this result physically?');
infotp('');
prompt('continue');
infotp('');
infotp('Note that we have now seen two distinct types of behaviour in');
infotp('the nonlinear oscillator as illustrated by the following ...');
infotp('');

[t thetalo omegalo]   =  pendulum(tmax, level, theta0, 1.9);
[t thetahi omegahi]   =  pendulum(tmax, level, theta0, 2.05);

clf; hold on; 
plot(t, thetalo, 'g'); plot(t, thetahi, 'r');
title('"lo(w)" and "hi(gh)" behaviour');
infotp('');
infotp('... so somewhere in the interval [1.9, 2.05] there should');
infotp('be a "critical value" of omega0 where the behaviour type');
infotp('changes.');
infotp('');
infotp('We will try to find this value using the technique of binary');
infotp('search, or bisection, but will leave this script to do so!');
infotp('');
