function [dvdlnd,dvdlnd_err] = PlotPSDeaOPC(Struct, Bin_D, D_dust_range, rou, i)
% PlotPSDeaOPC(Struct)
% Input number concentration (#/m^3), representiative geometric diameter
% Bin_D (unit: micron), dust density rou (kg/m^3), plot the norm. volume
% PSD at each of the i^{th} heights

z = Struct.z(:,i);
z_legend = mean(z);
ustar = Struct.SaltationSand.ustRe; 
cm = [];
cm_err = [];

% define colormap
colormap = {[1,0,0.2],[0.2,0.6,0.2],[0,0.6,0.8],[0.8,0.6,0.4]};
% red, green, blue, brown

% put 7 bins of the i^{th} together
for j = 1:7
    cn_temp = Struct.(['bin',num2str(j)]).y(:,i); % in unit of #/m^3
    cn_err_temp = Struct.(['bin',num2str(j)]).y_err(:,i);% in unit of #/m^3 
        
    cm_temp = NumberConcen2MassConcen(cn_temp, Bin_D(j), rou); % Number concentration (#/m^3) to mass concentration on kg/m^3
    cm_err_temp = NumberConcen2MassConcen(cn_err_temp, Bin_D(j), rou); % Number concentration (#/m^3) to mass concentration on kg/m^3
        
    cm = [cm, cm_temp]; 
    cm_err = [cm_err, cm_err_temp]; 
end

% calculate sum of the 7 bins, in unit of kg/m^3
for j = 1:length(ustar)
    cm(j,8) = sum(cm(j,1:7));
    cm_err(j,8) = sum(cm_err(j,1:7).^2);
end

% normalize mass flux
for j = 1:7
    cm_norm(:,j) = cm(:,j)./cm(:,8); 
    cm_err_norm(:,j) = sqrt((cm_err(:,j)./cm(:,8)).^2 + (-1*cm(:,j).*cm_err(:,8)./cm(:,8)./cm(:,8)).^2); % error propagation   
end

% prepare dlnd
dlnd = log(D_dust_range(2:end))-log(D_dust_range(1:(end-1)));

% calculate the mean PSD
for j = 1:7
    dvdlnd(:,j) = cm_norm(:,j)/dlnd(j);
    dvdlnd_err(:,j) = cm_err_norm(:,j)/dlnd(j);
end

% calculate the mean PSD
for j = 1:7
    dvdlnd_mean(j) = mean(dvdlnd(:,j));
    dvdlnd_err_mean(j) = sqrt(sum(dvdlnd_err(:,j).^2))/length(ustar);
end

% plot PSD
plot(Bin_D, dvdlnd_mean,'o','MarkerFaceColor',colormap{i},'MarkerEdgeColor','none','Markersize',7,'linewidth',1); hold on;
plot(Bin_D, dvdlnd_mean,'-','color',colormap{i},'linewidth',1); hold on;

for j = 1:7
    tempx = [Bin_D(j), Bin_D(j)];
    tempy1 = [dvdlnd_mean(j)-dvdlnd_err_mean(j),dvdlnd_mean(j)+dvdlnd_err_mean(j)];
    plot(tempx,tempy1,'-','color',colormap{i},'linewidth',1); hold on;
end

% plot setting
set(gca,'xscale','log');
set(gca,'yscale','log');
xlim([0.4 19]);
ylim([0.05 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