function D = divergence ( X, Y, U, V )

%*****************************************************************************80
%
%% divergence evaluates the divergence of a vector field.
%
%  Discussion:
%
%    We assume the grid spacings are uniform.
%
%    The value D(i,j) estimates du/dx + dv/dy at node (i,j).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(nr,nc), Y(nr,nc): the node coordinates.
%
%    real U(nr,nc), V(nr,nc): the node velocities.
%
%  Output:
%
%    real D(nr,nc): an estimate of the divergence of the velocity at each node.
%
  [ nr, nc ] = size ( U );

  dx = X(1,2) - X(1,1);
  dy = Y(1,1) - Y(2,1);
%
%  Estimate du/dx.
%  In column 1: use forward difference (column 2 is ahead)
%  In columns 2:nc-1, use centered difference (column c+1 - column c-1)
%  In column nc, use backward difference (column c-1 is behind)
%
  for c = 1 : nc
    if ( c == 1 )
      dUdX(1:nr,c) = ( U(1:nr,c+1) - U(1:nr,c) ) / dx;
    elseif ( c < nc )
      dUdX(1:nr,c) = ( U(1:nr,c+1) - U(1:nr,c-1) ) / dx / 2.0;
    else
      dUdX(1:nr,c) = ( U(1:nr,c) - U(1:nr,c-1) ) / dx;
    end
  end
%
%  Estimate dv/dy
%  In row 1: use backward difference (row r+1 is lower)
%  In rows 2:nr-1, use centered difference (row r-1 minus row r+1)
%  In row nr, use forward difference (row r-1 is higher)
%
  for r = 1 : nr
    if ( r == 1 )
      dVdY(r,1:nc) = ( V(r,1:nc) - V(r+1,1:nc) ) / dy;
    elseif ( r < nr )
      dVdY(r,1:nc) = ( V(r-1,1:nc) - V(r+1,1:nc) ) / dy / 2.0;
    else
      dVdY(r,1:nc) = ( V(r-1,1:nc) - V(r,1:nc) ) / dy;
    end
  end

  D = dUdX + dVdY;

  return
end

