function double_pendulum_movie_small ( )

%*****************************************************************************80
%
%% double_pendulum_movie_small runs the small perturbation test.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2020
%
%  Author:  
%
%    John Burkardt.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'double_pendulum_movie_small\n' );
  fprintf ( 1, '  Animate the double pendulum small perturbation problem.\n' );
%
%  Set initial conditions.
%
  phi1 = 0.25;
  dtphi1 = 0.0;
  phi2 = 0.0;
  dtphi2 = 0.0;
%
%  Set problem parameters.
%
  g = 9.81; 
  m1 = 1.0; 
  m2 = 1.0; 
  l1 = 2.0; 
  l2 = 1.0;
%
%  Set movie parameters.
%
  duration = 100;
  fps = 10;
  movie = true;
  label = 'double_pendulum_movie_small';
%
%  Pack initial conditions and parameters into one array.
%
  ivp = [ phi1; dtphi1; phi2; dtphi2; g; m1; m2; l1; l2 ];
%
%  Call for solution and movie.
%
  double_pendulum_movie ( ivp, duration, fps, movie, label );

  return
end
function double_pendulum_movie ( ivp, duration, fps, movie, label )

%*****************************************************************************80
%
%% double_pendulum_movie animates a double pendulum solution.
%
%  Discussion:
%
%    This function is usually called by double_pendulum_init.
%    It calls double_pendulum_deriv().
%
%    Example direct call:
%      double_pendulum_movie([pi;0;pi;5;9.81;1;1;2;1],100,10,false)
%    Or, to call with pre-chosen data:
%      double_pendulum_init
%
%  Licensing:
%
%    Copyright (c) 2010, Alexander Erlich
%    All rights reserved.
%    
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions are
%    met:
%    
%        * Redistributions of source code must retain the above copyright
%          notice, this list of conditions and the following disclaimer.
%        * Redistributions in binary form must reproduce the above copyright
%          notice, this list of conditions and the following disclaimer in
%          the documentation and/or other materials provided with the distribution
%    
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%    ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%    LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%    
%  Modified:
%
%    23 March 2020
%
%  Author:  
%
%    Original MATLAB version by Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt.
%
%  Input:
%   
%    real ivp9): [ phi1; dtphi1; phi2; dtphi2; g; m1; m2; l1; l2]
%    Data defining the initial value problem. 
%    phi1 and dtphi1 are the initial angle and angular velocity for mass 1.
%    phi2 and dtphi2 are the initial angle and angular velocity for mass 2. 
%    g is the gravitational coefficient;
%    m1 and m2 are the pendulum masses.
%    l1 and l2 are the pendulum lengths.
%  
%    real duration: the time interval over which the ODE is to be solved,
%    in seconds.
%
%    integer fps: Frames Per Second. The framerate is relevant both for 
%    normal real time animation and movie recording.
%
%    logical movie: if false, a normal realtime animation of the motion 
%    of the double pendulum is shown, the framerate being fps.
%    If true, an AVI movie is recorded. 
%
%    string label: a label that is used to construct the movie filename.
%    The movie filename will be label+'.avi'.
%

%
%  Use ode45 to compute the solution.
%
  f = @ double_pendulum_deriv;
  tspan = [ 0.0, duration ];
  sol = ode45 ( f, tspan, ivp );
%
%  Interpolate solution values at NFRAMES equally spaced times.
%
  nframes = duration * fps;
  t = linspace ( 0.0, duration, nframes );
  y = deval ( sol, t );
%
%  Extract solution data.
%
  phi1 = y(1,:)';
  dtphi1 = y(2,:)';
  phi2 = y(3,:)';
  dtphi2 = y(4,:)';
  l1 = ivp(8);
  l2 = ivp(9);
%
%  Open the movie file.
%
  if ( movie )
    filename = sprintf ( '%s.avi', label );
    vw = VideoWriter ( filename );
    open ( vw );
  end
%
%  Prepare the graphics.
% 
  clf ( );
  fig = plot ( 0.0, 0.0, 'MarkerSize', 30, 'Marker', '.', 'LineWidth', 2 );

  range = 1.1 * ( l1 + l2 );
  axis ( [ -range, range, -range, range ] );
  axis ( 'square' );
  grid ( 'on' );
  viscircles ( [ 0.0, 0.0; 0.0, 0.0; 0.0, 0.0 ], [ l1, l1+l2, abs(l1-l2) ] );
  set ( gca, 'nextplot', 'replacechildren' );
