%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Driver script for 'pendulume' which solves the nonlinear pendulum 
% equation using second-order finite difference techniques as discussed 
% in class, but which also computes kinetic, potential and total energy.
%
% Performs the following calculations
%
% 1) Basic test of pendulume, including plots of all three energy 
%    quantites.
% 
% 2) Convergence test using levels 6, 7, 8 and 9 that plots deviations in
%    total energy relative to t = 0.
% 
% 2) Convergence test using levels 6, 7, 8 and 9 that plots scaled 
%    deviations in total energy relative to t = 0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% Turn paging off.
more off;

tmax = 40.0, theta0 = 0, omega0 = 1.9;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) Basic test (levels 6, 7, 8)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform level 6 calculation and plot energy quantities.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t theta omega T V E] =  pendulume(tmax, 6, theta0, omega0);

clf;
hold on;
plot(t, T, 'r');
plot(t, V, 'g');
plot(t, E, 'b');
disp('');
disp('Level 6 calculation');
disp('Plot of energy quantities');
disp('Red: Kinetic');
disp('Green: Potential');
disp('Blue: Total');
disp('');
input('Type ENTER to continue: ');
disp('');

% Compute deviation in total energy, relative to initial value and plot.
dE6 = E - E(1);
clf;
plot(t, dE6, 'b');
disp('Level 6 calculation');
disp('Plot of deviation in total energy');
disp('');
input('Type ENTER to continue: ');
disp('');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform level 7 calculation and plot energy quantities.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t theta omega T V E] =  pendulume(tmax, 7, theta0, omega0);

clf;
hold on;
plot(t, T, 'r');
plot(t, V, 'g');
plot(t, E, 'b');
disp('');
disp('Level 7 calculation');
disp('Plot of energy quantities');
disp('Red: Kinetic');
disp('Green: Potential');
disp('Blue: Total');
disp('');
input('Type ENTER to continue: ');
disp('');

% Compute deviation in total energy, relative to initial value and plot.
dE7 = E - E(1);
clf;
plot(t, dE7, 'b');
disp('Level 7 calculation');
disp('Plot of deviation in total energy');
disp('');
input('Type ENTER to continue: ');
disp('');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Perform level 8 calculation and plot energy quantities.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t theta omega T V E] =  pendulume(tmax, 8, theta0, omega0);

clf;
hold on;
plot(t, T, 'r');
plot(t, V, 'g');
plot(t, E, 'b');
disp('');
disp('Level 8 calculation');
disp('Plot of energy quantities');
disp('Red: Kinetic');
disp('Green: Potential');
disp('Blue: Total');
disp('');
input('Type ENTER to continue: ');
disp('');

% Compute deviation in total energy, relative to initial value and plot.
dE8 = E - E(1);
clf;
plot(t, dE8, 'b');
disp('Level 8 calculation');
disp('Plot of deviation in total energy');
disp('');
input('Type ENTER to continue: ');
disp('');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2) Convergence test: Recompute at levels 6, 7, 8, 9, plot deviations in 
%    total energy  and plot.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t6 theta omega T V E] =  pendulume(tmax, 6, theta0, omega0);
dE6 = E - E(1);
[t7 theta omega T V E] =  pendulume(tmax, 7, theta0, omega0);
dE7 = E - E(1);
[t8 theta omega T V E] =  pendulume(tmax, 8, theta0, omega0);
dE8 = E - E(1);
[t9 theta omega T V E] =  pendulume(tmax, 9, theta0, omega0);
dE9 = E - E(1);

clf;
hold on;
plot(t6, dE6, 'b');
plot(t7, dE7, 'g');
plot(t8, dE8, 'r');
plot(t9, dE9, 'm');
disp('');
disp('Plot of deviations in total energy');
disp('Blue: Level 6');
disp('Green: Level 7');
disp('Red: Level 8');
disp('Magenta: Level 9');
disp('');
input('Type ENTER to continue: ');
disp('');


% Scale energy deviations and replot.

dE7 = 4 * dE7;
dE8 = 16 * dE8;
dE9 = 64 * dE9;

clf;
hold on;
plot(t6, dE6, 'r');
plot(t7, dE7, 'g');
plot(t8, dE8, 'b');
plot(t9, dE9, 'm');
disp('Plot of scaled deviations in total energy');
disp('Blue: Level 6');
disp('Green: 4 * Level 7');
disp('Red: 16 * Level 8');
disp('Magenta: 64 * Level 9');
disp('Increasing coincidence of curves as step size decreases');
disp('signals second order convergence');
disp('');
input('Type ENTER to continue: ');
disp('');
