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 MIT 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

