function gray_scott_pde ( nx, ny, nt, tmax )

%*****************************************************************************80
%
%% gray_scott_pde() solves the Gray-Scott PDE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2021
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Willem Hundsdorfer, Jan Verwer,
%    Numerical solution of time-dependent advection-diffusion-reaction equations,
%    Springer, 2003
%    ISBN: 978-3-662-09017-6
%
%  Input:
%
%    integer NX, NY, the number of nodes in the X and Y directions.
%
%    integer NT, the number of time steps.
%
%    real TMAX, the final time.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'gray_scott_pde:\n' );
  fprintf ( 1, '  Solve the Gray-Scott diffusion-reaction PDE\n' );

% nx = 251;
  xmin = 0.0;
  xmax = 2.5;
  x = linspace ( xmin, xmax, nx );
  dx = ( xmax - xmin ) / ( nx - 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', nx );
  fprintf ( 1, '  %g <= X <= %g\n', xmin, xmax );
  fprintf ( 1, '  dx = %g\n', dx );

% ny = 251;
  ymin = 0.0;
  ymax = 2.5;
  y = linspace ( ymin, ymax, ny );
  dy = ( ymax - ymin ) / ( ny - 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of Y nodes = %d\n', ny );
  fprintf ( 1, '  %g <= Y <= %g\n', ymin, ymax );
  fprintf ( 1, '  dy = %g\n', dy );

  [ X, Y ] = meshgrid ( x, y );

% nt = 150001;
  tmin = 0.0;
% tmax = 1500.0;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of T values = %d\n', nt );
  fprintf ( 1, '  %g <= T <= %g\n', tmin, tmax );
  fprintf ( 1, '  dt = %g\n', dt );
%
%  Physical parameters:
%
  d1 = 8.0E-05;
  d2 = 4.0E-05;
  gamma = 0.024;
  kappa = 0.06;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  U diffusion D1 = %g\n', d1 );
  fprintf ( 1, '  V diffusion D2 = %g\n', d2 );
  fprintf ( 1, '  Flow rate Gamma = %g\n', gamma );
  fprintf ( 1, '  Kill rate Kappa = %g\n', kappa );
%
%  Begin the time loop.
%
  for it = 1 : nt
%
%  On first step, use the initial condition.
%
    if ( it == 1 )

      V = 0.25 * ( sin ( 4.0 * pi * X ) ).^2 .* ( sin ( 4.0 * pi * Y ) ).^2;
      V = V .* ( 1.0 <= X ) .* ( X <= 1.5 ) .* ( 1.0 <= Y ) .* ( Y <= 1.5 );
      U = ones ( nx, ny ) - 2.0 * V;
%
%  On subsequent steps, use the Euler method to update the data.
%
    else

      delU = laplacian9_torus ( U, dx, dy );
      delV = laplacian9_torus ( V, dx, dy );

      dUdt = d1 * delU - U .* V.^2 + gamma * ( 1.0 - U );
      dVdt = d2 * delV + U .* V.^2 - ( gamma + kappa ) * V;

      U = U + dt * dUdt;
      V = V + dt * dVdt;

    end
%
%  Display the data at certain steps.
%
    if ( mod ( it - 1, 100 ) == 0 )
      figure ( 1 );
      clf ( );
      image ( U, 'CDataMapping', 'scaled' );
      label = sprintf ( 'Gray-Scott PDE, U, time = %g\n', t(it) );
      title ( label );
      axis ( 'equal' );
      axis ( 'tight' );
      set ( gca, 'xtick', [] );
      set ( gca, 'ytick', [] );
    end
%
%  Save images at times matching the textbook results.
%
    if ( it - 1 == 0 )
      filename = 'gray_scott_u0000.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( it - 1 == 10000 )
      filename = 'gray_scott_u0100.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( it - 1 == 20000 )
      filename = 'gray_scott_u0200.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( it - 1 == 50000 )
      filename = 'gray_scott_u0500.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( it - 1 == 100000 )
      filename = 'gray_scott_u1000.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    elseif ( it - 1 == 150000 )
      filename = 'gray_scott_u1500.png';
      print ( '-dpng', filename );
      fprintf ( 1, '  Graphics saved as "%s"\n', filename );
    end

  end

  return
end
function L = laplacian9_torus ( A, dx, dy )

%*****************************************************************************80
%
%% laplacian9_torus() uses a 9 point Laplacian stencil on a rectangular grid.
%
%  Discussion:
%
%    1) Assumes uniform, equal mesh in both directions.
%    2) Assumes the spacing is dx in both directions.
%    3) Assumes periodic boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2020
%
%  Input:
%
%    real A(nxp,nyp): the data array.
%
%    real dx, dy: the X and Y spacing between data locations.
%
%  Output:
%
%    real L(nxp,nyp): the 9 point estimate for the Laplacian.
%
  [ nxp, nyp ] = size ( A );
  L = zeros ( nxp, nyp );

  for j = 1 : nyp

    jm1 = i4_wrap ( j - 1, 1, nyp );
    jp1 = i4_wrap ( j + 1, 1, nyp );

    for i = 1 : nxp

      im1 = i4_wrap ( i - 1, 1, nxp );
      ip1 = i4_wrap ( i + 1, 1, nxp );

      L(i,j) = ...
        ( 1.0 * A(im1,jm1) +  4.0 * A(im1,j) + 1.0 * A(im1,jp1) ...
        + 4.0 * A(i,  jm1) - 20.0 * A(i,j)   + 4.0 * A(i,  jp1) ...
        + 1.0 * A(ip1,jm1) +  4.0 * A(ip1,j) + 1.0 * A(ip1,jp1) ) ...
        / 6.0 / dx^2;

    end
  end

  return
end
function value = i4_wrap ( value, lo, hi )

%*****************************************************************************80
%
%% i4_wrap forces an integer to lie between given limits by wrapping.
%
%  Example:
%
%    LO = 4, HI = 8
%
%    In   Out
%
%    -2     8
%    -1     4
%     0     5
%     1     6
%     2     7
%     3     8
%     4     4
%     5     5
%     6     6
%     7     7
%     8     8
%     9     4
%    10     5
%    11     6
%    12     7
%    13     8
%    14     4
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer VALUE, an integer value.
%
%    integer LO, HI, the desired bounds for the integer value.
%
%  Output:
%
%    integer VALUE, a "wrapped" version of IVAL.
%
  value = lo + mod ( value - lo, hi + 1 - lo );

  return
end

