function PlotMassConcenProfile(Struct, Bin_D, rou, N, path)
% PlotMassConcenProfile(Struct, Bin_D, rou, N, path)
% Input the number concentration (#/m^3) of N^{th}, representiative geometric 
% diameter Bin_D (unit: micron), dust density rou (kg/m^3), plot the 
% mass concentration (kg/m^3)profile, and print into the folder of the path

z_temp = Struct.z(:,1:4);
cn_temp = Struct.(['bin',num2str(N)]).y(:,1:4);
cn_err_temp = Struct.(['bin',num2str(N)]).y_err(:,1:4);

cm_temp = NumberConcen2MassConcen(cn_temp, Bin_D(N), rou); % Number concentration (#/m^3) to mass concentration on kg/m^3
cm_err_temp = NumberConcen2MassConcen(cn_err_temp, Bin_D(N), rou); % Number concentration (#/m^3) to mass concentration on kg/m^3

ustar_temp = Struct.SaltationSand.ustRe;
Q_temp = Struct.SaltationSand.Q;

% in the order of increasing of shear velocity
[ustar, index_ustar] = sort(ustar_temp, 'ascend');
z = z_temp(index_ustar, :);
c = cm_temp(index_ustar, :); % mass concentration (kg/m^3)
c_err = cm_err_temp(index_ustar, :); % mass concentration (kg/m^3)
Q = Q_temp(index_ustar);

% in preparation of power fitting
for i = 1:4
    z_z0_power(:,i) = log10(z(:,i)./z(:,1)); % normalization and log scale in preparation for plotting 
end
c_power = log10(c);
c_err_power = sqrt((c_err./c/log(10)).^2); % error propagation

% power fitting
for i = 1:length(ustar)
   [a, b, ~, ~, ~] =  Linearfit(z_z0_power(i,:), c_power(i,:), c_err_power(i,:));
   c_fit_power(i,:) = a + b*z_z0_power(i,:); 
   c_fit(i,:) = 10.^c_fit_power(i,:);
end

% plot figures
N_per_fig = 6; % plot 6 lines per figures
colormap = {[1,0,0.2],[1,0.4,0.2],[0.2,0.6,0.2],[0,0.6,0.8],[0.4,0.4,0.8],[0.8,0.6,0.4]};
% red, orange, green, blue, purple, brown
N_fig = floor(length(ustar)/N_per_fig); % number of plot with 6 lines
N_line = length(ustar) - N_fig*N_per_fig; % number of lines in the last plot

for i = 1:N_fig
    figure(i)
    index_N = (i-1)*6;
    for j = 1:6
        plot(c(index_N+j,:), z(index_N+j,:),'o','MarkerFaceColor',colormap{j},'MarkerEdgeColor','none','Markersize',7); hold on;
        plot(c_fit(index_N+j,:), z(index_N+j,:),'-','Color',colormap{j},'linewidth',1); hold on;
        for k = 1:4
            tempx = [c(index_N+j,k)-c_err(index_N+j,k), c(index_N+j,k)+c_err(index_N+j,k)];
            tempy = [z(index_N+j,k), z(index_N+j,k)];
            plot(tempx, tempy, '-','color',colormap{j},'linewidth',1); hold on;
        end
    end
    
    % plot legend
    p_1 = plot(c(index_N+1,1), z(index_N+1,1),'o','MarkerFaceColor',colormap{1},'MarkerEdgeColor','none','Markersize',7); hold on;
    p_2 = plot(c(index_N+2,1), z(index_N+2,1),'o','MarkerFaceColor',colormap{2},'MarkerEdgeColor','none','Markersize',7); hold on;
    p_3 = plot(c(index_N+3,1), z(index_N+3,1),'o','MarkerFaceColor',colormap{3},'MarkerEdgeColor','none','Markersize',7); hold on;
    p_4 = plot(c(index_N+4,1), z(index_N+4,1),'o','MarkerFaceColor',colormap{4},'MarkerEdgeColor','none','Markersize',7); hold on;
    p_5 = plot(c(index_N+5,1), z(index_N+5,1),'o','MarkerFaceColor',colormap{5},'MarkerEdgeColor','none','Markersize',7); hold on;
    p_6 = plot(c(index_N+6,1), z(index_N+6,1),'o','MarkerFaceColor',colormap{6},'MarkerEdgeColor','none','Markersize',7); hold on;
    
    thelegend = legend([p_1 p_2 p_3 p_4 p_5 p_6],...
        {['u* = ',num2str(ustar(index_N+1),'%.3f'),' m/s, Q = ',num2str(Q(index_N+1),'%.2f'),' g/m/s'],...
        ['u* = ',num2str(ustar(index_N+2),'%.3f'),' m/s, Q = ',num2str(Q(index_N+2),'%.2f'),' g/m/s'],...
        ['u* = ',num2str(ustar(index_N+3),'%.3f'),' m/s, Q = ',num2str(Q(index_N+3),'%.2f'),' g/m/s'],...
        ['u* = ',num2str(ustar(index_N+4),'%.3f'),' m/s, Q = ',num2str(Q(index_N+4),'%.2f'),' g/m/s'],...
        ['u* = ',num2str(ustar(index_N+5),'%.3f'),' m/s, Q = ',num2str(Q(index_N+5),'%.2f'),' g/m/s'],...
        ['u* = ',num2str(ustar(index_N+6),'%.3f'),' m/s, Q = ',num2str(Q(index_N+6),'%.2f'),' g/m/s']});
    set(thelegend,'location','northeast','FontSize',9,'FontWeight','bold');
    
    % plot setting
    set(gca,'xscale','log');
    set(gca,'yscale','log');
    ylim([0.7,7]);    
    xlabel('Aerosol mass concentration \itC_{m}\rm\bf, kg/m^{3}','FontSize',12,'FontWeight','bold');
    ylabel('Sensor height \itz\rm\bf, m','FontSize',12,'FontWeight','bold');
    title(['Bin - ',num2str(N),' (figure #',num2str(i),' of total ',num2str(N_fig+1),' ones)'],'FontSize',12,'FontWeight','bold');    
    print([path,'ConcenProfile_Bin_',num2str(N),'_',num2str(i),'.png'],'-dpng'); 
    clf;
end

if N_line > 0
figure(N_fig+1);
index_N = N_fig*6;
for j = 1:N_line
    plot(c(index_N+j,:), z(index_N+j,:),'o','MarkerFaceColor',colormap{j},'MarkerEdgeColor','none','Markersize',7); hold on;
    plot(c_fit(index_N+j,:), z(index_N+j,:),'-','Color',colormap{j},'linewidth',1); hold on;
    for k = 1:4
        tempx = [c(index_N+j,k)-c_err(index_N+j,k), c(index_N+j,k)+c_err(index_N+j,k)];
        tempy = [z(index_N+j,k), z(index_N+j,k)];
        plot(tempx, tempy, '-','color',colormap{j},'linewidth',1); hold on;
    end
end
% plot legend
p_1 = plot(c(index_N+1,1), z(index_N+1,1),'o','MarkerFaceColor',colormap{1},'MarkerEdgeColor','none','Markersize',7); hold on;
p_2 = plot(c(index_N+2,1), z(index_N+2,1),'o','MarkerFaceColor',colormap{2},'MarkerEdgeColor','none','Markersize',7); hold on;
p_3 = plot(c(index_N+3,1), z(index_N+3,1),'o','MarkerFaceColor',colormap{3},'MarkerEdgeColor','none','Markersize',7); hold on;
p_4 = plot(c(index_N+4,1), z(index_N+4,1),'o','MarkerFaceColor',colormap{4},'MarkerEdgeColor','none','Markersize',7); hold on;
p_5 = plot(c(index_N+5,1), z(index_N+5,1),'o','MarkerFaceColor',colormap{5},'MarkerEdgeColor','none','Markersize',7); hold on;

thelegend = legend([p_1 p_2 p_3 p_4 p_5],...
    {['u* = ',num2str(ustar(index_N+1),'%.3f'),' m/s, Q = ',num2str(Q(index_N+1),'%.2f'),' g/m/s'],...
    ['u* = ',num2str(ustar(index_N+2),'%.3f'),' m/s, Q = ',num2str(Q(index_N+2),'%.2f'),' g/m/s'],...
    ['u* = ',num2str(ustar(index_N+3),'%.3f'),' m/s, Q = ',num2str(Q(index_N+3),'%.2f'),' g/m/s'],...
    ['u* = ',num2str(ustar(index_N+4),'%.3f'),' m/s, Q = ',num2str(Q(index_N+4),'%.2f'),' g/m/s'],...
    ['u* = ',num2str(ustar(index_N+5),'%.3f'),' m/s, Q = ',num2str(Q(index_N+5),'%.2f'),' g/m/s']});
set(thelegend,'location','northeast','FontSize',9,'FontWeight','bold');

% plot setting
set(gca,'xscale','log');
set(gca,'yscale','log');
ylim([0.7,7]);
xlabel('Aerosol mass concentration \itC_{m}\rm\bf, kg/m^{3}','FontSize',12,'FontWeight','bold');
ylabel('Sensor height \itz\rm\bf, m','FontSize',12,'FontWeight','bold');
title(['Bin - ',num2str(N),' (figure #',num2str(N_fig+1),' of total ',num2str(N_fig+1),' ones)'],'FontSize',12,'FontWeight','bold');   
       
print([path,'ConcenProfile_Bin_',num2str(N),'_',num2str(N_fig+1),'.png'],'-dpng'); 
clf;
end

end