function vector_stream_grid ( node_file, vector_file, grid_num, sample_num )

%% VECTOR_STREAM_GRID contour plots vector streamlines from scattered data.
%
%  Discussion:
%
%    The prime feature of this program is that it uses interpolation
%    from possibly scattered data to get a uniform grid of data,
%    from which the MATLAB STREAMLINE command can do its work.
%
%    This program was written as a quick and convenient way to
%    view finite element data.  It can also be used for any situation
%    in which a collection of points and velocities are available.
%
%    However, a small "sin" is committed when we use this program
%    on finite element data, since we take only the node locations 
%    and velocities, but not the elements.  We simply use MATLAB's grid 
%    data feature which constructs a cubic spline interpolant to 
%    scattered data.  This is OK, and probably doesn't distort the 
%    data too much, but in fact, we are not actually viewing 
%    values of the finite element function, and there is no guarantee
%    that this surrogate function will actually satisfy the boundary
%    conditions or other characteristics of the finite element
%    function.
%
%  Modified:
%
%    02 October 2006
%
%  Author:
%
%    John Burkardt
%    Hyung-Chun Lee
%
%  Usage:
%
%    vector_stream_grid ( node_file, vector_file, contour_num, sample_num )
%
%  Parameters:
%
%    Input, string NODE_FILE, the name of the node file.
%
%    Input, string VECTOR_FILE, the name of the vector file.
%
%    Input, integer GRID_NUM, the number of grid points to use
%    when generating the interpolating function.  The grid points
%    form a 2D grid. GRID_NUM is actually the number of such points
%    along one dimension.
%
%    Input, integer SAMPLE_NUM, the number of sample points to use
%    when drawing streamlines.  The sample points
%    form a 2D grid. SAMPLE_NUM is actually the number of such points
%    along one dimension.
%
  fprintf ( 1, '\n' );
  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'VECTOR_STREAM_GRID:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Display a plot of vector streamlines,\n' );
  fprintf ( 1, '  based on scattered data.\n' );
%
%  First argument is the node file.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    node_file = input ( '  Enter the name of the node file.' );
  end
%
%  Second argument is the vector file.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    vector_file = input ( '  Enter the name of the vector file.' );
  end
%
%  Third argument is the number of grid points.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    grid_num = input ( '  Enter the number of grid points in one dimension.' );
  end
%
%  Fourth argument is the number of sample points.
%
  if ( nargin < 4 )
    fprintf ( 1, '\n' );
    sample_num = input ( '  Enter the number of sample points in one dimension.' );
  end
%
%  Read the node data.
%
  [ dim_num, node_num ] = dtable_header_read (  node_file );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', node_file );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension DIM_NUM = %d\n', dim_num );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  if ( dim_num ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VECTOR_STREAM_GRID - Fatal error!\n' );
    fprintf ( 1, '  Node data must have spatial dimension 2.\n' );
    error ( 'VECTOR_STREAM_GRID - Fatal error!' );
  end

  node_xy = dtable_data_read ( node_file, dim_num, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', node_file );

  r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, 2, 5, ...
    '  2 by 5 portion of data read from file:' );
%
%  Read the vector data.
%
  [ vector_order, node_num2 ] = dtable_header_read ( vector_file );

  if ( vector_order ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VECTOR_STREAM_GRID - Fatal error!\n' );
    fprintf ( 1, '  The vector file must list pairs of numbers.\n' );
    error ( 'VECTOR_STREAM_GRID - Fatal error!' );
  end

  if ( node_num2 ~= node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VECTOR_STREAM_GRID - Fatal error!\n' );
    fprintf ( 1, '  The number of vectors and nodes do not match.\n' );
    error ( 'VECTOR_STREAM_GRID - Fatal error!' );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', vector_file );

  uv = dtable_data_read ( vector_file, vector_order, node_num );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the data in "%s".\n', vector_file );

  r8mat_transpose_print_some ( 2, node_num, uv, ...
    1, 1, 2, 10, '  Portion of vector array:' );
%
%  Determine the rectangle that contains the data.
%
  xlo = min ( node_xy(1,:) );
  xhi = max ( node_xy(1,:) );
  ylo = min ( node_xy(2,:) );
  yhi = max ( node_xy(2,:) );
%
%  Determine the coordinates of the grid points.
%
  if ( grid_num <= 1 )

    x_vec = 0.5 * ( xlo + xhi );
    y_vec = 0.5 * ( ylo + yhi );

  elseif ( yhi - ylo < xhi - xlo ) then

    grid_num_x = grid_num;

    inc = ( xhi - xlo ) / ( grid_num_x - 1 );

    grid_num_y = ( round ( yhi - ylo ) / inc );
    yfrac = 0.5 * ( ( yhi - ylo ) - grid_num_y * inc );

    x_vec = xlo         : inc : xhi;
    y_vec = ylo + yfrac : inc : yhi - yfrac;

  else

    grid_num_y = grid_num;

    inc = ( yhi - ylo ) / ( grid_num_y - 1 );

    grid_num_x = ( round ( xhi - xlo ) / inc );
    xfrac = 0.5 * ( ( xhi - xlo ) - grid_num_x * inc );

    x_vec = xlo + xfrac : inc : xhi - xfrac;
    y_vec = ylo         : inc : yhi;

  end
%
%  Create the grid points from the X and Y coordinates.
%  Create the U and V grid values by interpolation of the scattered data.
%
  [ x_grid, y_grid ] = meshgrid ( x_vec, y_vec );

  u_grid = griddata ( node_xy(1,:), node_xy(2,:), uv(1,:), ...
    x_grid, y_grid, 'cubic' );

  v_grid = griddata ( node_xy(1,:), node_xy(2,:), uv(2,:), ...
    x_grid, y_grid, 'cubic' );
%
%  Determine the coordinates of the sample points.
%
  if ( sample_num <= 1 )

    x_vec = 0.5 * ( xlo + xhi );
    y_vec = 0.5 * ( ylo + yhi );

  elseif ( yhi - ylo < xhi - xlo ) then

    sample_num_x = sample_num;

    inc = ( xhi - xlo ) / ( sample_num_x - 1 );

    sample_num_y = ( round ( yhi - ylo ) / inc );
    yfrac = 0.5 * ( ( yhi - ylo ) - sample_num_y * inc );

    x_vec = xlo         : inc : xhi;
    y_vec = ylo + yfrac : inc : yhi - yfrac;

  else

    sample_num_y = sample_num;

    inc = ( yhi - ylo ) / ( sample_num_y - 1 );

    sample_num_x = ( round ( xhi - xlo ) / inc );
    xfrac = 0.5 * ( ( xhi - xlo ) - sample_num_x * inc );

    x_vec = xlo + xfrac : inc : xhi - xfrac;
    y_vec = ylo         : inc : yhi;

  end
  
  [ x_sample, y_sample ] = meshgrid ( x_vec, y_vec );
%
%  Make a vector streamline contour plot using the grid data.
%
  streamline ( x_grid, y_grid, u_grid, v_grid, x_sample, y_sample, ...
    [ 0.1, 100 ] )

  xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16, 'Rotation', 0 );

  title ( 'Vector Streamlines', 'FontName', 'Helvetica', 'FontWeight', ...
    'bold', 'FontSize', 16 );

  axis square

  fprintf ( 1, '\n' );
  fprintf ( 1, 'VECTOR_STREAM_GRID:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp;