function lp = laplacian3_circle ( up, dx )

%*****************************************************************************80
%
%% laplacian3_circle() approximates the laplacian on a periodic 1D domain.
%
%  Discussion:
%
%    If the domain is represented by nx equally spaced points, it is
%    assumed that only the data for the first nxp=nx-1 points is supplied.
%    In other words, the last entry of up() is at the node just BEFORE
%    the node where the periodic condition takes effect.
%
%    This is the standard laplacian estimate on a uniform grid, 
%    except that, by periodicity, and suppressing the last node, we have
%    the east neighbor of node 1 is the last node, and
%    the west neighbor of the last node is node 1.
%
%    This function offers four equivalent ways of computing the Laplacian.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    27 June 2020
%
%  Input:
%
%    real up(nxp): the values of a function on an equally spaced grid
%    of points in a periodic domain.
%
%    real dx: the spacing between any two consecutive nodes.
%
%  Output:
%
%    real lp(nxp): the estimate for the laplacian.
%
  nxp = length ( up );
  lp = zeros ( size ( up ) );
%
%  Simplest formula.
%
  if ( false )

    ip = 1;
    lp(ip) = ( up(ip+1) - 2.0 * up(ip) + up(nxp) ) / dx^2;

    for ip = 2 : nxp - 1
      lp(ip) = ( up(ip+1) - 2.0 * up(ip) + up(ip-1) ) / dx^2;
    end

    ip = nxp; 
    lp(ip) = ( up(1) - 2.0 * up(ip) + up(ip-1) ) / dx^2;
%
%  Use the mod() function to take care of indexing.
%
  elseif ( false )

    for ip = 1 : nxp
      east = 1 + mod ( ip,     nxp );
      west = 1 + mod ( ip - 2, nxp );
      lp(ip) = ( up(east) - 2.0 * up(ip) + up(west) ) / dx^2;
    end
%
%  Use "circshift" to force indices to wrap around.
%
  elseif ( true )
    lp = (         circshift ( up,  1 ) ... 
           - 2.0 * circshift ( up,  0 ) ...
                 + circshift ( up, -1 ) ) / dx^2;
%
%  Just make a list of the indexes.
%
  else

    im1 = [ nxp, 1:nxp-2, nxp-1 ];
    i =   [ 1,   2:nxp-1, nxp   ];
    ip1 = [ 2,   3:nxp,   1     ];
    lp = ( up(ip1) - 2.0 * up(i) + up(im1) ) / dx^2;

  end

  return
end

