function pendulum_double_movie ( tspan, y0, fps, movie, label )

%*****************************************************************************80
%
%% pendulum_double_movie animates a double pendulum solution.
%
%  Discussion:
%
%    Example direct call:
%      pendulum_double_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 ] = pendulum_double_parameters ( );
%
%  Use ode45 to compute the solution.
%
  f = @ pendulum_double_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 );
    V = VideoWriter ( filename );
    V.FrameRate = 30;
    open ( V );
  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 ( V, frame );
    if ( ~ movie )
      pause ( t(i+1) - t(i) );
    end

  end

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

  return
end
