function PlotBulkAlpha(Struct, bin_N, D_dust_range, D_rou)
% PlotBulkAlpha(Struct, bin_N, D_dust_range, D_rou)
% [b, SE_b] = PlotAlpha(Struct, N, D_dust_range, D_rou), if output value of
% power and its uncertainty
% plot figures for the bulk sandblasting efficiency summed over first bin_N
% number of bins
% plot the data in exponential scale to better represent the errorbars, and
% include the bin's geometric "D_dust_range" and aerodynamic diamter range

% convert dust geometric diameter to aerodynamic diameter
[Dae_dust_range] = GeometricD2AerodynamicD(D_dust_range, D_rou);

%% calculate bulk alpha
Q = Struct.Q*10^(-3); % saltation flux, in unit of kg/m/s
Q_err = Struct.sigma_Q*10^(-3); % in unit of kg/m/s
[Fd, Fd_err] = CalculateBulkMassFlux(Struct, bin_N); % Fd, and Fd_err are in units of kg/m^2/s
bulk_alpha = Fd./Q; % sandblasting efficient, in unit of m^(-1)
bulk_alpha_err = sqrt((-1*Fd./Q./Q.*Q_err).^2 + (Fd_err./Q).^2); % in unit of m^(-1)     

%% arrange data in the order of increasing shear velocity
temp_x = Struct.ustar;
[x, index_ustar] = sort(temp_x,'ascend');

y = bulk_alpha(index_ustar);
y_err = bulk_alpha_err(index_ustar);

%% power fit on the sandblasting efficient as a function of shear velocity
x_power = log10(x);
y_power = log10(y);
yerr_power = sqrt((y_err./y/log(10)).^2); % error propagation

[a,b,SE_a,SE_b,cov_ab2] = Linearfit(x_power, y_power, yerr_power);

y_fit_power = a + b*x_power; % fitted y under log10 scale
yerr_fit_power = sqrt(SE_a^2 + SE_b^2*x_power.^2 + 2*cov_ab2*x_power); % error of fitted y under log10 scale

y_fit = 10.^y_fit_power;
yerr_fit = sqrt((10.^y_fit_power*log(10).*yerr_fit_power).^2);

%% plot data under log10 scale
plot(x_power, y_power,'o','MarkerFaceColor','none','MarkerEdgeColor','k','Markersize',7,'linewidth',1); hold on;

%% plot error bar under log10 scale
for i = 1:length(x)
    tempx = [x_power(i), x_power(i)];
    tempy = [y_power(i)-yerr_power(i), y_power(i)+yerr_power(i)];
    plot(tempx,tempy,'k-'); hold on;
end

%% plot fit lines under log10 scale
se_fit_p = plot(x_power, y_fit_power,'k-'); hold on;
se_fitdn_p = plot(x_power, y_fit_power-yerr_fit_power,'k:'); hold on;
se_fitup_p = plot(x_power, y_fit_power+yerr_fit_power,'k:'); hold on;

%% plot legend ---
% thethe = legend(se_fit_p,['\alpha \propto u_{*}^{',num2str(b,'%.2f'),'}']);
thethe = legend(se_fit_p,['\alpha \propto u_{*}^{',num2str(b,'%.2f'),'+/-',num2str(SE_b,'%.2f'),'}']);
xlabel('Shear velocity \itu_{*}\rm\bf, m/s','FontSize',12,'FontWeight','bold');

%% plot setting ---
xlim([log10(0.29),log10(0.43)]);
set(gca,'xtick',log10(0.29:0.01:0.43));
set(gca,'xticklabel',{'0.29','0.3','','','','','0.35','','','','','0.4','','','0.43'});
    
ylim([-6.4,-5.4]);
set(gca,'ytick',-6.4:0.1:-5.4);
set(gca,'yticklabel',{'6E-7','','','','10^{-6}','','','','','','7E-6'});
title(['\itD_{g}\rm\bf = ',num2str(D_dust_range(1),'%.2f'),' - ',num2str(D_dust_range(bin_N+1),'%.2f'),' \mum, ',...
    '\itD_{a}\rm\bf = ',num2str(Dae_dust_range(1),'%.2f'),' - ',num2str(Dae_dust_range(bin_N+1),'%.2f'),' \mum'],...
    'FontSize',12,'FontWeight','bold'); 
ylabel('Sandblasting efficienty \it\alpha = F_{d}/Q\rm\bf, m^{-1}','FontSize',12,'FontWeight','bold');
set(thethe,'location','northwest','FontSize',12,'FontWeight','bold');

end

