function hw10 ( )

%*****************************************************************************80
%
%% hw10 solves the spiral reaction-diffusion equation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 June 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'hw10:\n' );
  fprintf ( 1, '  Solve the spiral reaction-diffusion equation.\n' );
%
%  Set N, the number of spatial nodes in each direction.
%
  n = 401;
%
%  Set xmin, xmax, vector x, spacing dx.
%
  nx = n;
  xmin = 0.0;
  xmax = 80.0;
  dx = ( xmax - xmin ) / ( nx - 1 );
%
%  Set ymin, ymax, vector y, spacing dy.
%
  ny = n;
  ymin = 0.0;
  ymax = 80.0;
  dy = ( ymax - ymin ) / ( ny - 1 );
%
%  Set number of steps nt, tmin, tmax, t, dt
%
  nt = 5001;
  tmin = 0.0;
  tmax = 5.0;
  t = linspace ( tmin, tmax, nt );
  dt = ( tmax - tmin ) / ( nt - 1 );
%
%  Set problem parameters alpha, beta, delta, epsilon.
%
  alpha =   0.25;
  beta =    0.001;
  delta =   0.00001;
  epsilon = 0.002;
%
%  Create the X and Y node arrays.
%
  [ X, Y ] = xy_spiral ( nx, ny );
%
%  Time loop using k as counter.
%
  for k = 1 : nt
%
%  Set initial condition for U and V.
%
    if ( k == 1 )

      [ U, V ] = uv_spiral ( X, Y, alpha );
%
%  Take a time step.
%
    else
%
%  Compute the laplacians of U and V, using either a 5 or 9 point stencil.
%
      if ( false )
        del2U = laplacian5_torus ( U, dx, dy );
        del2V = laplacian5_torus ( V, dx, dy );
      else
        del2U = laplacian9_torus ( U, dx, dy );
        del2V = laplacian9_torus ( V, dx, dy );
      end
%
%  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 ) == 1 )

      figure ( 1 )
      label = sprintf ( 'Spiral, step = %d, time = %f\n', k, t(k) );
      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)' );

    end

  end
%
%  Save the last plot.
%
  filename = 'hw10.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

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

%*****************************************************************************80
%
%% laplacian5_torus uses a 5 point Laplacian stencil on a torus.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2020
%
%  Author:
%
%    John Burkardt
%

%
%  Vectorized version:
%
  if ( true )

    f_xx = (        circshift ( A, [  1,  0 ] ) ... 
           - 2.0 .* circshift ( A, [  0,  0 ] ) ...
                  + circshift ( A, [ -1,  0 ] ) ) / dx^2;

    f_yy = (        circshift ( A, [  0,  1 ] ) ...
           - 2.0 .* circshift ( A, [  0,  0 ] ) ...
                  + circshift ( A, [  0, -1 ] ) ) / dy^2;

    L = f_xx + f_yy;
%
%  Scalar version.
%
  else

    [ nxp, nyp ] = size ( A );
    L = zeros ( nxp, nyp );

    for j = 1 : nyp

      jp1 = j + 1;
      if ( nyp < jp1 )
        jp1 = 1;
      end
      jm1 = j - 1;
      if ( jm1 < 1 )
        jm1 = nyp;
      end

      for i = 1 : nxp

        ip1 = i + 1;
        if ( nxp < ip1 )
          ip1 = 1;
        end
        im1 = i - 1;
        if ( im1 < 1 )
          im1 = nxp;
        end

        L(i,j) = ( A(ip1,j) - 2.0 * A(i,j) + A(im1,j) ) / dx^2 ...
               + ( A(i,jp1) - 2.0 * A(i,j) + A(i,jm1) ) / dy^2;
      end
    end

  end

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

%*****************************************************************************80
%
%% laplacian9_torus uses a 9 point Laplacian stencil on a torus.
%
%  Discussion:
%
%    1) Assumes uniform, equal mesh in both directions.
%    2) Assumes the spacing is dx in both directions.
%
%  Modified:
%
%    28 June 2020
%
  [ nxp, nyp ] = size ( A );
  L = zeros ( size ( A ) );

  for j = 1 : nyp

    jp1 = j + 1;
    if ( nyp < jp1 )
      jp1 = 1;
    end
    jm1 = j - 1;
    if ( jm1 < 1 )
      jm1 = nyp;
    end

    for i = 1 : nxp

      ip1 = i + 1;
      if ( nxp < ip1 )
        ip1 = 1;
      end
      im1 = i - 1;
      if ( im1 < 1 )
        im1 = nxp;
      end

      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 [ X, Y ] = xy_spiral ( nx, ny )

%*****************************************************************************80
%
%% xy_spiral: node coordinates for the spiral reaction-diffusion equation.
%
%  Modified:
%
%    26 June 2020
%
%  Input:
%
%    integer nx, ny: the number of rows and columns of nodes.
%
%  Output:
%
%    real X(nx,ny), Y(nx,ny): the node coordinates.
%
  x = linspace ( 0.0, +80.0, nx );
  y = linspace ( 0.0, +80.0, ny );

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

  return
end
function [ U, V ] = uv_spiral ( X, Y, alpha )

%*****************************************************************************80
%
%% uv_spiral: initial velocity vectors for spiral reaction-diffusion equation.
%
%  Modified:
%
%    26 June 2020
%
%  Input:
%
%    real X(nr,nc), Y(nr,nc): the node coordinates.
%
%    real alpha: a problem parameter.
%
%  Output:
%
%    real U(nr,nc), V(nr,nc): the node velocities.
%
  U = ( 40.0 <= X );

  V = 0.5 * alpha * ( 40.0 <= Y );

  return
end

