function spiral_pde_movie ( nx, ny, nt, dx, dy, dt, laplacian )

%*****************************************************************************80
%
%% spiral_pde_movie() creates a movie of a reaction-diffusion PDE solution.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 February 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer nx, ny: the number of nodes in the X and Y directions.
%
%    integer nt: the number of time steps.
%
%    real dx, dy: the spacing in the X and Y directions.
%
%    real dt: the size of the time step.
%
%    function handle laplacian: a MATLAB function of the form
%    L = laplacian ( U, X, Y )
%    which estimates the laplacian of the array U at (X,Y).
%    Initialize by a command like
%      laplacian = @ laplacian5_torus;
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'spiral_pde_movie:\n' );
%
%  Get parameters:
%
  [ alpha, beta, delta, epsilon ] = spiral_parameters ( );
%
%  Create the X and Y node arrays.
%
  [ X, Y ] = spiral_xy ( nx, ny );
%
%  Initialize the movie.
%
  colormap ( 'jet' );
  moviename = 'spiral_pde_movie.avi';
  Video = VideoWriter ( moviename, 'Motion JPEG AVI' );
  Video.FrameRate = 10;
  open ( Video );
%
%  Time loop using k as counter.
%
  for k = 0 : nt

    t = k * dt;
%
%  Set initial condition for U and V.
%
    if ( k == 0 )

      [ U, V ] = spiral_initial_condition ( X, Y, alpha );
%
%  Take a time step.
%
    else
%
%  Compute the laplacians of U and V using a user-specified function.
%
      del2U = laplacian ( U, dx, dy );
      del2V = laplacian ( V, dx, dy );
%
%  Take an explicit Euler time step for U and V.
%
      dUdt = del2U + ( 1.0 / epsilon ) * ...
        U .* ( 1.0 - U ) .* ( U - ( V + beta ) / alpha );
      dVdt = delta * del2V + U - V;

      U = U + dt * dUdt;
      V = V + dt * dVdt;
%
%  Apply the boundary conditions dU/dn = 0, dV/dn = 0.
%
      U(1:nx,1)  = U(1:nx,2);
      U(1:nx,ny) = U(1:nx,ny-1);
      U(1,1:ny)  = U(2,1:ny);
      U(nx,1:ny) = U(nx-1,1:ny);

      V(1:nx,1)  = V(1:nx,2);
      V(1:nx,ny) = V(1:nx,ny-1);
      V(1,1:ny)  = V(2,1:ny);
      V(nx,1:ny) = V(nx-1,1:ny);

    end
%
%  Plot once every 25 time steps.
%
    if ( mod ( k, 25 ) == 0 )

      figure ( 1 );
      clf ( );
      label = sprintf ( 'Spiral PDE, step = %d, time = %f\n', k, t );
      sgtitle ( label );

      subplot ( 1, 2, 1 );
      contourf ( X, Y, U, 'edgecolor', 'none' );
      axis ( 'equal' );
      set ( gca, 'xtick', [] );
      set ( gca, 'ytick', [] );
      title ( 'U(x,y)' );

      subplot ( 1, 2, 2 );
      contourf ( X, Y, V, 'edgecolor', 'none' );
      axis ( 'equal' );
      set ( gca, 'xtick', [] );
      set ( gca, 'ytick', [] );
      title ( 'V(x,y)' );

      frame = getframe ( 1 );
      writeVideo ( Video, frame );

    end

  end

  close ( Video );
  fprintf ( 1, '  Movie saved as "%s"\n', moviename );

  return
end

