function PlotPSDMean(PSD, D_rep, colormap)
% PlotPSDMean(PSD, D_rep)
% Given a series of PSD, and their size
% Plot the mean PSd of them
N_dvdlnd = length(PSD.ustar);

for i = 1:size(PSD.dvdlnd,2)
    temp_dvdlnd = PSD.dvdlnd(:,i);
    temp_dvdlnd_err = PSD.dvdlnd_err(:,i);
    
    dvdlnd_p(i) = mean(temp_dvdlnd);
    dvdlnd_err_p(i) = sqrt(sum(temp_dvdlnd_err.^2))/N_dvdlnd; % error propogation
end

plot(D_rep, dvdlnd_p,'o','MarkerFaceColor','none','MarkerEdgeColor',colormap,'Markersize',7,'linewidth',1); hold on;
plot(D_rep, dvdlnd_p,'-','color',colormap,'linewidth',1); hold on;

for i = 1:size(PSD.dvdlnd,2)
    tempx = [D_rep(i), D_rep(i)];
    tempy = [dvdlnd_p(i)-dvdlnd_err_p(i),dvdlnd_p(i)+dvdlnd_err_p(i)];
    plot(tempx, tempy, '-','color',colormap); hold on;
end
%--- add Kok (2011) theory ---
Ds = linspace(0.1,20,1000);
Ds_avge = 3.4;
Sigma = 3.0;
lam = 12;
cv = 12.62;
JK_dvdlnd = Ds/cv.*(1+erf(log(Ds/Ds_avge)/sqrt(2)/log(Sigma))).*exp(-(Ds/lam).^3);
JK2011 = plot(Ds,JK_dvdlnd,'b-.'); hold on;
%--- plot setting---
set(gca,'xscale','log');
set(gca,'yscale','log');

axis([0.4,19,5*10^(-4),1.9]);
set(gca,'xtick',[0.4,1,10,19]);
set(gca,'xticklabel',{'0.4','1','10',''});
set(gca,'ytick',[10^(-3),10^(-2),10^(-1),10^(0)]);
set(gca,'yticklabel',{'10^{-3}','10^{-2}','10^{-1}','10^{0}'});

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