function vector_arrows_grid2 ( node_file, velocity_file, scale, number )

%% VELOCITY_ARROWS_GRID2 plots a velocity field grid based on scattered data.
%
%  Discussion:
%
%    VELOCITY_ARROWS_GRID2 is a MATLAB routine which reads 
%    a single file listing node coordinates and velocities,
%    possibly at scattered points, or on a grid that is too dense
%    or not dense enough; it creates a vector plot
%    of the velocity field on a uniform grid of user-specified density.
%
%    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:
%
%    04 June 2007
%
%  Author:
%
%    John Burkardt
%    Hyung-Chun Lee
%
%  Usage:
%
%    velocity_arrows_grid2 ( xyuv_file, scale, number )
%
%  Parameters:
%
%    Input, string XYUV_FILE, the name of the node and 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_GRID2:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Display a grid of vectors based on scattered data.\n' );
%
%  First argument is the xyuv file.
%
  if ( nargin < 1 )
    fprintf ( 1, '\n' );
    xyuv_file = input ( '  Enter the name of the node and velocity file.' );
  end
%
%  Second argument is the scale.
%
  if ( nargin < 2 )
    fprintf ( 1, '\n' );
    scale = input ( '  Enter the vector scale factor or 0 for autoscaling.' );
  end
%
%  Third argument is the number of grid points.
%
  if ( nargin < 3 )
    fprintf ( 1, '\n' );
    number = input ( '  Enter the number of grid points in one dimension.' );
  end
%
%  Read the xyuv data.
%
  [ dim_two, node_num ] = dtable_header_read (  xyuv_file );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read the header of "%s".\n', xyuv_file );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Data per row =              %d\n', dim_two );
  fprintf ( 1, '  Number of points NODE_NUM = %d\n', node_num );

  if ( dim_two ~= 4 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VELOCITY_ARROWS_GRID2 - Fatal error!\n' );
    fprintf ( 1, '  XYUV file must have spatial dimension 4.\n' );
    error ( 'VELOCITY_ARROWS_GRID2 - Fatal error!' );
  end

  xyuv = dtable_data_read ( xyuv_file, dim_two, node_num );

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

  r8mat_transpose_print_some ( dim_two, node_num, xyuv, ...
    1, 1, dim_two, 10, '  Some of the XY and UV data:' );

  dim_num = 2;
%
%  Determine the rectangle that contains the data.
%
  x_min = min ( xyuv(1,1:node_num) );
  x_max = max ( xyuv(1,1:node_num) );
  y_min = min ( xyuv(2,1:node_num) );
  y_max = max ( xyuv(2,1:node_num) );
%
%  Determine the coordinates of the grid points.
%
  if ( number <= 1 )

    x_vec = 0.5 * ( x_min + x_max );
    y_vec = 0.5 * ( y_min + y_max );

  elseif ( y_max - y_min < x_max - x_min ) then

    x_num = number;

    inc = ( x_max - x_min ) / ( x_num - 1 );

    y_num = ( round ( y_max - y_min ) / inc );
    y_frac = 0.5 * ( ( y_max - y_min ) - y_num * inc );

    x_vec = x_min          : inc : x_max;
    y_vec = y_min + y_frac : inc : y_max - y_frac;

  else

    y_num = number;

    inc = ( y_max - y_min ) / ( y_num - 1 );

    x_num = ( round ( x_max - x_min ) / inc );
    x_frac = 0.5 * ( ( x_max - x_min ) - x_num * inc );

    x_vec = x_min + x_frac : inc : x_max - x_frac;
    y_vec = y_min          : inc : y_max;

  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 ( xyuv(1,:), xyuv(2,:), xyuv(3,:), ...
    x_grid, y_grid, 'cubic' );

  v_grid = griddata ( xyuv(1,:), xyuv(2,:), xyuv(4,:), ...
    x_grid, y_grid, 'cubic' );
%
%  Make a vector plot using the grid data.
%
  h = quiver ( x_grid, y_grid, u_grid, v_grid, scale );
%
%  You can change the width of the arrows here.
%  You could also change their color.
%
  set ( h, 'LineWidth', 1 )
%
%  Provide text labels for the plot and the axes.
%
  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 );
%
%  For some reason, if you let MATLAB choose the axes, it's too generous.
%
  axis ( [ x_min, x_max, y_min, y_max ] );
  axis square
%
%  To turn off the axis labeling, tick marks and background,
%  uncomment the command "AXIS OFF"
%
% axis off

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

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