function [Flux] = NumConcenProfile2Flux(Struct, Bin_D, rou, N)
% [Flux] = NumConcenProfile2Flux(Struct, Bin_D, rou, N)
% Input a number concentration, OPC representative bin size, dust density,
% and which sensors of the 4 OPCs used
% Output the mass flux (unit: kg/m^2/s) using the input sensors
% calculation based on method in Kind (1992)


karman = 0.387; % von karman constant, Andreas et al. (2006) from Journal of Fluid Mechanics 

Fd = [];
Fd_err = [];

%--- gradient method (Kind, 1992) to obtain the mass flux
for i = 1:7
   temp_y = Struct.(['bin',num2str(i)]).y(:,N(1):N(2));
   temp_y_err = Struct.(['bin',num2str(i)]).y_err(:,N(1):N(2));
   
   y = NumberConcen2MassConcen(temp_y, Bin_D(i), rou); % Mass concentration on kg/m^3
   y_err = NumberConcen2MassConcen(temp_y_err, Bin_D(i), rou); % Mass concentration on kg/m^3
   
   temp_z0 = Struct.z(:,N(1));
   temp_z = Struct.z(:,N(1):N(2)); 
   
   for j = 1:size(temp_z0,1)
       temp_z2z0(j,:) = temp_z(j,:)/temp_z0(j);
   end
   
   z_power = log(temp_z2z0); %log(z_i/z_r)
   
   ustar = Struct.SaltationSand.ustRe; %shear velocity, in unit of m/s
   ustar_err = Struct.SaltationSand.sigma_ustRe; %shear velocity, in unit of m/s
   
   for j = 1:size(temp_z0,1)
       [~, b, ~, sigma_b, ~] = Linearfit(z_power(j,:)', y(j,:)',y_err(j,:)');
       fd(j,1) = -1*b*karman*ustar(j); % Mass concentration, in the unit of kg/m^2/s, from Kind (1992)
       fd_err(j,1) = sqrt((-1*sigma_b*karman*ustar(j)).^2 + (-1*b*karman*ustar_err(j)).^2); % from error propogation
   end
   
   Fd = [Fd, fd];
   Fd_err = [Fd_err, fd_err];   
end
% --- obtain sum of the 7 bins' mass flux and store it in the Fd(:,8)
for j = 1:size(Fd,1)
    Fd(j,8) = sum(Fd(j,1:7)); 
    Fd_err(j,8) = sqrt(sum(Fd_err(j,1:7).^2)); % error propogation
end
%--- store the matrix---

Flux.Fd = Fd;
Flux.Fd_err = Fd_err;
Flux.Fd_unit = 'kg/m^2/s';
Flux.z = temp_z; %store z, in unit of m/s

Flux.StartTimes = Struct.SaltationSand.StartTimes;
Flux.EndTimes = Struct.SaltationSand.EndTimes;
Flux.Q = Struct.SaltationSand.Q;
Flux.sigma_Q = Struct.SaltationSand.sigma_Q;
Flux.ustRe = Struct.SaltationSand.ustRe;
Flux.sigma_ustRe = Struct.SaltationSand.sigma_ustRe;
Flux.theta = Struct.SaltationSand.theta;

display('%--- function ConcenProfile2Flux() complete ----%');

end