%
%  Display each frame.
%
  for i = 1 : nframes

    x1 =   l1 * sin ( phi1(i) );
    y1 = - l1 * cos ( phi1(i) );
    x2 =   l1 * sin ( phi1(i) ) + l2 * sin ( phi2(i) );
    y2 = - l1 * cos ( phi1(i) ) - l2 * cos ( phi2(i) );

    set ( fig, 'XData', [ 0.0, x1, x2 ], 'YData', [ 0.0, y1, y2 ] );
    label2 = s_escape_tex ( label );
    plot_label = sprintf ( '%s frame %d', label2, i );
    title ( plot_label );

    frame = getframe ( gcf );
    writeVideo ( vw, frame );
    if ( ~ movie )
      pause ( t(i+1) - t(i) );
    end

  end

  if ( movie )
    close ( vw );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Animation saved as "%s"\n', filename );
  end

  return
end
function xdot = double_pendulum_deriv ( t, x )

%*****************************************************************************80
%
%% double_pendulum_deriv defines the double pendulum ODE.
%
%  Discussion:
%
%    This function uses a simple trick to pass in 5 problem constants
%    along with the 4 problem variables.
%
%    The input vector x contains:
%    * x(1): phi1, the angle of the first pendulum
%    * x(2): dphi1dt, the angular velocity of the first pendulum
%    * x(3): phi2, the angle of the second pendulum
%    * x(4): dphi2dt, the angular velocity of the second pendulum
%    * x(5): g, the constant gravitational coefficient
%    * x(6): m1, the constant mass 1
%    * x(7): m2, the constant mass 2
%    * x(8): l1, the constant length 1
%    * x(9): l2, the constant length 2
%
%    x(5) through x(9) are given dummy ODE's of the form dxdt = 0.
%
%  Licensing:
%
%    Copyright (c) 2010, Alexander Erlich
%    All rights reserved.
%    
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions are
%    met:
%    
%        * Redistributions of source code must retain the above copyright
%          notice, this list of conditions and the following disclaimer.
%        * Redistributions in binary form must reproduce the above copyright
%          notice, this list of conditions and the following disclaimer in
%          the documentation and/or other materials provided with the distribution
%    
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%    ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%    LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%    
%  Modified:
%
%    23 March 2020
%
%  Author:  
%
%    Original MATLAB version by Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt.
%
%  Input:
%
%    real T: the time. 
%
%    real X(9): the current solution.
%
%  Output:
%
%    real xdot(9): the right hand side of the ODEs.
%

%
%  Retrieve the constants.
%
  g = x(5);
  m1 = x(6);
  m2 = x(7);
  l1 = x(8);
  l2 = x(9);

  xdot = zeros ( 9, 1 );

  xdot(1) = x(2);

  xdot(2) = - ...
  ( ...
    ( ...
      g * ( 2.0 * m1 + m2 ) * sin ( x(1) ) + m2 * ...
      ( ...
        g * sin ( x(1) - 2.0 * x(3) ) + 2.0 * ...
        ( ...
          l2 * x(4)^2 + l1 * x(2)^2 * cos ( x(1) - x(3) ) ...
        ) ...
        * sin ( x(1) - x(3) ) ...
      ) ...
    ) ...
    / ...
    ( ...
      2.0 * l1 * ( m1 + m2 - m2 * cos ( x(1) - x(3) )^2 ) ...
    ) ...
  );

  xdot(3) = x(4);

  xdot(4) = ...
  ( ...
    ( ...
      ( m1 + m2 ) * ...
      ( ...
        l1 * x(2)^2 + g * cos ( x(1) ) ...
      ) ...
      + l2 * m2 * x(4)^2 * cos ( x(1) - x(3) ) ...
    ) ...
    * sin ( x(1) - x(3) ) ...
  ) ...
  / ...
  ( ...
    l2 * ...
    ( ...
      m1 + m2 - m2 * cos ( x(1) - x(3) )^2 ...
    ) ...
  );

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% S_ESCAPE_TEX de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string S1: the string to be de-escaped.
%
%  Output:
%
%    string S2: a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end

