function PlotPSD_in3Groups(Struct, D_dust)
% Input the PSD struct, OPC bin boundary sizes 
% Output the averaged PSD in the N groups

% set plot scale
axisset = {[0.5, 10],[4*10^(-2),2]};
% first order the group with increaing shear velocity order
[ustar, index_ustar] = sort(Struct.ustar,'ascend');

for i = 1:size(Struct.dvdlnd,2) % loop over size bins
    dvdlnd(:,i) = Struct.dvdlnd(index_ustar,i);
    dvdlnd_err(:,i) = Struct.dvdlnd_err(index_ustar,i);
end
    
% decide size of the N groups
N_ustar = length(ustar);

N = 3;
N1 = floor(N_ustar/N);
N2 = floor((N_ustar-N1)/(N-1));
N3 = N - N1 - N2;

N1_index = 1:1:N1;
N2_index = (N1+1):1:(N1+N2);
N3_index = (N1+N2+1):1:N_ustar;

% distribute data into the groups
Group1.dvdlnd = dvdlnd(N1_index,:);
Group1.dvdlnd_err = dvdlnd_err(N1_index,:);
Group1.ustar = ustar(N1_index);

Group2.dvdlnd = dvdlnd(N2_index,:);
Group2.dvdlnd_err = dvdlnd_err(N2_index,:);
Group2.ustar = ustar(N2_index);

Group3.dvdlnd = dvdlnd(N3_index,:);
Group3.dvdlnd_err = dvdlnd_err(N3_index,:);
Group3.ustar = ustar(N3_index);

% mean PSD within each of the 3 groups
for i = 1:size(Struct.dvdlnd,2) % loop over size bins
    Group1.dvdlnd_mean(i) = mean(Group1.dvdlnd(:,i));
    Group2.dvdlnd_mean(i) = mean(Group2.dvdlnd(:,i));
    Group3.dvdlnd_mean(i) = mean(Group3.dvdlnd(:,i));
    
    Group1.dvdlnd_err_mean(i) = sqrt(sum(Group1.dvdlnd_err(:,i).^2))/length(N1_index);
    Group2.dvdlnd_err_mean(i) = sqrt(sum(Group2.dvdlnd_err(:,i).^2))/length(N2_index);
    Group3.dvdlnd_err_mean(i) = sqrt(sum(Group3.dvdlnd_err(:,i).^2))/length(N3_index);    
end
Group1.ustar_mean = mean(Group1.ustar);
Group2.ustar_mean = mean(Group2.ustar);
Group3.ustar_mean = mean(Group3.ustar);

% plot the three PSDs
g1_p = plot(D_dust, Group1.dvdlnd_mean,'>','MarkerFaceColor','none','MarkerEdgeColor',[1,0,0.4],'Markersize',7,'linewidth',1); hold on;
plot(D_dust, Group1.dvdlnd_mean,'-','Color',[1,0,0.4],'linewidth',1); hold on;

g2_p = plot(D_dust, Group2.dvdlnd_mean,'s','MarkerFaceColor','none','MarkerEdgeColor',[0,0.6,0],'Markersize',7,'linewidth',1); hold on;
plot(D_dust, Group2.dvdlnd_mean,'-','Color',[0,0.6,0],'Markersize',7,'linewidth',1); hold on;

g3_p = plot(D_dust, Group3.dvdlnd_mean,'<','MarkerFaceColor','none','MarkerEdgeColor',[0.4,0.6,1],'Markersize',7,'linewidth',1); hold on;
plot(D_dust, Group3.dvdlnd_mean,'-','Color',[0.4,0.6,1],'Markersize',7,'linewidth',1); hold on;

%plot errorbar
for i = 1:size(Struct.dvdlnd,2) % loop over size bins
    tempx = [D_dust(i), D_dust(i)];
    tempy1 = [Group1.dvdlnd_mean(i)-Group1.dvdlnd_err_mean(i), Group1.dvdlnd_mean(i)+Group1.dvdlnd_err_mean(i)];
    tempy2 = [Group2.dvdlnd_mean(i)-Group2.dvdlnd_err_mean(i), Group2.dvdlnd_mean(i)+Group2.dvdlnd_err_mean(i)];
    tempy3 = [max(Group3.dvdlnd_mean(i)-Group3.dvdlnd_err_mean(i),axisset{2}(1)), Group3.dvdlnd_mean(i)+Group3.dvdlnd_err_mean(i)]; % to prevent a negative value under log10 scale
    plot(tempx, tempy1,'-','color',[1,0,0.4],'linewidth',1); hold on;
    plot(tempx, tempy2,'-','color',[0,0.6,0],'linewidth',1); hold on;
    plot(tempx, tempy3,'-','color',[0.4,0.6,1],'linewidth',1); hold on;
end

%set the legend
thelegend = legend([g1_p g2_p g3_p],['Oceano Dunes (u_{*} = ',num2str(Group1.ustar_mean,'%.2f'),' m/s)'],...
        ['Oceano Dunes (u_{*} = ',num2str(Group2.ustar_mean,'%.2f'),' m/s)'],...
        ['Oceano Dunes (u_{*} = ',num2str(Group3.ustar_mean,'%.2f'),' m/s)']);    
set(thelegend,'location','northwest','FontSize',7,'FontWeight','bold');

%plot setting
set(gca,'xscale','log');
set(gca,'yscale','log'); 
xlim([axisset{1}(1),axisset{1}(2)]);
ylim([axisset{2}(1),axisset{2}(2)]);

set(gca,'xtick',[0.5,1,10]);
set(gca,'xticklabel',{'0.5','1','10'});
set(gca,'ytick',[4e-2,10^(-1),10^(0),2]);
set(gca,'yticklabel',{'4E-2','10^{-1}','1','2'});

xlabel('Dust geometric diameter \itD_{g}, \rm\bf\mum','FontSize',14,'FontWeight','bold');
ylabel('Norm. volume size distr. \itdV/dlnD_{g}','FontSize',14,'FontWeight','bold');

end