function output = interpBedmachineAntarctica(X,Y,string,method,ncdate),

if nargin<3, string = 'bed'; end
if nargin<4
	if strcmp(string,'mask') | strcmp(string,'source')
      method='nearest'; % default method
   else
      method='cubic'; % default method
   end
end
if nargin<5
	%ncdate='2014-03-12';
	%ncdate='2014-04-28';
	%ncdate='2014-07-31';
	%ncdate='2014-09-23';
	%ncdate='2015-11-06';
	%ncdate='2015-12-09';
	ncdate='2016-06-10';
	ncdate='2016-08-26';
	ncdate='2016-11-17';
	ncdate='2017-01-05';
	ncdate='2017-03-10';
	ncdate='2018-02-02';
	ncdate='2018-03-20';
	ncdate='2018-05-18';
	ncdate='2018-06-08';
	ncdate='2018-09-26';
	ncdate='2018-11-14';
	ncdate='2019-01-30';
	ncdate='2019-03-12';
	ncdate='2019-04-02';
	ncdate='2019-04-15'; % no need of firn correction
	ncdate='2020-06-18';
	ncdate='2020-07-15'; %BedMachine v2?
end

basename = 'BedMachineAntarctica';

%read data
switch (oshostname()),
	case {'ronne'}
		nc=['/home/ModelData/Antarctica/BedMachine/' basename '-' ncdate '.nc'];
	case {'thwaites','murdo','astrid'}
		nc=['/u/astrid-r1b/ModelData/BedMachine/' basename '-' ncdate '.nc'];
	otherwise
		error('hostname not supported yet');
end

if nargout==2,
	string = 'bed';
end

disp(['   -- BedMachine Antarctica version: ' ncdate]);
xdata = double(ncread(nc,'x'));
ydata = double(ncread(nc,'y'));

offset=2;

xmin=min(X(:)); xmax=max(X(:));
posx=find(xdata<=xmax);
id1x=max(1,find(xdata>=xmin,1)-offset);
id2x=min(numel(xdata),posx(end)+offset);

ymin=min(Y(:)); ymax=max(Y(:));
posy=find(ydata>=ymin);
id1y=max(1,find(ydata<=ymax,1)-offset);
id2y=min(numel(ydata),posy(end)+offset);

if strcmp(string,'icemask'),
	disp(['   -- BedMachine Antarctica: loading ' string]);
	%data  = double(ncread(nc,'mask'))';
	data  = double(ncread(nc,'mask',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';
	xdata=xdata(id1x:id2x);
	ydata=ydata(id1y:id2y);
	%ice ocean interface is between 0 and 3, so we might get some 1 by interpolating
	data(find(data==3))=0;
else
	disp(['   -- BedMachine Antarctica: loading ' string]);
	%data  = double(ncread(nc,string))';
	data  = double(ncread(nc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';
	xdata=xdata(id1x:id2x);
	ydata=ydata(id1y:id2y);
end

disp(['   -- BedMachine Antarctica: interpolating ' string]);
disp(['       -- Interpolation method: ' method]);
if strcmp(string,'mask') | strcmp(string,'source'),
	%Need nearest neighbor to avoid interpolation between 0 and 2
	tic
	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest');
	toc
	%tic
	%output = FastInterp(xdata,ydata,data,X,Y,'nearest');
	%toc
else
	%disp('InterpFromGrid');
	%tic
	%output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'cubic'); 
	output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),method); % now the interpolation method can be defined by the user
	%toc
	%disp('FastInterp');
	%tic
	%output = FastInterp(xdata,ydata,data,X,Y,'bilinear');
	%toc
end

end
function zi = FastInterp(x,y,data,xi,yi,method)

	%get data size
	[M N] = size(data);

	% Get X and Y library array spacing
	ndx = 1/(x(2)-x(1));    ndy = 1/(y(2)-y(1));
	% Begin mapping xi and yi vectors onto index space by subtracting library
	% array minima and scaling to index spacing

	xi = (xi - x(1))*ndx;       yi = (yi - y(1))*ndy;

	% Fill Zi with NaNs
	zi = NaN(size(xi));

	if strcmpi(method,'nearest'),
		% Find the nearest point in index space
		rxi = round(xi)+1;  ryi = round(yi)+1;
		% Find points that are in X,Y range
		flag = rxi>0 & rxi<=N & ~isnan(rxi) & ryi>0 & ryi<=M & ~isnan(ryi);
		% Map subscripts to indices
		ind = ryi + M*(rxi-1);
		zi(flag) = data(ind(flag));

	else %Bilinear

		% Transform to unit square
		fxi = floor(xi)+1;  fyi = floor(yi)+1; % x_i and y_i
		dfxi = xi-fxi+1;    dfyi = yi-fyi+1;   % Location in unit square

		% flagIn determines whether the requested location is inside of the data arrays
		flagIn = fxi>0 & fxi<N & ~isnan(fxi) & fyi>0 & fyi<M & ~isnan(fyi);

		%Toss all out-of-bounds variables now to save time
		fxi  = fxi(flagIn);  fyi  = fyi(flagIn);
		dfxi = dfxi(flagIn); dfyi = dfyi(flagIn);

		%Find bounding vertices
		ind1 = fyi + M*(fxi-1);     % indices of (  x_i  ,  y_i  )
		ind2 = fyi + M*fxi;         % indices of ( x_i+1 ,  y_i  )
		ind3 = fyi + 1 + M*fxi;     % indices of ( x_i+1 , y_i+1 )
		ind4 = fyi + 1 + M*(fxi-1); % indices of (  x_i  , y_i+1 )

		% Bilinear interpolation
		zi(flagIn) = ...
			data(ind1).*(1-dfxi).*(1-dfyi) + ...
			data(ind2).*dfxi.*(1-dfyi) + ...
			data(ind4).*(1-dfxi).*dfyi + ...
			data(ind3).*dfxi.*dfyi;
	end
end
