%         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

% Get the pendulum configuration from the user.
disp('Type in the two lengths and two masses, enclosed in [ ]s ');
disp('The sum of the lengths should be 1, as should the sum of the masses.');
disp('For example, you could enter:  [ .6 .4 .8 .2 ]');  
reply = input('  ');

L1 = reply(1);
L2 = reply(2);
M1 = reply(3);
M2 = reply(4);

% 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;

% Set the initial conditions.  The user specifies alpha and beta, and
% we'll assume that the pendulum starts at rest, so the "dot" variables are 0.
disp('Type in the initial values for alpha and beta, enclosed in [ ]s ');
disp('Values should be entered in radians.  For example, [ pi/3 pi/4 ]');
reply = input('  ');

alpha0   = reply(1);
beta0    = reply(2);
alphadot0= 0.;
betadot0 = 0.;

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

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('pendulum', 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);

% 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 );

x2 = x0 - L2*sin( beta );
y2 = y0 + L2*cos( beta );

% Compute the initial position.
initx = [ 0, x0(1), x1(1), x2(1) ];
inity = [ 0, y0(1), y1(1), y2(1) ];

plot( initx, inity, 'o', initx, inity, x0, y0, x1, y1 );

% 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 mass M1 in time.');
disp('Also shown is the initial position of the double pendulum.');

% ask the user if we should show the movie

yesno = input('See the movie for this motion?  (0 = no, 1 = yes): ');
if yesno == 1

   disp('Wait, making movies takes some space and time ' )
   
   % Free up some space for the movie storage
   clear z;

   % Number of time steps.  
   N = length(t);

   % Allocate space for the movie frames.
   % We will take only every third time step.
   MM = floor((N+2)/3 + 0.1 );
   M = moviein(MM);
   fprintf('Movie length = %d frames. ', MM);

   % Make the movie.  Each kk produces a frame.
   % The frames are stored by M(:,kk) = getframe;
   kk = 1;
   for k = 1:3:N
      x = [ 0, x0(k), x1(k), x2(k) ];
      y = [ 0, y0(k), y1(k), y2(k) ];
      plot( x, y, 'o', x, y, x1(1:k), y1(1:k));
      axis( [ -2.0, 2.0, -2.5, 1.5 ] );
      axis( 'square' );
      text(0.0,1.3,sprintf('t = %5.3f',t(k)));
      M(:,kk) = getframe;
      kk = kk + 1;
   end

   % display it 2 times
   disp(' ' );
   disp('The movie will be shown twice, hit control C to kill it. ' );
   disp('Type movie(M) to see the movie again ');

   movie(M,2);

end
