function PlotRemoveSeaSalt(Struct_Salta_fit, Struct_NoSalta_fit, N, type)
% PlotRemoveSeaSalt(Struct_fit, Struct_Salta, Struct_NoSalta, N)
% plot figures to diagnose the sea salt removal process, for bin N
    x_raw_s = Struct_Salta_fit.(['bin',num2str(N)]).ustar;
    
    y_raw_s = Struct_Salta_fit.(['bin',num2str(N)]).Fd_raw*10^9; % ug/m^2/s
    yerr_raw_s = Struct_Salta_fit.(['bin',num2str(N)]).Fd_err_raw*10^9; % ug/m^2/s    
    y_fit_s = Struct_Salta_fit.(['bin',num2str(N)]).Fd_fit*10^9; % ug/m^2/s
    yerr_fit_s = Struct_Salta_fit.(['bin',num2str(N)]).Fd_err_fit*10^9; % ug/m^2/s
    
    x_raw_n = Struct_NoSalta_fit.(['bin',num2str(N)]).ustar;   
    
    y_raw_n = abs(Struct_NoSalta_fit.(['bin',num2str(N)]).Fd_raw)*10^9; % ug/m^2/s
    yerr_raw_n = Struct_NoSalta_fit.(['bin',num2str(N)]).Fd_err_raw*10^9; % ug/m^2/s
    y_fit_n = Struct_NoSalta_fit.(['bin',num2str(N)]).Fd_fit*10^9; % ug/m^2/s  %fit the fitted absolute value of deposition flux
    yerr_fit_n = Struct_NoSalta_fit.(['bin',num2str(N)]).Fd_err_fit*10^9; % ug/m^2/s
    
    %--- order series with increasing ustar order, for saltation occurance
    [x_raw_s_sort, index_salta] = sort(x_raw_s,'ascend');
    
    y_raw_s_sort = y_raw_s(index_salta);
    yerr_raw_s_sort = yerr_raw_s(index_salta);
    y_fit_s_sort = y_fit_s(index_salta);
    yerr_fit_s_sort = yerr_fit_s(index_salta);
    
    %--- order series with increasing ustar order, for saltation occurance    
    [x_raw_n_sort, index_nosalta] = sort(x_raw_n,'ascend');
    
    y_raw_n_sort = y_raw_n(index_nosalta);
    yerr_raw_n_sort = yerr_raw_n(index_nosalta);
    y_fit_n_sort = y_fit_n(index_nosalta);
    yerr_fit_n_sort = yerr_fit_n(index_nosalta);
    
    %--- plot saltation and nosaltation onto the same figure for each bin N
    raw_s = plot(x_raw_s_sort, y_raw_s_sort,'<','markerfacecolor','r','markeredgecolor','none','markersize',4); hold on;
    raw_n = plot(x_raw_n_sort, y_raw_n_sort,'o','markerfacecolor','none','markeredgecolor',[0,0.6,0.8],'markersize',4); hold on;
    
    fit_s = plot(x_raw_s_sort, y_fit_s_sort,'r-','Linewidth',1.5); hold on;
    fit_dn_s = plot(x_raw_s_sort, y_fit_s_sort - yerr_fit_s_sort,'r:','Linewidth',1.5); hold on;
    fit_up_s = plot(x_raw_s_sort, y_fit_s_sort + yerr_fit_s_sort,'r:','Linewidth',1.5); hold on;
    
    fit_n = plot(x_raw_n_sort, y_fit_n_sort,'-','color',[0,0.6,0.8],'Linewidth',1.5); hold on;
    fit_dn_n = plot(x_raw_n_sort, y_fit_n_sort - yerr_fit_n_sort,':','color',[0,0.6,0.8],'Linewidth',1.5); hold on;
    fit_up_n = plot(x_raw_n_sort, y_fit_n_sort + yerr_fit_n_sort,':','color',[0,0.6,0.8],'Linewidth',1.5); hold on;    
    
    %--- plot error bar
    for i = 1:length(x_raw_s_sort)
        tempx = [x_raw_s_sort(i),x_raw_s_sort(i)];
        tempy = [y_raw_s_sort(i), y_raw_s_sort(i) + yerr_raw_s_sort(i)];
        plot(tempx, tempy, 'r-'); hold on;
    end
    
    for i = 1:length(x_raw_n_sort)
        tempx = [x_raw_n_sort(i), x_raw_n_sort(i)];
        tempy = [y_raw_n_sort(i), y_raw_n_sort(i) + yerr_raw_n_sort(i)];
        plot(tempx, tempy,'-','color',[0,0.6,0.8]); hold on; 
    end
    
    set(gca,'xscale','log');
    set(gca,'yscale','log');
    
    xlim([0.07, 0.6]);
    xlabel('Shear velocity \itu_{*}\rm\bf, m/s','FontSize',12,'FontWeight','bold');
    ylabel('Vertical aerosol mass flux \itF_{d}\rm\bf, \mug/m^{2}/s','FontSize',12,'FontWeight','bold');
    title(['Bin - ',num2str(N)],'FontSize',12,'FontWeight','bold');   
    
    thelegend = legend([raw_s raw_n],'Q > 0, F_{d} > 0','Q = 0, F_{d} < 0');
    set(thelegend,'location','northwest','FontSize',12,'FontWeight','bold');
end