function [] = modifyELCOMbathyGridSize(bathFileIn,bathFileOut,dxNew,dyNew,method)

% function [] = modifyELCOMbathyGridSize(bathFileIn,bathFileOut,dx,dy,method)
%
% Inputs:
%		bathFileIn  : filename of bathymetry file to convert
%		bathFileOut : output file name
%   dx          : new dx values (can be single number or array for plaid grid)
%   dy          : new dy values (can be single number or array for plaid grid)
%   method(optional) : method used in griddata for interpolating data
%
% Uses:
%		readELCOMbathy.m
%		writeELCOMbathy.m
%
% Written by C. Dallimore 16 June 06
%

if nargin == 4
	method = 'linear';
elseif nargin ==5
else
	error('Incorrect number of arguments to modifyELCOMbathyGridSize')
end

% Get Bathymetry from file
bathStruct = readELCOMbathy(bathFileIn);

bath = (bathStruct.bathData);

% Set land and open bc cells to the maximum height plus a small number
pts = find(bath == bathStruct.land_value | ...
					 bath == bathStruct.open_value);

% First set land cells to NaN in case Land value is greater than the maximum height
bath(pts) = NaN;
minHeight = min(min(bath));
maxHeight = max(max(bath));
maxHeight+(maxHeight-minHeight)/100
bath(pts) = maxHeight+(maxHeight-minHeight)/100;


% set up x and y array for old file
x_old = zeros(bathStruct.x_rows,1);
y_old = zeros(bathStruct.y_columns,1);

if ~isempty(bathStruct.dx)
	x_old(1) = bathStruct.dx(1)/2
	for i =2:bathStruct.x_rows
		x_old(i) = x_old(i-1)+(bathStruct.dx(i)+bathStruct.dx(i-1))/2;
	end
else
	x_old(:) = bathStruct.x_grid/2+(0:bathStruct.x_rows-1)*bathStruct.x_grid;
end
if ~isempty(bathStruct.dy)
	y_old(1) = bathStruct.dy(1)/2
	for i =2:bathStruct.y_columns
		y_old(i) = y_old(i-1)+(bathStruct.dy(i)+bathStruct.dy(i-1))/2;
	end
else
	y_old(:) = bathStruct.y_grid/2+(0:bathStruct.y_columns-1)*bathStruct.y_grid;
end

% set up x and y array for new file
maxX = max(x_old);

if length(dxNew) ==1
	i = 1;
	x_new(i) = dxNew/2;
	while x_new(i) < maxX
		i = i+1;
		x_new(i) = x_new(i-1)+dxNew;
	end
else
	x_new(1) = dyNew(1)/2
	for i =2:length(dxNew)
		x_new(i) = x_new(i-1)+(dyNew(i)+dyNew(i-1))/2;
	end
end

if length(dyNew) ==1
	maxY = max(y_old);
	i = 1;
	y_new(i) = dyNew/2;
	while y_new(i) < maxY
		i = i+1;
		y_new(i) = y_new(i-1)+dyNew;
	end
else
	y_new(1) = dyNew(1)/2
	for i =2:length(dyNew)
		y_new(i) = y_new(i-1)+(dyNew(i)+dyNew(i-1))/2;
	end
end

% Interpolate the data using meshgrid
[X,Y]=meshgrid(x_new,y_new);
bath_new=griddata(x_old,y_old,bath',X,Y,method);

% Set point above the old max to land value
pts = find(bath_new > maxHeight | isnan(bath_new));
bath_new(pts) = bathStruct.land_value;

% Setup the new structure
bathStructNew = bathStruct;
bathStructNew.bathData=bath_new';

bathStructNew.x_grid=dxNew;
bathStructNew.y_grid=dyNew;

bathStructNew.x_rows = length(x_new);
bathStructNew.y_columns = length(y_new);

% Write the output file
writeELCOMbathy(bathFileOut,bathStructNew)