function [BinsX,BinsY,BinsFitY]=binscatterplot(YY,XX,ZZ,nbins,OrderPol,ShallIplot,NamesYX)

if nargin<6
    ShallIplot=0;
end

%% Orthogonalize ZZ
Const = ones(size(YY,1),1);
ZZ = [Const ZZ];

if size(ZZ,2)>1
    betasY = ZZ\YY; 
    OrthY = YY - ZZ(:,2:end)*betasY(2:end);

    betasX = ZZ\XX; 
    OrthX = XX - ZZ(:,2:end)*betasX(2:end);
else
    OrthY = YY;
    OrthX = XX;
end

%% 
% OrthY = OrthY(:);
% OrthX = OrthX(:);

%% Take linear regression YY on XX
OrthX_all = [ones(size(OrthX))];
for pol = 1:OrderPol; 
    OrthX_all = [OrthX_all OrthX.^pol];
end

beta = OrthX_all\OrthY;
FitOrthY = OrthX_all*beta;

%% Calculate Bins
[~,Sorter] = sort(OrthX);
OrthY = OrthY(Sorter,1);
OrthX = OrthX(Sorter,:);
FitOrthY = FitOrthY(Sorter,1);

ObsPerBin = ceil(size(OrthY,1)/nbins);

for jj = 1:nbins
    if ObsPerBin*jj<=size(OrthX,1)
    BinsX(jj,1) = mean(OrthX(ObsPerBin*(jj-1)+1:ObsPerBin*jj,1));
    BinsY(jj,1) = mean(OrthY(ObsPerBin*(jj-1)+1:ObsPerBin*jj));
    BinsFitY(jj,1) = mean(FitOrthY(ObsPerBin*(jj-1)+1:ObsPerBin*jj));
    end
end


% GridX = [min(XX):(max(XX)-min(XX))/(nbins):max(XX)]';
% 
% for jj = 2:size(GridX,1)
%     Select = (XX>=GridX(jj-1) & XX<=GridX(jj));
%     BinsX(jj-1,1) = mean(OrthX(Select,2));
%     BinsY(jj-1,1) = mean(OrthY(Select));
%     BinsFitY(jj-1,1) = mean(FitOrthY(Select));
%     
% end

%% Plot?
if ShallIplot ==1
    figure
    scatter(BinsX,BinsY,'MarkerFaceColor','b','MarkerEdgeColor','b'); 
    hold on; 
    plot(BinsX,BinsFitY,'color','g','linestyle','-','linewidth',2);
    xlabel(NamesYX{1}); 
    ylabel(NamesYX{2}); 
    xlim([min(BinsX)-.01 max(BinsX)+.01]);
    grid on
end
