% Integrate modeled backscattering efficiencies values (Qbb) 
% to compute the modeled backscattering coefficinets, using band-averaged
% Qbb at common bands of multispectral sensors. 
%
% (c)Tihomir Kostadinov, California State University San Marcos 
% 2019-2021

base_dir = '';

input_Qs_fname = fullfile(base_dir, 'MonteCarlo_Qbb_Bandavg_EAP.mat');
Qs = load(input_Qs_fname);

input_Qs_fname_det = fullfile(base_dir, 'MonteCarlo_Qbb_Bandavg_NAP.mat');
Qs_det = load(input_Qs_fname_det);

J_slope = (2.5:0.05:6)';

wave_all_sensors = Qs.wave_all_sensors;

Do = 2e-6; %meters
No_phi = 5e+16;
No_det = 2*No_phi; %This key assumption results in the assumption that phytoC = (1/3)*POC, 
% if the limits of integration do not change between POC and phyto C estimation. 
%No_total = No_phi + No_det; %No_total is the parameter retrieved by the inverse PSD model

if isempty(gcp('nocreate'))
    parpool(4);
end

[bbp_phi, bbp_det] =deal(NaN(size(Qs.Qbb_bandavg,3),numel(Qs.Qbb_bandavg,2),numel(J_slope)));

tic
for mm = 1:size(Qs.Qbb_bandavg,3) %MC simulations loop
    
    Qbb = Qs.Qbb_bandavg(:,:,mm);
    D = Qs.D_for_int(:,mm);
    
    Qbb_det = Qs_det.Qbb_bandavg_homogen(:,:,mm);
    D_det = Qs_det.D_for_int(:,mm);
    
    for k = 1:size(Qbb,2) %Wavelength loop
        if ~mod(mm,5)
            [mm wave_all_sensors(k)]
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% To limit integration to >= 0.5 um only! %%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Qbbl05 = Qbb(:,k);
        Dl05 = D;
        qq = D>=0.5e-6;
        Dl05(~qq) = [];
        Qbbl05(~qq) = [];
                
        Qbbl_det = Qbb_det(:,k);
        Dl_det = D_det;

        Dl05(isnan(Qbbl05)) = [];
        Qbbl05(isnan(Qbbl05)) = [];
        
        Dl_det(isnan(Qbbl_det)) = [];
        Qbbl_det(isnan(Qbbl_det)) = [];
        
        %Integrate Qbb & PSD to get modeled bbp for phytoplankton, and for NAP
        %See, for example, Kostadinov et al. (2009) 
        parfor jj = 1:numel(J_slope)
            xi = J_slope(jj);
            bbp_phi(mm,k,jj) = (pi/4)*Do^xi*No_phi.*trapz(Dl05,Qbbl05.*(Dl05.^(2-xi)));
            bbp_det(mm,k,jj) = (pi/4)*Do^xi*No_det.*trapz(Dl_det,Qbbl_det.*(Dl_det.^(2-xi)));
        end
    end
end

bbp_tot = bbp_phi + bbp_det; 
bbp_phi_frac = bbp_phi./bbp_tot;

etime = toc;
mfile=mfilename('fullpath');

save('-mat','-v7.3',fullfile(base_dir,'bbp_2pop_from_bandavg_Qs_allMC.mat'),...
    'bbp*','J_slope','input_Qs_fname*','No_phi','No_det','wave_all_sensors','mfile','etime');