%         This file contains the matlab code to solve the
%         system of ODEs describing the double pendulum.

%         Also need file pendulum.m

%              Prof. J. Strikwerda  Oct. 31, 1995
%              modified by Kirk Hogenson, Nov. 11, 1997

% Clear variables.
clear;

% Global variables are needed to get information to the function pendulum.
global B
global GG
global Lbar
global L0 
global L1

L1 = 0.5;
L2 = 0.5;
M1 = 0.7;
M2 = 0.3;

% Set the global variables.
L0   = 1;
B    = ( L1*M1 - L2*M2 ) / ( M1 + M2 );
Lbar = ( M1*L1*L1 + M2*L2*L2 ) / ( M1 + M2 ) / L0;

% Acceleration of gravity in meters/sec^2
GG   = 9.8;

% Set the initial and final time.
t0       =  0.0;
tfinal   = 10.0;

alpha0   = pi/2;
beta0    = pi/4;
gamma0   = 0;
alphadot0= 0.;
betadot0 = 0.;
gammadot0 = 0.;

% Put in standard form  z' = f( t, z ) 
z0 = [ alpha0, beta0, gamma0, alphadot0, betadot0, gammadot0 ];

disp('Wait, this may take a while...' );

% Now solve it!  Use Matlab's automatic ODE solver, ode45.
% The function  pendulum  is in the file pendulum.m

tspan = [t0, tfinal];
[t,z] = ode45('pendulum5', tspan, z0);

% Plot the variables alpha, beta, alpha', beta' as functions of t .
plot( t, z );

disp(' ');
disp('This first plot shows the variation of alpha, beta, alphadot, and betadot');
disp('as functions of t. ');
disp(' ');
disp('Hit any key to continue...');
pause;

% Plot the initial conditions and the solution in the (x,y) plane.

% Get the two angles.
alpha = z(:,1);
beta  = z(:,2);
gamma = z(:,3);

% Compute the positions of the pendulum.  (x0,y0) is the position
% of the point P1, (x1,y1) is the position of M1, and (x2,y2) is the
% position of M2.

x0 =  L0* sin( alpha );
y0 = -L0* cos( alpha );

x1 = x0 + L1*sin( beta );
y1 = y0 - L1*cos( beta );

x3 = x1 + L1*sin( gamma );
y3 = y1 - L1*cos( gamma );

% Compute the initial position.
initx = [ x1(1), x3(1) ];
inity = [ y1(1), y3(1) ];

plot( initx, inity, 'o', initx, inity, x1, y1, x3, y3, 'o' );

% The axes are set to have all the plots be on the same axes.
axis( [ -2.0, 2.0, -2.5, 1.5 ] );

% Have the axes have the same length scales.
axis( 'square' );

disp(' ');
disp('This second plot shows the path of the pendulum and masses M3 ');
disp('(circles) and M1 in time.');
disp(' ');
disp('Also shown is the initial position of the double pendulum.');

