function hw8 ( )

%*****************************************************************************80
%
%% hw8 plots a stream function psi from a velocity data file.
%
%  Discussion:
%
%    A fluid flow problem was solved on a 41x41 grid.
%
%    The file 'xy_inout.txt' lists pairs of (x,y) values.
%
%    The file 'uv_inout.txt' lists pairs of (u,v) values.
%
%    The ordering of the nodes is suggested by this array:
%
%      1641 1642 1643 1644 ... 1680 1681
%      .... .... .... .... ... .... ....
%        83   84   85   85 ...  122  123
%        42   43   44   45 ...   81   82
%         1    2    3    4 ...   40   41
%
%    In other words, the first node has the lowest values of X and Y.
%
%    This program reads the data, computes the stream function, and plots it.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 June 2020
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'hw8:\n' );
  fprintf ( 1, '  Read files of XY coordinates and UV velocities.\n' );
  fprintf ( 1, '  Compute the stream function, and plot it.\n' );
%
%  Read the XY coordinates, two vectors of 41*41 pairs of (x,y) values.
%
  xy = load ( 'xy_inout.txt' );
  x = xy(:,1);
  y = xy(:,2);
%
%  Reshape the vectors into arrays.
%
  X = reshape ( x, 41, 41 );
  Y = reshape ( y, 41, 41 );
%
%  Read the UV values, two vectors of 41*41 pairs of (u,v) values.
%
  uv = load ( 'uv_inout.txt' );
  u = uv(:,1);
  v = uv(:,2);
%
%  Reshape the vectors into arrays.
%
  U = reshape ( u, 41, 41 );
  V = reshape ( v, 41, 41 );
%
%  Compute the stream function.
%
  PSI = psi_from_uv ( X, Y, U, V );
%
%  Plot the stream function.
%
  figure ( );
  contourf ( X, Y, PSI );
  xlabel ( '<-- x -->' );
  ylabel ( '<-- y -->' );
  title ( 'INOUT Stream function \psi(x,y)' );
  
  filename = 'hw8.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  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 = 1 : nr

    if ( r == 1 )
      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

  return
end
