function plot_fig = grid_tension ( datafile, ....
                                   variable, ...
				   data_grid, ...
				   xgrid, ...
				   ygrid, ... 
				   anisotropy, ...
				   tolerance, ...
			           omega, ...
				   nrange )
% GRID_TENSION:  Plots results of gridding at different tensions.
% 
% Usage:  grid_tension ( datafile,
%                        variable,
%                        data_grid,
%                        xgrid,
%                        ygrid,
%                        [anisotropy] )
% 
% Input:
%          datafile:  Name of file which contains the variables.
%          variable:  Name of the variable we want to have gridded.
%          data_grid:  What is to be gridded.
%          xgrid:  x values of data_grid.
%          ygrid:  y values of data_grid.
%          anisotropy:  dx/dy.  Optional scalar, default = 1.0.
%          tolerance:  when to stop gridding.  Optional, default = 5.0.
%          omega:  SOR parameter.  Optional, default = 1.4.
%          nrange:  spread parameter.  Optional, default = 10.
%
% Output:
%          handle of figure window containing plots
%         
% If the anisotropy is not specified when this function is invoked, it is
% set to 1.
%

if ( nargin < 5 )
  help rgscript
  return;
  end
if ( nargin < 6 )
  anisotropy = 1.0;
end
if ( nargin < 7 ) 
  tolerance = 5.0;
end
if ( nargin < 8 )
  omega = 1.4;
end
if ( nargin < 9 )
  nrange = 10;
end

% load in the datafile which should contain the variable in question
eval ( sprintf ( 'load %s', datafile ) );


plot_name = sprintf ( 'grid_tension:  anisotropy = %.3f', anisotropy );
plot_fig = figure ('name', plot_name, ...
                   'position', [100 100 800 800] );

for i = 1:4
  subplot(2,2,i)
  tension = (i-1)/3; 
  new_grid = ngrid ( data_grid, ...
                     tension, ...
		     anisotropy, ...
		     tolerance, ...
		     omega, ...
		     nrange );
  grid_handle = surf(xgrid,ygrid,new_grid);
  shading interp
  set(gca,'view',[0 90]);
  hold on 
  colormap(cool)
  set(gca,'clim',[0 100]);
  set(gca, 'aspectratio', [1 1] );
  
  % Need to make sure that the original data points are on top of surface.
  % 1e10 z values ought to do it.
  plot3 ( x, y, 1e10*ones(size(x)), '+k' );
  title( sprintf('tension = %f', tension) );
end

return

