function [ U, V ] = uv_from_psi ( X, Y, PSI )

%*****************************************************************************80
%
%% uv_from_psi computes a velocity field (u,v) from a stream function.
%
%  Discussion:
%
%    We assume the x and y grid spacings are each uniform.
%
%    U =   d PSI/dy
%    V = - d PSI/dx
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(nr,nc), Y(nc,nc): the node coordinates.
%
%    real PSI(nr,nc): the stream function.
%
%  Output:
%
%    real U(nr,nc), V(nr,nc): the velocity.
%
  [ nr, nc ] = size ( PSI );

  dx = X(1,2) - X(1,1);
  dy = Y(1,1) - Y(2,1);
  
  U = zeros ( nr, nc );
  V = zeros ( nr, nc );
%
%  For U, use centered differences except in first and last rows.
%  For Y derivatives, the larger Y value is at the SMALLER r value.
%
  for r = 1 : nr
    if ( r == 1 )
      U(r,1:nc) = ( PSI(r,1:nc) - PSI(r+1,1:nc) ) / dy;
    elseif ( r < nr )
      U(r,1:nc) = ( PSI(r-1,1:nc) - PSI(r+1,1:nc) ) / dy / 2.0;
    else
      U(r,1:nc) = ( PSI(r-1,1:nc) - PSI(r,1:nc) ) / dy;
    end
  end
%
%  For V, use centered differences except in first and last columns.
%  For X derivatives, the larger X value is at the larger c value.
%
  for c = 1 : nc
    if ( c == 1 )
      V(1:nr,c) = - ( PSI(1:nr,c+1) - PSI(1:nr,c) ) / dx;
    elseif ( c < nc )
      V(1:nr,c) = - ( PSI(1:nr,c+1) - PSI(1:nr,c-1) ) / dx / 2.0;
    else
      V(1:nr,c) = - ( PSI(1:nr,c) - PSI(1:nr,c-1) ) / dx;
    end
  end

  return
end

