function flow_vector ( )

%*****************************************************************************80
%
%% flow_vector() makes a vector plot of a flow velocity field.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2022
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'flow_vector():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s:\n', version ( ) );
  fprintf ( 1, '  Make a vector plot of a velocity flow field.\n' );
%
%  Get the data.
%
  n = 16;
  [ xmat, ymat, umat, vmat ] = flow_field ( n );
%
%  Get ready to plot.
%
  figure ( );
  clf ( );
%
%  Plot the vectors.
%
  scale = 1.0;
  quiver ( xmat, ymat, umat, vmat, scale );
%
%  Finish the plot.
%
  axis ( 'equal' );
  axis ( 'tight' );
  xlabel ( '<--X-->' );
  ylabel ( '<--Y-->' );
  title ( 'Velocity flow field' );
  filename = 'flow_vector.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'flow_vector():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function [ xmat, ymat, umat, vmat ] = flow_field ( n )

%*****************************************************************************80
%
%% flow_field() returns velocity data for a flow.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    24 May 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real N: the number of nodes in the X and Y directions.
%
%  Output:
%
%    real XMAT(N,N), YMAT(N,N): the node coordinates.
%
%    real UMAT(N,N), VMAT(N,N): the velocity components.
%
  xvec = linspace ( 0.0, 1.0, n );
  yvec = linspace ( 0.0, 1.0, n );

  [ xmat, ymat ] = meshgrid ( xvec, yvec );

  umat = - (       xmat.^4 - 2.0 * xmat.^3 + xmat.^2 ) ...
        .* ( 2.0 * ymat.^3 - 3.0 * ymat.^2 + ymat    );

  vmat =   ( 2.0 * xmat.^3 - 3.0 * xmat.^2 + xmat    ) ...
        .* (       ymat.^4 - 2.0 * ymat.^3 + ymat.^2 );

  return
end
