function [Centroids] = plot_density_scatter(IARwin, quietwin_index, Xdata, xabsoption, xlogoption, Ydata, yabsoption, ylogoption, xdatalabel, ydatalabel, s, datestring, stationstring)
global fntsize xstr titlestr fstr
% Function for plotting density plot, as coloured scatter, contour or
% surface, of two variables, for example the distance a lightning strike is
% away from the ICM station vs the lag time dt that the matching code
% selects for a best match between lightning strike and B peak
timeindexvec = [quietwin_index, 1/IARwin(3)*IARwin(1)+1:1/IARwin(3)*IARwin(2)];
 f1=figure('Position', [0 0 s(3) s(4)]);%('Position', [0 0 s(3) s(4)]); % New figure
 %f2=figure('Position', [0 0 s(3) s(4)]);%('Position', [0 0 s(3) s(4)]); % New figure
for i = 1:length(timeindexvec)
    if strcmp(xabsoption,'abs')
        x = abs(Xdata{timeindexvec(i)});
    else
        x = Xdata{timeindexvec(i)};
    end
    if strcmp(yabsoption,'abs')
        y = abs(Ydata{timeindexvec(i)});
    else
        y = Ydata{timeindexvec(i)};
    end
    if length(x) < 10 || length(y) < 10
        continue
    end
    figure(f1);
    ax(i)=subplot(3,4,i);
    if strcmp(xlogoption,'') && strcmp(ylogoption,'')
        [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y);
        hold on
        dscatter(x,y,'PLOTTYPE','contour');
    elseif strcmp(xlogoption,'log') && strcmp(ylogoption,'')
        [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'logx',true);
        hold on
        dscatter(x,y,'PLOTTYPE','contour','logx',true);
    elseif strcmp(xlogoption,'') && strcmp(ylogoption,'log')
        [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'logy',true);
        hold on
        dscatter(x,y,'PLOTTYPE','contour','logy',true);
    elseif strcmp(xlogoption,'log') && strcmp(ylogoption,'log')
        [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'logx',true,'logy',true);
        hold on
        dscatter(x,y,'PLOTTYPE','contour','logx',true,'logy',true);
    end
    xlabel(xdatalabel);
    ylabel(ydatalabel);
    title(strcat(datestring, stationstring,'\_',num2str(timeindexvec(i)*IARwin(3)-IARwin(3)),'-',num2str(timeindexvec(i)*IARwin(3))));
    
    
   % Median magnetic field B for distance bins
   %medians_B_d(ldist,timeindexvec(i),0,1000,10000)
   
   % Median lEnergy and magnetic field B for LEnergy quartile bins in each
   % distance bin
   
    %Fit = polyfit(x,y,1);
    
%     figure(f2);
%     ax2(i)=subplot(3,4,i);
%     figure(f2);
%     ax2(i)=subplot(3,4,i);
%     if strcmp(xlogoption,'') && strcmp(ylogoption,'')
%         [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'PLOTTYPE','surf');
%     elseif strcmp(xlogoption,'log') && strcmp(ylogoption,'')
%         [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'PLOTTYPE','surf','logx',true);
%     elseif strcmp(xlogoption,'') && strcmp(ylogoption,'log')
%         [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'PLOTTYPE','surf','logy',true);
%     elseif strcmp(xlogoption,'log') && strcmp(ylogoption,'log')
%         [hAxes, F, col, ctrs1, ctrs2] = dscatter(x,y,'PLOTTYPE','surf','logx',true,'logy',true);
%     end
%     xlabel(xdatalabel);
%     ylabel(ydatalabel);
%     title(strcat(datestring, stationstring,'\_',num2str(timeindexvec(i)*IARwin(3)-IARwin(3)),'-',num2str(timeindexvec(i)*IARwin(3))));
%     view([0 0 1]);
%     shading flat;
%     surzmax = zlim;
%     hold on
    
    %Calculate centroids
    % colpeaks = col(col>quantile(col,0.95));
    % TF = islocalmax(colpeaks);
    % ncentroid = size(colpeaks(TF),1);
    %     ncentroid = input('Number of centroids in each subplot including zero:  ');
    %     if ncentroid > 0
    %         opts = statset('Display','final');
    %         % for n=2:11
    %         % IDXf{n-1}=kmeans([x, y],n,'Distance','cityblock','Replicates',5,'Options',opts);
    %         % S(n-1)=mean(silhouette(F, IDXf{n-1}));
    %         % end
    %         [IDX, C]=kmeans([x, y],ncentroid,'Distance','cityblock','Replicates',5,'Options',opts);
    %         plot3(C(:,1),C(:,2),[surzmax(2)*ones(length(C(:,1)),1)],'kx','MarkerSize',15)
    %         Centroids{i} = [C surzmax(2)*ones(length(C(:,1)),1)];
    %         %plot(C(:,1),C(:,2),'kx','MarkerSize',15); (for contour plot)
    %         %plot(x,polyval(Fit,x));
    %         hold off
    %     else
    %         hold off
    %     end
    
    % Calculate Maxima
    nmaxima = input('Number of maxima in each subplot including zero:  ');
    [zmax,imax] = extrema2(F);
    if nmaxima > 0
        [I J] = ind2sub(size(F),imax(1:nmaxima));
        plot3(ctrs1(J),ctrs2(I),zmax(1:nmaxima),'kx','MarkerSize',15)
        Centroids{i} = [ctrs1(J)',ctrs2(I)'];
    else
        Centroids{i} = [NaN, NaN];
        hold off
    end
    view([0 0 zmax(1)]);

%      figure(f2);
%      ax2(i)=subplot(3,4,i);
%      ax2(i)=subplot(3,4,i);
%      quartvec = quantile(y(x<1000),[0.25 0.5 0.75]);
% quartvec(2,:) = quantile(y(x>1000 & y<2000),[0.25 0.5 0.75]);
% quartvec(3,:) = quantile(y(x>2000 & y<3000),[0.25 0.5 0.75]);
% quartvec(4,:) = quantile(y(x>3000 & y<4000),[0.25 0.5 0.75]);
% quartvec(5,:) = quantile(y(x>4000 & y<5000),[0.25 0.5 0.75]);
% quartvec(6,:) = quantile(y(x>5000 & y<6000),[0.25 0.5 0.75]);
% quartvec(7,:) = quantile(y(x>6000 & y<7000),[0.25 0.5 0.75]);
% quartvec(8,:) = quantile(y(x>7000 & y<8000),[0.25 0.5 0.75]);
% quartvec(9,:) = quantile(y(x>8000 & y<9000),[0.25 0.5 0.75]);
% quartvec(10,:) = quantile(y(x>9000 & y<10000),[0.25 0.5 0.75]);
% xl = {'0-1k','1k-2k','2k-3k','3k-4k','4k-5k','5k-6k','6k-7k','7k-8k','8k-9k','9k-10k'};
% plot(quartvec(:,1),'-')
% set(gca,'xticklabel',xl)
% hold on
% plot(quartvec(:,2),'-')
% plot(quartvec(:,3),'-')
end
linkaxes(ax,'xy')
%linkaxes(ax2,'xy')

end
