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:
%
%    16 November 2020
%
%  Author:  
%
%    Original MATLAB version by Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'double_pendulum_movie_small\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Animate the double pendulum ODE solution.\n' );
%
%  Set the initial conditions.
%
  theta1 = 0.25;
  dtheta1dt = 0.0;
  theta2 = 0.0;
  dtheta2dt = 0.0;
%
%  Set problem parameters.
%
  [ g, m1, m2, l1, l2 ] = double_pendulum_parameters ( );
%
%  Set movie parameters.
%
  t0 = 0.0;
  tstop = 100;
  tspan = [ t0, tstop ];

  fps = 10;
  movie = true;
  label = 'double_pendulum_movie_small';
%
%  Pack initial conditions and parameters into one array.
%
  y0 = [ theta1; dtheta1dt; theta2; dtheta2dt ];
%
%  Call for solution and movie.
%
  double_pendulum_movie ( tspan, y0, fps, movie, label );

  return
end
function double_pendulum_movie ( tspan, y0, fps, movie, label )

%*****************************************************************************80
%
%% double_pendulum_movie animates a double pendulum solution.
%
%  Discussion:
%
%    Example direct call:
%      double_pendulum_movie([pi;0;pi;5],100,10,false)
%
%  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:
%
%    16 November 2020
%
%  Author:  
%
%    Original MATLAB version by Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt.
%
%  Input:
%
%    real tspan(2): the endpoints of the time interval, in seconds.
%
%    real y0(4): [ phi1; dtphi1; phi2; dtphi2 ]
%    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. 
%  
%    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'.
%
  [ g, l1, l2, m1, m2 ] = double_pendulum_parameters ( );
%
%  Use ode45 to compute the solution.
%
  f = @ double_pendulum_deriv;
  t0 = tspan(1);
  tstop = tspan(2);

  sol = ode45 ( f, tspan, y0 );
%
%  Interpolate solution values at NFRAMES equally spaced times.
%
  nframes = tstop * fps;
  t = linspace ( t0, tstop, nframes );
  y = deval ( sol, t );
%
%  Extract solution data.
%
  phi1 = y(1,:)';
  dtphi1 = y(2,:)';
  phi2 = y(3,:)';
  dtphi2 = y(4,:)';
%
%  Open the movie file.
%
  if ( movie )
    filename = sprintf ( '%s.avi', label );
    vw = VideoWriter ( filename );
    open ( vw );
  end
%
%  Prepare the graphics.
% 
  clf ( );
  hold ( 'on' );
  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' );
  for i = 1 : 32
    xp = cos ( i * 2 * pi / 32 );
    yp = sin ( i * 2 * pi / 32 );
    plot ( l1 * xp, l1 * yp, ...
      'markersize', 10, 'marker', '.' );
    plot ( ( l1 + l2 ) * xp, ( l1 + l2 ) * yp, ...
      'markersize', 10, 'marker', '.' );
    plot ( abs ( l1 -l2 ) * xp, abs ( l1 - l2 ) * yp, ...
      'markersize', 10, 'marker', '.' );
  end
  hold ( 'off' )

  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 dydt = double_pendulum_deriv ( t, y )

%*****************************************************************************80
%
%% double_pendulum_deriv defines the double pendulum ODE.
%
%  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:
%
%    16 November 2020
%
%  Author:  
%
%    Original MATLAB version by Alexander Erlich (alexander.erlich@gmail.com)
%    Modifications by John Burkardt.
%
%  Input:
%
%    real t: the time. 
%
%    real y(4): the current solution.
%
%  Output:
%
%    real dydt(4): the right hand side of the ODEs.
%
  [ g, m1, m2, l1, l2 ] = double_pendulum_parameters ( );

  dydt = zeros ( 4, 1 );

  dydt(1) = y(2);

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

  dydt(3) = y(4);

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

  return
end
function [ g, m1, m2, l1, l2 ] = double_pendulum_parameters ( )

%*****************************************************************************80
%
%% double_pendulum_parameters: parameters for the double pendulum ODE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 March 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real G: the gravitational force coefficient.
%
%    real M1, M2: the masses of pendulums 1 and 2.
%
%    real L1, L2: the lengths of pendulums 1 and 2.
%
  g = 9.81;
  m1 = 1.0;
  m2 = 1.0;
  l1 = 2.0;
  l2 = 1.0;

  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


