function spiral_pde_test ( )

%*****************************************************************************80
%
%% spiral_pde_test() tests spiral_pde().
%
%  Discussion:
%
%    I had to modify the original MATLAB code in order to run under Octave.
%    The sgtitle() function was not available.
%    The problem size was reduced from 401x401 to 201x201.
%    The contourf() function seemed to have trouble with the initial data,
%    which was all 0 and 1 values.  I decided to randomly perturb the data
%    a bit, which seemed to help.
%
%    After all these odd problems, I got the code to run to the end and
%    produce the spiral solution.  I don't know whether I encountered serious 
%    or trivial issues in Octave.  This doesn't happen often, but I pity a
%    user who needs a result and can't think of an effective workaround.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    05 February 2021
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../spiral_pde' );

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'spiral_pde_test():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  spiral_pde() solves a reaction-diffusion\n' );
  fprintf ( 1, '  partial differential equation with spiral solution.\n' );
%
%  Set N, the number of spatial nodes in each direction.
%  We need to reduce from n = 401, because Octave can't handle big sizes.
%
% n = 401;
  n = 201;
%
%  Set the X values.
%
  nx = n;
  xmin = 0.0;
  xmax = 80.0;
  dx = ( xmax - xmin ) / ( nx - 1 );
  x = linspace ( xmin, xmax, n );
%
%  Set the Y values.
%
  ny = n;
  ymin = 0.0;
  ymax = 80.0;
  dy = ( ymax - ymin ) / ( ny - 1 );
  y = linspace ( ymin, ymax, n );
%
%  Set the X, Y grid.
%
  [ X, Y ] = meshgrid ( x, y );
%
%  Set the T values.
%
  nt = 5001;
  tmin = 0.0;
  tmax = 5.0;
  dt = ( tmax - tmin ) / ( nt - 1 );
  t = linspace ( tmin, tmax, nt );
%
%  Report discretization values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  discretization values:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    nx = %d\n', nx );
  fprintf ( 1, '    xmin = %g\n', xmin );
  fprintf ( 1, '    xmax = %g\n', xmax );
  fprintf ( 1, '    dx = %g\n', dx );
  fprintf ( 1, '    ny = %d\n', ny );
  fprintf ( 1, '    ymin = %g\n', ymin );
  fprintf ( 1, '    ymax = %g\n', ymax );
  fprintf ( 1, '    dy = %g\n', dy );
  fprintf ( 1, '    nt = %d\n', nx );
  fprintf ( 1, '    tmin = %g\n', tmin );
  fprintf ( 1, '    tmax = %g\n', tmax );
  fprintf ( 1, '    dt = %g\n', dy );
%
%  Set problem parameters alpha, beta, delta, epsilon.
%
  [ alpha, beta, delta, epsilon ] = spiral_parameters ( );
%
%  Report parameter values.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  parameters:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    alpha = %g\n', alpha );
  fprintf ( 1, '    beta = %g\n', beta );
  fprintf ( 1, '    delta = %g\n', delta );
  fprintf ( 1, '    epsilon = %g\n', epsilon );
%
%  Desperate remedy.
%  Octave contourf() seems confused if data consists of only a few
%  distinct values.
%
  R = 0.95 + 0.05 * rand ( 1:nx, 1:ny );

  for k = 1 : nt
%
%  Set initial condition for U and V.
%
    if ( k == 1 )

      [ U, V ] = spiral_initial_condition ( X, Y );
%
%  Use the Forward Euler method to take a time step.
%
    else

      [ dUdt, dVdt ] = spiral_deriv ( t, U, V, dx, dy );
      U = U + dt * dUdt;
      V = V + dt * dVdt;
%
%  Apply zero Neumann boundary conditions.
%
      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 - 1, 25 ) == 0 )

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

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

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

      if ( mod ( k - 1, 500 ) == 0 )
        filename = sprintf ( 'spiral_%04d.png', k );
        print ( '-dpng', filename );
        fprintf ( 1, '  Graphics saved as "%s"\n', filename );
      end

    end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'spiral_pde_test():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../spiral_pde' );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp() prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

