function varargout = avg_ascii2(files, x, y, varargin)
%  1.  NAME
%
%  avg_ascii2 - Average 3-Dimensional data from different input files
%
%  2.  SYNOPSIS
%
%  [new, grid] = avg_ascii2(files, x, y, grid, method )
%
%  3.  DESCRIPTION
%
%  This Matlab function calculates averages of data between different
%  ascii columnar files. Each column specifies the values of a variable
%  at different grid points.
%
%  Parameters x and y specify the columns to be used as x/y points. A
%  common grid to all datasets is determined, and data are regridded onto
%  this new grid. If grid is specified, it is taken as the new grid.
%
%  The interpolation method for regridding can be specified by method.
%
%  Since regridding is done by regrid2(3m), bad data along the regrid
%  axes is allowed, and duplicate grid points on input will be averaged
%  together.
%
%  4.  PARAMETERS
%
%     files
%        Cell array of file names containing columnar ascii data to be
%        averaged together. Files may also be a cell array of already
%        processed data.
%
%     x  Column number to use for regridding data along a common x grid
%        before averaging.
%
%     y  Column number to use for regridding data along a common y grid
%        before averaging.
%
%     grid optional
%        Grid on which data should be averaged (default is to build a
%        default grid).  It must be a cell array with 2 column vectors
%        specifying the new grid points where x=grid{1} and y=grid{2}.
%
%     method optional
%        Method to use for regridding data (default is linear). See
%        Matlab function griddata for more information.
%
%  5.  RETURNS
%
%     new
%        Array with averaged data on the new grid.
%
%     grid optional
%        The new grid points.
%
%  6.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

%--- Default values:
method = 'linear';
grid   = [];

%--- Arguments:
if (nargin < 3) | (nargin > 5),
	error([ 'Wrong number of arguments' ]);
end
if nargin > 3,
	grid = varargin{1};
end
if nargin > 4,
	method = varargin{2};
end

%--- Initialization:
n_files = length(files);
if n_files < 2,
	error([ 'Too few files/data for averaging data' ]);
end

if isstr(files{1}),
	%--- Read the data:
	for f=1:n_files,
		data{f} = tonan(read_asc(files{f}));	
	end
else
	data = files;
end

%--- Check if all files have the same number of columns:
nc = size(data{1}, 2);
for f=2:n_files,
	if nc ~= size(data{f},2),
		error([ 'Files ' int2str(1) ' and ' int2str(f) ...
					' have different column number' ]);
	end
end
vv = find( ((1:nc) ~= x) & ((1:nc) ~= y) );

%--- Check that regrid columns are within the columns of data:
if x > nc | y > nc,
	error([ 'Regrid column outside range' ]);
end

%--- Get the grid values from input files and removes undefined positions:
for f=1:n_files,
	[ grd{f}, data{f} ] = ascii_grd(data{f}, [x,y]);
end

%--- Get the most included ranges:
if isempty(grid),
	[xmin, xmax, n] = com_grd(data, x);
	grid{1} = new_grid(xmin, xmax, n);
	[xmin, xmax, n] = com_grd(data, y);
	grid{2} = new_grid(xmin, xmax, n);
end
nx = length(grid{1});
ny = length(grid{2});

%--- Regrid data and reshape to have it in form of a contour data:
new = NaN * ones(ny*nx, nc, n_files);
for f=1:n_files,
	tmp = regrid2(data{f}, x, y, grid, method);
	new(:,:,f) = reshape(tmp, nx*ny, nc);
end
clear data;
data = new;

%--- Average data:
nl = size(data,1);
new = NaN * ones(nl, nc);
for n=1:nl,
	for c=1:nc,
		igood = find(~isnan(data(n,c,:)));
		if ~isempty(igood),
			new(n,c) = mean(data(n,c,igood));
		end
	end
end

%--- Remove lines with only bad data:
ibad = ones(nl, 1);
for v=vv,
	ibad = ibad & isnan(new(:,v));
end
igood = find(~ibad);
new = new(igood,:);

%--- Return values:
if nargout > 0,
	varargout{1} = new;
end
if nargout > 1,
	varargout{2} = grid;
end

return;

	






