function [betap_corr, bbp_corr, flags] = betap_Kcorrection(beta_total, beta_filt_interp, ...
  lambda_beta, k_exp, Chi, ac_p, ac_g, lambda_a, lambda_c, fdom_ag_correlation)
  % This function performs the attenuation correction on beta measured with backscattering sensors
  % This method is based on the HS4 correction devised by Doxaran et al., 2016 (10.1364/OE.24.003615)
  % but rather than iterating, we expand the exponential in that paper in a two terms Taylor series.
  % Author: Emmanuel Boss and Guillaume Bourdin
  % Date: 2024-08-29
  %
  % Inputs:
  % - beta_total: uncalibrated measured beta total
  % - beta_filt_interp: uncalibrated measured beta filtered
  % - lambda_beta: beta wavelength
  % - k_exp: optical pathlength
  % - Chi: factor needed to convert the particulate VSF to particulate backscattering from Sullivan et al. 2013
  % - ac_p: measured calibrated absorption and attenuation of particulate matter
  % - ac_g: measured calibrated absorption and attenuation of dissolved matter
  % - lambda_a: absorption channel wavelength
  % - lambda_c: absorption channel wavelength
  % - fdom_ag_correlation (option) if ag and cg are not compute from underway data
  %
  %%
  if nargin < 10
    % TODO add case when fdom ag/cg is computed out of inlineanalysis
  end
  % TODO add salinity correction forwarding TSG to the function
  
  % round all timestamps to flag interpolated data
  beta_total = round_timestamp(beta_total, minutes(1));
  beta_filt_interp = round_timestamp(beta_filt_interp, minutes(1));
  ac_p = round_timestamp(ac_p, minutes(1));
  ac_g = round_timestamp(ac_g, minutes(1));
  
  %% create flag table
  betap_approx = beta_total.beta - beta_filt_interp.beta;
  bbp_approx = 2 * pi * Chi .* betap_approx;
  flags = table(false(size(beta_total.dt)), false(size(beta_total.dt)), false(size(beta_total.dt)), ...
    false(size(beta_total.dt)), false(size(beta_total.dt)), false(size(beta_total.dt)), ...
    false(size(beta_total.dt)), false(size(beta_total.dt)), false(size(beta_total.dt)), false(size(beta_total.dt)), ...
    'VariableNames', {'no_attenuation_correction', 'no_ag_attenuation_correction', ...
    'simplified_Doxaran_correction', 'ap_interpolated', 'cp_interpolated', ...
    'ag_interpolated', 'cg_interpolated', 'ag_input_interpolated_linearly', ...
    'cg_input_interpolated_linearly','fdom_input_interpolated_linearly'});

  % Format ap, cp and compute bp if particulate data available
  if ~isempty(ac_p)
    [ap_interp, flag.ap_interpolated] = interp_on_bb_lambda(ac_p.dt, ac_p.ap, lambda_a, beta_total.dt, lambda_beta);
    [cp_interp, flag.cp_interpolated] = interp_on_bb_lambda(ac_p.dt, ac_p.cp, lambda_c, beta_total.dt, lambda_beta);
    % compute bp
    bp_interp = cp_interp - ap_interp;
  end

  % Format ag and cg if dissolved data available, otherwise use a function derived independently
  if ~isempty(ac_g)
    % rebuild ag from exponential fit parameters
    ag_rebuilt = NaN(size(ac_g.ag));
    cg_rebuilt = NaN(size(ac_g.cg));
    % define exponential function
    expfun = @(p, xd) p(1) * exp(p(2) * (xd - 440));
    for i = 1:size(ac_g, 1)
      ag_rebuilt(i, :) = expfun([ac_g.y_intercp_fit_ag(i) ac_g.base_fit_ag(i)], lambda_a);
      cg_rebuilt(i, :) = expfun([ac_g.y_intercp_fit_cg(i) ac_g.base_fit_cg(i)], lambda_c);
    end
    % merge ag and rebuilt ag
    ag_merged = ac_g.ag;
    ag_merged(isnan(ag_merged)) = ag_rebuilt(isnan(ag_merged));
    ag_merged(all(isnan(ac_g.ag), 2), :) = NaN;
    % merge cg and rebuilt cg
    cg_merged = ac_g.cg;
    cg_merged(isnan(cg_merged)) = cg_rebuilt(isnan(cg_merged));
    cg_merged(all(isnan(ac_g.ag), 2), :) = NaN;

    % Reconstruct ag and cg based on correlation with fdom if available otherwise interpolate linearly between filter events
    if any(strcmp(ac_p.Properties.VariableNames, 'fdom')) && any(strcmp(ac_g.Properties.VariableNames, 'fdom'))
      [beta_total.fdom, flags.fdom_input_interpolated_linearly] = interp_on_bb_lambda(ac_p.dt, ac_p.fdom, [], beta_total.dt, []);
      [a_reg, c_reg] = regress_ac_fdom(ag_merged, cg_merged, ac_g.fdom);
      % interpolate slope on lambda_beta
      a_slope_interp = interp1(lambda_a, a_reg.slope, lambda_beta, 'linear');
      c_slope_interp = interp1(lambda_c, c_reg.slope, lambda_beta, 'linear');
      % linear interpolation of a, c, and fdom
      lin_ag = interp_on_bb_lambda(ac_g.dt, ag_merged, lambda_a, beta_total.dt, lambda_beta);
      lin_cg = interp_on_bb_lambda(ac_g.dt, cg_merged, lambda_c, beta_total.dt, lambda_beta);
      lin_fdom = interp_on_bb_lambda(ac_g.dt, ac_g.fdom, [], beta_total.dt, []);
      % Re-build ag and cg interpolated between filter events
      ag_interp = lin_ag + (beta_total.fdom - lin_fdom) .* a_slope_interp; 
      cg_interp = lin_cg + (beta_total.fdom - lin_fdom) .* c_slope_interp;
      
      % visProd3D(lambda_beta, beta_total.dt, beta_total.ag_interp, false, 'Wavelength', false, 2);
      % title('ag fdom interpolated')
      % visProd3D(lambda_beta, beta_total.dt, beta_total.cg_interp, false, 'Wavelength', false, 3);
      % title('cg fdom interpolated')
    
      % Fill missing values (start/end of analyzed period) with nearest filt_interp data
      ag_interp = fillmissing(ag_interp, 'linear', 'SamplePoints', beta_total.dt, 'EndValues', 'nearest');
      cg_interp = fillmissing(cg_interp, 'linear', 'SamplePoints', beta_total.dt, 'EndValues', 'nearest');
    else
      % interpolate ag on BB lambda
      [ag_interp, flag.ag_interpolated] = interp_on_bb_lambda(ac_g.dt, ag_merged, lambda_a, beta_total.dt, lambda_beta);
      % interpolate cg on BB lambda
      [cg_interp, flag.cg_interpolated] = interp_on_bb_lambda(ac_g.dt, cg_merged, lambda_c, beta_total.dt, lambda_beta);
      flags.ag_input_interpolated_linearly = true;
      flags.cg_input_interpolated_linearly = true;
    end

  elseif ~isempty(fdom_ag_correlation)
    % Case when the correlation between fdom and ag is computed from data from a benchtop spectrophotometer
    error('Not implemented')
  end
  
  
  if isempty(ac_p)
    % Just use beta_total - beta_filt as betap and flag data
    warning('No ap data input for attenuation correction of beta: attenuation correction not applied')
    flags.no_attenuation_correction = true(size(beta_total, 1), 1);
    betap_corr = betap_approx;
    bbp_corr = bbp_approx;
  elseif isempty(ac_g)
    % Approximate ag if dissolved not available and apply simplified Doxaran's correction and flag data
    warning('No ag data input for attenuation correction of beta: attenuation correction applied without ag')
    flags.no_ag_attenuation_correction = true(size(beta_total, 1), 1);
    flags.simplified_Doxaran_correction = true(size(beta_total, 1), 1);
    % TODO estimate ag from ap or chlorophyll
    error('Not implemented')
    ag = 99999999999;
    % Adapted Doxaran's method to work with only ap and bp as input
    k_total = ag + ap_interp + 3.3 .* bp_interp;
    % apply correction
    betap_corr = (beta_total.beta - beta_filt_interp.beta) ./ (1 - k_total * k_exp);
    bbp_corr = 2 * pi * Chi .* betap_corr;
    
    % visProd3D(lambda_beta, beta_total.dt, k_total, false, 'Wavelength', false, 72);
    % title('k total')

    % betap_notcorr = beta_total.beta - beta_filt_interp.beta;
    % visProd3D(lambda_beta, beta_total.dt, betap_notcorr, false, 'Wavelength', false, 73);
    % title('betap approx')

    % visProd3D(lambda_beta, beta_total.dt, bbp_corr, false, 'Wavelength', false, 74);
    % title('bbp corr')

    % visProd3D(lambda_beta, beta_total.dt, (bbp_corr-betap_notcorr)./betap_notcorr*100, false, 'Wavelength', false, 75);
    % title('% change Doxaran')

  else
    % apply Zhang's sigma correction if ag and ap are available  
    load('Zhang_hyper_bb_sigma_correction.mat', 'B', 'bbr', 'c', 'w')
    % correct beta_filt
    beta_filt_corr = beta_filt_interp.beta .* exp(ag_interp * k_exp);
    % Estimate c
    c_estimated = cp_interp + cg_interp;
    % Estimate single scattering albedo (w)
    w_estimated = bp_interp ./ c_estimated;
    % Force w to be between 0 and 1 when 1 < w_estimated < 1.005 (uncertainty of ACS == 0.005 m^-1)
    w_estimated(w_estimated > 1 & w_estimated <= 1.005) = 1;
    
    % Find sigma correction by iteration starting with betap approximation for backscattering ratio computation
    keep_going = true;
    bbp_approx_change = NaN(5, size(beta_total.beta, 2));
    bbp_approx_new = bbp_approx;
    i = 1;
    while keep_going
      % Estimate backscattering ratio
      bbr_estimated = bbp_approx_new ./ bp_interp ;
      % isolate bbp_approx to compare after iteration
      bbp_approx_old = bbp_approx_new;
      % Find sigma in Zhang lookup table
      sigma_corr = interp3(bbr, w, c, B, bbr_estimated, w_estimated, c_estimated, "linear");
      sigma_corr = fillmissing(sigma_corr, 'nearest', 'SamplePoints', datetime(beta_total.dt, 'ConvertFrom', 'datenum'), 'MaxGap',minutes(40));
      betap_approx_new = beta_total.beta .* sigma_corr - beta_filt_corr;
      bbp_approx_new = 2 * pi * Chi .* betap_approx_new;
      bbp_approx_change(i, :) = mean((bbp_approx_new - bbp_approx_old) ./ bbp_approx_old, 1, 'omitmissing');
      if max(bbp_approx_change(i, :)) < 0.000001 || i == size(bbp_approx_change, 1)
        keep_going = false;
        betap_corr = betap_approx_new;
        bbp_corr = bbp_approx_new;
      else
        i = i +1;
      end
    end

    % visProd3D(lambda_beta, beta_total.dt, beta_total.ag_interp, false, 'Wavelength', false, 68);
    % zlabel('a_g input')
    % 
    % visProd3D(lambda_beta, beta_total.dt, beta_total.cg_interp, false, 'Wavelength', false, 69);
    % zlabel('c_g input')
    % 
    % visProd3D(lambda_beta, beta_total.dt, beta_total.ap_interp, false, 'Wavelength', false, 70);
    % zlabel('a_p input')
    % 
    % visProd3D(lambda_beta, beta_total.dt, beta_total.cp_interp, false, 'Wavelength', false, 71);
    % zlabel('c_p input')
    % 
    % visProd3D(lambda_beta, beta_total.dt, bbp_approx, false, 'Wavelength', false, 72);
    % zlabel('bbp approximation (2 * pi * Chi .* (beta_{total} - beta_{filt}))'); title('uncorrected')
    % 
    % visProd3D(lambda_beta, beta_total.dt, sigma_corr, false, 'Wavelength', false, 73);
    % zlabel('sigma correction'); title("Zhang's lookup table")
    % 
    % visProd3D(lambda_beta, beta_total.dt, bbp_corr, false, 'Wavelength', false, 74);
    % zlabel('bbp corrected'); title("Zhang's lookup table")
    % 
    % visProd3D(lambda_beta, beta_total.dt, (bbp_approx_new - bbp_approx) ./ bbp_approx*100, false, 'Wavelength', false, 75);
    % zlabel('% change bbp'); title("Zhang's lookup table")

  end
