function exercise4 ( )

%*****************************************************************************80
%
%% exercise4 determines a stream function psi from a velocity field (u,v).
%
%  Discussion:
%
%    We assume the x and y grid spacings are each uniform.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X(nr,nc), Y(nc,nc): the node coordinates.
%
%    real U(nr,nc), V(nr,nc): horizontal and vertical velocities.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'exercise4\n' );
  fprintf ( 1, '  For each example flow, determine the stream function PSI\n' );
  fprintf ( 1, '  by integrating the velocity field (U,V).\n' );

  [ X, Y ] = xy_channel ( 5, 11 );
  [ U, V ] = uv_channel ( X, Y );
  PSI = psi_from_uv ( X, Y, U, V );

  [ X, Y ] = xy_corner ( 11, 11 );
  [ U, V ] = uv_corner ( X, Y );
  PSI = psi_from_uv ( X, Y, U, V );

  [ X, Y ] = xy_shear ( 11, 11 );
  [ U, V ] = uv_shear ( X, Y );
  PSI = psi_from_uv ( X, Y, U, V );

  [ X, Y ] = xy_vortex ( 11, 11 );
  [ U, V ] = uv_vortex ( X, Y );
  PSI = psi_from_uv ( X, Y, U, V );

  return
end
function PSI = psi_from_uv ( X, Y, U, V )

%*****************************************************************************80
%
%% psi_from_uv computes a stream function psi from a velocity field (u,v).
%
%  Discussion:
%
%    We assume the x and y grid spacings are each uniform.
%
%  Modified:
%
%    13 June 2020
%
%  Input:
%
%    real X(nr,nc), Y(nc,nc): the node coordinates.
%
%    real U(nr,nc), V(nr,nc): horizontal and vertical velocities.
%
%  Output:
%
%    real PSI(nr,nc): the stream function.
%
  [ nr, nc ] = size ( U );

  dx = X(1,2) - X(1,1);
  dy = Y(1,1) - Y(2,1);
  
  PSI = zeros ( nr, nc );

  for r = nr : -1 : 1

    if ( r == nr )
      PSI(r,1) = 0.0;
    else
      PSI(r,1) = PSI(r+1,1) + 0.5 * ( U(r+1,1) + U(r,1) ) * dy;
    end

    for c = 2 : nc
      PSI(r,c) = PSI(r,c-1) - 0.5 * ( V(r,c-1) + V(r,c) ) * dx;
    end

  end

  figure ( );
  contourf ( X, Y, PSI );
  xlabel ( '<-- x -->' );
  ylabel ( '<-- y -->' );
  title ( 'Stream function \psi(x,y)' );

  return
end
