function [D,X,Y] = binning_diff(xin,yin,zin,xin2,yin2,kin,xres,yres,westl,eastl,southl,northl)

% Bins data into a grid of equal resolution gres in x and y
% xin, yin and zin are column vector all with the same size
% global northl southl westl eastl xres yres
% Aug 2, 2011. Adapted for the entire Southern Ocean. Gets rid of profiles
% or data points outside the grid (such that we have nonzero values in
% xbin, ybin and hence accumarray doesn't give error)
% March 11, 2019: Modified two compute the difference between 2 sample means
% for purpo))

% % To test out:
% % Input grid paraeters:
% xres=2; % in degrees
% yres=2;
% westl=-180; eastl=180;
% northl=90; southl=-90;

% % Input data to be gridded (column vectors of same size)
% obsMLD=load('/glade/work/magdac/argo_copernicus/ARGO_MLDbmax_TEOS10_Copernicus_PF_2000-2017_all.mat');

% xin=obsMLD.LONr(:);
% yin=obsMLD.LATr(:);
% zin=obsMLD.MLD_bmax(:);
% % kin=obsMLD.MLD_bmax(:);
% % if using the same dataset for zin and kin, difference should be zero
% % test subsamplig zin, to check it will work for zin and kin of different dimensions
% kin=zin(1:5:end);
% xin2=xin(1:5:end); % then, dimension should be consistent.. need to input xin and yin subsampled to use with kin
% yin2=yin(1:5:end);

%%  Define Grid
xedges=westl:xres:eastl;
yedges=southl:yres:northl;
lnx=length(xedges); lny=length(yedges);
%xborder=westl*ones(lny,1);
%yborder=southl*ones(lnx,1);
gsX=abs(westl-eastl)/xres; % Grid Size will be 151 x 76
gsY=abs(northl-southl)/yres;


% Get rid of all profiles outside the grid boundaries
        % Define polygon
        RLat=[northl northl southl southl northl];
        RLon=[westl eastl eastl westl westl]; % W long

        % Extract region
        in=inpolygon(xin,yin,RLon,RLat);
        
         xin=xin(in);
         yin=yin(in);
         zin=zin(in);
         
         in2=inpolygon(xin2,yin2,RLon,RLat);
         xin2=xin2(in2);
         yin2=yin2(in2);
         kin=kin(in2);


% Create a border with nan such that in every bin there's a value (or nan),
% otherwise accumarray doesn't work (i.e. a dummy border that then we crop)
y = [yin; yedges']; % to add border with NaN values in z
x = [xin; xedges'];
    
if length(y)==length(x)
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(y)>length(x) % adjust x to the length of y
    nrep=ceil((length(y)-length(x))./length(xedges));
    xfaked=repmat(xedges,1,nrep);
    x = [x; xfaked']; 
    x = x(1:length(y));
    z = [zin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(x)>length(y) % adjust y to the length of x
    nrep=ceil((length(x)-length(y))./length(yedges));
    yfaked=repmat(yedges,1,nrep);
    y = [y; yfaked'];
    y = y(1:length(x));
    z = [zin; nan(lnx,1)]; % need to modify the faked brder accordingly if grid size changes!
end
% Note: 
% % x,y,z and k have to be all the same size
% % need to to add a border of nan values to the grid so that at list there's
% % 1 value (or nan) in each bin

% Use histc to get the index for each bin.
[nx, xbin]=histc(x,xedges);
[ny, ybin]=histc(y,yedges);

% Repeat for K values
y2 = [yin2; yedges']; % to add border with NaN values in z
x2 = [xin2; xedges'];
    
if length(y2)==length(x2)
    k = [kin; nan(length(y2),1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(y2)>length(x2) % adjust x to the length of y
    nrep=ceil((length(y2)-length(x2))./length(xedges));
    xfaked=repmat(xedges,1,nrep);
    x2 = [x2; xfaked']; 
    x2 = x2(1:length(y2));
    k = [kin; nan(lny,1)]; % need to modify the faked brder accordingly if grid size changes!
elseif length(x2)>length(y2) % adjust y to the length of x
    nrep=ceil((length(x2)-length(y2))./length(yedges));
    yfaked=repmat(yedges,1,nrep);
    y2 = [y2; yfaked'];
    y2 = y2(1:length(x2));
    k = [kin; nan(lnx,1)]; % need to modify the faked brder accordingly if grid size changes!
end
% Note: 
% % x,y,z and k have to be all the same size
% % need to to add a border of nan values to the grid so that at list there's
% % 1 value (or nan) in each bin

% Use histc to get the index for each bin.
[nx, xbin2]=histc(x2,xedges);
[ny, ybin2]=histc(y2,yedges);


% accumarray computes nanmean on z values corresponding to each bin.
A = accumarray([xbin ybin],z,[],@nanmean);
A=A';
A(A==0)=NaN;
% Remove faked edge
A=A(1:end-1,1:end-1);

% accumarray computes nanmean on k values corresponding to each bin.
B = accumarray([xbin2 ybin2],k,[],@nanmean);
B=B';
B(B==0)=NaN;
% Remove faked edge
B=B(1:end-1,1:end-1);

% Save the difference between z and k values in each bin
D=A-B; % returns difference between two sample means in each bin

% Could implement the use of hist3, instead of histc twice
% edges{1} = xedges;
% edges{2} = yedges;
% N = hist3([x y],'Edges',edges); %hist3 ignores NaN
% N=N';
% meanZ = A./N;


% Jan 27, 2016: Added grid as output:
% Grid to plot in a map
[X,Y]=meshgrid(westl+xres/2:xres:eastl-xres/2,southl+yres/2:yres:northl-yres/2);
X=X-xres/2;
Y=Y-yres/2;