function L = laplacian9_torus ( A, dx, dy )

%*****************************************************************************80
%
%% laplacian9_torus uses a 9 point Laplacian stencil on a torus.
%
%  Discussion:
%
%    1) Assumes the region is a rectangular torus;
%    2) Assumes the mesh is uniform in both X and Y directions.
%    3) Assumes the spacing is dx in both directions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 December 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(NXP,NYP): function values, sampled on a regular grid.
%
%    real DX, DY: the grid spacing in the X and Y directions.
%
%  Output:
%
%    real L(NXP,NYP): the values of the Laplacian, as estimated by the
%    9 point periodic stencil.
%
  [ 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

