function L = laplacian5_torus ( A, dx, dy )

%*****************************************************************************80
%
%% laplacian5_torus() uses a 5 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.
%
%    The computation is done using for loops.  It could be speeded up
%    by using vectorized MATLAB statements.
%
%    Because of the periodic boundary conditions, the first and last rows
%    and columns are identified, that is, equivalent.  Therefore, when
%    computing indices im1, ip1, jm1, jp1, we have to skip over the last
%    row or column.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    22 September 2024
%
%  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
%    5 point periodic stencil.
%
  [ nxp, nyp ] = size ( A );
  L = zeros ( nxp, nyp );

  for j = 1 : nyp

    if ( j == nyp - 1 )
      jp1 = 1;
    elseif ( j == nyp )
      jp1 = 2;
    else
      jp1 = j + 1;
    end

    if ( j == 1 )
      jm1 = nyp - 1;
    else
      jm1 = j - 1;
    end

    for i = 1 : nxp

      if ( i == nxp - 1 )
        ip1 = 1;
      elseif ( i == nxp )
        ip1 = 2;
      else
        ip1 = i + 1;
      end

      if ( i == 1 )
        im1 = nxp - 1;
      else
        im1 = i - 1;
      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

  return
end

