function vector_arrows_grid ( node_file, velocity_file, scale, number )

%% VELOCITY_ARROWS_GRID plots a velocity field grid based on scattered data.
%
%  Discussion:
%
%    The prime feature of this program is that it displays the data
%    on a uniform grid of specified density, even though the original 
%    data is scattered or of a different density.
%
%    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:
%
%    29 September 2006
%
%  Author:
%
%    John Burkardt
%    Hyung-Chun Lee
%
%  Usage:
%
%    velocity_arrows_grid ( node_file, velocity_file, scale, number )
%
%  Parameters:
%
%    Input, string NODE_FILE, the name of the node file.
%
%    Input, string VELOCITY_FILE, the name of the velocity file.
%
%    Input, real SCALE, a scale factor.  If SCALE = 0, MATLAB
%    will do the scaling automatically.
%
%    Input, integer NUMBER, the number of grid points to use along
%    the widest direction.  (If the region is square, an
%    NUMBERxNUMBER grid will be used, otherwise the smaller spatial
%    dimension will have fewer grid points.)
%
  fprintf ( 1, '\n' );
  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'VELOCITY_ARROWS_GRID:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Display a grid of vectors 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 velocity file.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    velocity_file = input ( '  Enter the name of the velocity file.' );
  end
%
%  Third argument is the scale.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    scale = input ( '  Enter the vector scale factor or 0 for autoscaling.' );
  end
%
%  Fourth argument is the number of grid points.
%
  if ( nargin < 4 )
    fprintf ( 1, '\n' );
    number = input ( '  Enter the number of grid 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, 'VELOCITY_ARROWS_GRID - Fatal error!\n' );
    fprintf ( 1, '  Node data must have spatial dimension 2.\n' );
    error ( 'VELOCITY_ARROWS_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 velocity data.
%
  [ velocity_order, node_num2 ] = dtable_header_read ( velocity_file );

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

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

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

  uv = dtable_data_read ( velocity_file, velocity_order, node_num );

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

  r8mat_transpose_print_some ( 2, node_num, uv, ...
    1, 1, 2, 10, '  Portion of velocity 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 ( number <= 1 )

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

  elseif ( yhi - ylo < xhi - xlo ) then

    number_x = number;

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

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

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

  else

    number_y = number;

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

    number_x = ( round ( xhi - xlo ) / inc );
    xfrac = 0.5 * ( ( xhi - xlo ) - number_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' );
%
%  Make a vector plot using the grid data.
%
  quiver ( x_grid, y_grid, u_grid, v_grid, scale )

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

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

  title ( 'Velocity Vectors', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
    'FontSize', 16 );

  axis square

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

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