end

function [data_interp, flag] = interp_on_bb_lambda(data_dt, data, lambda, bb_dt, lambda_bb)
  flag = ~ismember(bb_dt, data_dt);
  % interpolate ap on BB lambda
  if ~isempty(lambda)
    data_bbwl = interp1(lambda, data', lambda_bb, 'linear')';
  else
    data_bbwl = data;
  end
  % interpolate ag on beta_total datetime
  data_interp = interp1(data_dt, data_bbwl, bb_dt, 'linear');
  data_interp = fillmissing(data_interp, 'nearest', 'SamplePoints', bb_dt);
end

%% regression between a/c and fdom
function [a_reg, c_reg] = regress_ac_fdom(a_filt, c_filt, fdom)
  % robust linear regression between filt a&c and fdom
  a_reg = array2table(NaN(size(a_filt, 2), 5), 'VariableNames', {'slope', 'intercept', 'RMSE', 'nRMSE', 'R2'});
  c_reg = array2table(NaN(size(c_filt, 2), 5), 'VariableNames', {'slope', 'intercept', 'RMSE', 'nRMSE', 'R2'});
  for i = 1:size(a_filt, 2)
    if sum(~isnan(a_filt(:, i))) > 3
      % regress variable with filtered water absorption
      [stats.b, stats.stats] = robustfit(fdom, a_filt(:, i));
      [~, MSGID] = lastwarn();
      warning('off', MSGID)
      a_reg.slope(i) = stats.b(2);
      a_reg.intercept(i) = stats.b(1);
      a_reg.RMSE(i) = stats.stats.robust_s;
      % calculated normalize RMSE by average and R2
      id_nonan = all(~isnan([a_filt(:, i) fdom]), 2);
      a_reg.nRMSE(i) = a_reg.RMSE(i)/abs(mean(a_filt(id_nonan, i), 'omitnan'))*100;
      a_reg.R2(i) = corr(a_filt(id_nonan, i), stats.b(1) + stats.b(2) * fdom(id_nonan))^2;
    end
    if sum(~isnan(c_filt(:, i))) > 3
      % regress variable with filtered water attenuation
      [stats.b, stats.stats] = robustfit(fdom, c_filt(:, i));
      [~, MSGID] = lastwarn();
      warning('off', MSGID)
      c_reg.slope(i) = stats.b(2);
      c_reg.intercept(i) = stats.b(1);
      c_reg.RMSE(i) = stats.stats.robust_s;
      % calculated normalize RMSE by average and R2
      id_nonan = all(~isnan([c_filt(:, i) fdom]), 2);
      c_reg.nRMSE(i) = c_reg.RMSE(i)/abs(mean(c_filt(id_nonan, i), 'omitnan'))*100;
      c_reg.R2(i) = corr(c_filt(id_nonan, i), stats.b(1) + stats.b(2) * fdom(id_nonan))^2;
    end
  end
end
