function vector_magnitude_grid ( node_file, vector_file, grid_num, contour_num )

%% VECTOR_MAGNITUDE_GRID contour plots a vector magnitude from 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:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%    Hyung-Chun Lee
%
%  Usage:
%
%    vector_magnitude_grid ( node_file, vector_file, contour_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.
%
%    Input, integer CONTOUR_NUM, the number of contour lines to use.
%
  fprintf ( 1, '\n' );
  timestamp;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'VECTOR_MAGNITUDE_GRID:\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Display a contour plot of vector magnitude,\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 contour lines.
%
  if ( nargin < 4 )
    fprintf ( 1, '\n' );
    contour_num = input ( '  Enter the number of contour lines.' );
  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_MAGNITUDE_GRID - Fatal error!\n' );
    fprintf ( 1, '  Node data must have spatial dimension 2.\n' );
    error ( 'VECTOR_MAGNITUDE_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_MAGNITUDE_GRID - Fatal error!\n' );
    fprintf ( 1, '  The vector file must list pairs of numbers.\n' );
    error ( 'VECTOR_MAGNITUDE_GRID - Fatal error!' );
  end

  if ( node_num2 ~= node_num )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'VECTOR_MAGNITUDE_GRID - Fatal error!\n' );
    fprintf ( 1, '  The number of vectors and nodes do not match.\n' );
    error ( 'VECTOR_MAGNITUDE_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' );
%
%  Compute the vector magnitude.
%
  uv_norm_grid(:,:) = sqrt ( u_grid(:,:).^2 + v_grid(:,:).^2 );
%
%  Make a vector magnitude contour plot using the grid data.
%
  contour ( x_grid, y_grid, uv_norm_grid, contour_num )

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

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

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

  axis square

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

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