function [p, g, bad, regression_stats, DIW_biofouling_correction] = processACS(lambda, tot, filt, param, modelG50, modelmphi, di, ...
  CDOM, FTH, fth_constants, interpolation_method, di_method, scattering_correction, compute_ad_aphi, TSG, days2run)
  % NOTE: wavelength of c are interpolated to wavelength of a
  % load psi_s & psi_t from Sullivan et al. 2006 for salinity & temprature correction
  psi = table();
  psi.wl = [400;402;404;406;408;410;412;414;416;418;420;422;424;426;428;430;432;434;436;438;440;442;444;446;448;450;452;454;456;458;460;462;464;466;468;470;472;474;476;478;480;482;484;486;488;490;492;494;496;498;500;502;504;506;508;510;512;514;516;518;520;522;524;526;528;530;532;534;536;538;540;542;544;546;548;550;552;554;556;558;560;562;564;566;568;570;572;574;576;578;580;582;584;586;588;590;592;594;596;598;600;602;604;606;608;610;612;614;616;618;620;622;624;626;628;630;632;634;636;638;640;642;644;646;648;650;652;654;656;658;660;662;664;666;668;670;672;674;676;678;680;682;684;686;688;690;692;694;696;698;700;702;704;706;708;710;712;714;716;718;720;722;724;726;728;730;732;734;736;738;740;742;744;746;748;750];
  psi.psiT = [0.0001;0.0001;0.0001;0.0001;0;0;0;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0002;0.0003;0.0003;0.0004;0.0005;0.0006;0.0006;0.0007;0.0008;0.0009;0.001;0.001;0.001;0.001;0.001;0.0009;0.0009;0.0008;0.0007;0.0006;0.0006;0.0005;0.0004;0.0003;0.0003;0.0002;0.0001;0.0001;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0002;0.0002;0.0002;0.0001;0.0001;0.0001;0;0;-0.0001;-0.0001;-0.0001;-0.0001;-0.0001;-0.0001;0;0;0.0001;0.0002;0.0003;0.0005;0.0007;0.0009;0.0013;0.0017;0.0021;0.0026;0.0032;0.0038;0.0045;0.0054;0.0063;0.0073;0.0083;0.0094;0.0104;0.0113;0.0121;0.0128;0.0133;0.0136;0.0136;0.0133;0.0129;0.0124;0.0116;0.0107];
  psi.a_psiS = [3.0e-05;3.0e-05;3.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;3.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;0;0;0;0;0;0;0;0;0;0;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;0;1.0e-05;2.0e-05;3.0e-05;3.0e-05;4.0e-05;5.0e-05;5.0e-05;6.0e-05;6.0e-05;6.0e-05;6.0e-05;6.0e-05;6.0e-05;6.0e-05;5.0e-05;5.0e-05;5.0e-05;5.0e-05;4.0e-05;4.0e-05;4.0e-05;4.0e-05;3.0e-05;3.0e-05;3.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;1.0e-05;1.0e-05;1.0e-05;0;0;-1.0e-05;-2.0e-05;-3.0e-05;-4.0e-05;-6.0e-05;-7.0e-05;-8.0e-05;-9.0e-05;-0.00011;-0.00012;-0.00014;-0.00015;-0.00016;-0.00017;-0.00018;-0.00019;-0.00020;-0.00020;-0.00020;-0.00021;-0.00021;-0.00021;-0.00021;-0.00021;-0.00020;-0.00017;-0.00013;-8.0e-05;-1.0e-05;7.0e-05;0.00016;0.00026;0.00037;0.00046;0.00054;0.00061;0.00067];
  psi.c_psiS = [-1.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-3.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-4.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-5.0e-05;-4.0e-05;-3.0e-05;-3.0e-05;-2.0e-05;-1.0e-05;0;1.0e-05;1.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;2.0e-05;1.0e-05;1.0e-05;1.0e-05;1.0e-05;0;0;0;-1.0e-05;-1.0e-05;-1.0e-05;-1.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-2.0e-05;-3.0e-05;-4.0e-05;-5.0e-05;-6.0e-05;-6.0e-05;-8.0e-05;-9.0e-05;-0.00010;-0.00011;-0.00013;-0.00014;-0.00016;-0.00017;-0.00018;-0.00019;-0.00020;-0.00021;-0.00022;-0.00022;-0.00023;-0.00023;-0.00023;-0.00024;-0.00024;-0.00024;-0.00024;-0.00022;-0.00021;-0.00017;-0.00012;-6.0e-05;2.0e-05;0.00012;0.00022;0.00031;0.00041;0.00049;0.00056;0.00062];
  psi.sigma_psiT = [0.0002;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0;0;0;0;0.0001;0.0001;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0;0;0;0.0001;0;0;0;0;0;0;0;0;0;0;0;0;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0001;0.0002;0.0002;0.0003;0.0003;0.0004;0.0004;0.0004;0.0004;0.0005;0.0005;0.0006;0.0006;0.0007;0.0007;0.0007;0.0006;0.0005;0.0004;0.0003;0.0003;0.0004;0.0005;0.0006;0.0007;0.0008;0.0009];
  psi.c_sigma_psiS = [4e-005;4e-005;4e-005;4e-005;4e-005;4e-005;4e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;0;-1e-005;-2e-005;-3e-005;-4e-005;-6e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;2e-005;2e-005;2e-005;3e-005;3e-005;3e-005];
  psi.a_sigma_psiS = [3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;3e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;2e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;NaN;NaN;NaN;NaN;NaN;NaN;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;1e-005;2e-005;2e-005;2e-005;3e-005];
  
  %% ap & cp
  % check FTH data
  if ~exist('fth_constants', 'var')
    % Assume most recent FlowControl software
    SWITCH_FILTERED = 1;
    SWITCH_TOTAL = 0;
  else
    SWITCH_FILTERED = fth_constants.SWITCH_FILTERED;
    SWITCH_TOTAL = fth_constants.SWITCH_TOTAL;
  end
  % check scattering correction method
  if ~any(strcmp(scattering_correction, {'Zaneveld1994_proportional', 'Rottgers2013_semiempirical', 'ZaneveldRottgers_blended'}))
    warning('only scattering correction supported: %s', strjoin({'Zaneveld1994_proportional', 'Rottgers2013_semiempirical', 'ZaneveldRottgers_blended'}, ', '))
    error('%s residual temperature and scattering correction not supported', scattering_correction)
  end
  flow_data = FTH.qc.tsw;

  % round time stamp and remove time duplicates
  tot = round_timestamp(tot, minutes(1));
  filt = round_timestamp(filt, minutes(1));
  flow_data = round_timestamp(flow_data, minutes(1));

  % check if fCDOM data loaded if fCDOM interpolation
  if ~isempty(CDOM)
    cdom_tbl_name = fieldnames(CDOM.prod);
    if isempty(cdom_tbl_name)
      error('No fDOM prod data loaded')
    end
    % Require both CDOM & Switch position
    if ~isempty(CDOM.prod.(cdom_tbl_name{1}))
      cdom_base = CDOM.prod.(cdom_tbl_name{1});
      % smooth fDOM if fDOM sensor is WSCD
      if strcmp(CDOM.model, 'WSCD') || strcmp(CDOM.prefix, 'WSCD')
        ancillary = table();
        foo_dt = [datetime(min([tot.dt; filt.dt]), 'ConvertFrom', 'datenum') datetime(max([tot.dt; filt.dt]), 'ConvertFrom', 'datenum')];
        ancillary.dt = datenum(dateshift((foo_dt(1):minutes(1):foo_dt(2))','start','minute'));
        id_nan = isnan(cdom_base.fdom);
        % create new fdom table interpolating for gaps < 60 min
        ancillary.fdom = interp1(cdom_base.dt, cdom_base.fdom, ancillary.dt, 'linear');
        ancillary.fdom = fillmissing(ancillary.fdom, 'nearest');
        % smooth WSCD cdom data removing frequencies higher than 15min
        d = designfilt('lowpassiir', 'FilterOrder', 1, 'PassbandFrequency', 1/(60*15), 'SampleRate', 1/60);
        ancillary.fdom = filtfilt(d, ancillary.fdom);
        cdom_base.fdom = interp1(ancillary.dt, ancillary.fdom, cdom_base.dt, 'linear');
        cdom_base.fdom(id_nan) = NaN;
      end
      if ~any(cdom_base.dt >= min([tot.dt; filt.dt]) & cdom_base.dt <= max([tot.dt; filt.dt]))
        warning('fCDOM dates do not correspond to ACS dates: interpolation switched to "linear"')
        interpolation_method = 'linear';
      else
        % round time stamp and remove time duplicates
        cdom_base = round_timestamp(cdom_base, minutes(1));
      end
    else
      warning('fCDOM data not loaded: interpolation switched to "linear"')
      interpolation_method = 'linear';
    end
    if ~isfield(param, 'min_nb_pts_per_cluster')
      fprintf('Minimum number of points per cluster set to 200\n')
      param.min_nb_pts_per_cluster = 200;
    end
  else
    cdom_base = [];
  end
  % check if TSG data loaded
  if ~isempty(TSG)
    tsg_tbl_name = fieldnames(TSG.prod);
    if ~isempty(tsg_tbl_name)
      if ~isempty(TSG.prod.(tsg_tbl_name))
        tsg_loaded = true;
        tsg_data = TSG.prod.(tsg_tbl_name);
      else
        tsg_loaded = false;
      end
    else
      tsg_loaded = false;
    end
    if ~tsg_loaded
      if ~isempty(TSG.qc.tsw)
        tsg_loaded = true;
        tsg_data = TSG.qc.tsw;
      else
        tsg_loaded = false;
      end
    end
    if tsg_loaded
      if ~any(tsg_data.dt >= min([tot.dt; filt.dt]) & tsg_data.dt <= max([tot.dt; filt.dt]))
        warning(['TSG dates do not correspond to ACS dates: T/S correction not applied before removing dissolved part from total.' newline ...
          'T/S assumed to be constant between filter and total events.' newline ...
          'Remaining T variation will be corrected using the residual temperature correction\n'])
        tsg_loaded = false;
      else
        tsg_loaded = true;
        % round time stamp and remove time duplicates
        tsg_data = round_timestamp(tsg_data, minutes(1));
      end
    end
  else
    tsg_loaded = false;
  end
  if ~tsg_loaded
    Warning('No TSG qc or prod data loaded')
  end

  % interpolate flow_data.swt onto binned data to fill missing flow data
  fth_interp = table([tot.dt; flow_data.dt; filt.dt], 'VariableNames', {'dt'});
  fth_interp = round_timestamp(fth_interp, minutes(1));
  fth_interp = sortrows(fth_interp, 'dt'); % sort dates
  fth_interp.swt = interp1(flow_data.dt, flow_data.(FTH.view.swt_variable), fth_interp.dt, 'previous');%, 'linear', 'extrap');
  fth_interp.swt = fth_interp.swt > 0;
  % Find switch events from total to filtered
  sel_start = find(fth_interp.swt(1:end-1) == SWITCH_TOTAL & fth_interp.swt(2:end) == SWITCH_FILTERED);
  % Find switch events from filtered to total
  sel_end = find(fth_interp.swt(1:end-1) == SWITCH_FILTERED & fth_interp.swt(2:end) == SWITCH_TOTAL);
  % Verify selections of filtered period
  if sel_start(1) > sel_end(1); sel_end(1) = []; end
  if sel_start(end) > sel_end(end); sel_end(end+1) = size(fth_interp.swt,1); end
  if size(sel_start,1) ~= size(sel_end,1); error('Inconsistent switch start/end data'); end

  %%%%%%%%%%%%%%%%%%%%%%% Compute filtered average %%%%%%%%%%%%%%%%%%%%%%%
  % TODO maybe: write a filter_average function that can be use with every instrument

  % create filter average data table
  filt_avg = table(NaN(size(sel_start)), 'VariableNames', {'dt'});
  filt_avg.start = fth_interp.dt(sel_start);
  filt_avg.end = fth_interp.dt(sel_end);
  filt_avg.a = NaN(size(filt_avg,1), size(lambda.a, 2));
  filt_avg.c = NaN(size(filt_avg,1), size(lambda.c, 2));
  filt_avg.a_avg_sd = NaN(size(filt_avg,1), size(lambda.a, 2));
  filt_avg.c_avg_sd = NaN(size(filt_avg,1), size(lambda.c, 2));
  filt_avg.a_avg_n = NaN(size(filt_avg,1), 1);
  filt_avg.c_avg_n = NaN(size(filt_avg,1), 1);
  % add flag for T/S correction (default true: data NOT T/S corrected)
  filt_avg.flag_deltaTSnotcorrected_filt = true(size(filt_avg,1), 1);

  % create filter interpolation table
  filt_interp = table(tot.dt, 'VariableNames', {'dt'});
  filt_interp.a = NaN(size(tot.a));
  filt_interp.c = NaN(size(tot.c));
  filt_interp.flag_deltaTSnotcorrected_total = true(size(tot.dt)); % by default set flag to true until data are corrected
  filt_interp.flag_deltaTSnotcorrected_filt0 = true(size(tot.dt)); % by default set flag to true until data are corrected
  filt_interp.flag_deltaTSnotcorrected_filt1 = true(size(tot.dt)); % by default set flag to true until data are corrected

  % add flag for T/S correction (default true: data NOT T/S corrected)
  filt.flag_deltaTSnotcorrected_filt = true(size(filt.dt));
  % add T/S to filt and filt_interp if TSG loaded
  if tsg_loaded
    % interpolate linear T and S on filt_interp table
    filt_interp = addvars(filt_interp, interp_extrap(tsg_data, filt_interp.dt, TSG.temperature_variable, 30, true, 'linear', 'nearest'), ...
      interp_extrap(tsg_data, filt_interp.dt, TSG.salinity_variable, 30, true, 'linear', 'nearest'), 'NewVariableNames', {'t','s'}, 'After', 'dt');
    % interpolate T and S on filtered data
    filt = addvars(filt, interp_extrap(tsg_data, filt.dt, TSG.temperature_variable, 30, true, 'linear', 'nearest'), ...
      interp_extrap(tsg_data, filt.dt, TSG.salinity_variable, 30, true, 'linear', 'nearest'), 'NewVariableNames', {'t','s'}, 'After', 'dt');
    % add T and S variables on filt_avg table
    filt_avg = addvars(filt_avg, NaN(size(filt_avg,1),1), NaN(size(filt_avg,1),1), 'NewVariableNames', {'t','s'}, 'After', 'dt');
    % remove psiS and psiT from tot and filt with Tref = median(T) of whole processed section
    Tref = median([filt.t; filt_interp.t], 'omitnan');
    a_psiS = interp1(psi.wl, psi.a_psiS, lambda.a, 'spline');
    c_psiS = interp1(psi.wl, psi.c_psiS, lambda.c, 'spline');
    a_psiT = interp1(psi.wl, psi.psiT, lambda.a, 'spline');
    c_psiT = interp1(psi.wl, psi.psiT, lambda.c, 'spline');
    % Correct independently only when t and s are available to avoid getting NaN when t or s are NaN
    %%%%%%%%%%%% TODO add flag when filt or tot are not corrected (t or s not available) %%%%%%%%%%%%
    % correct T/S in filt
    filt_ts_nan = isnan(filt.t) | isnan(filt.s);
    filt.a(~filt_ts_nan, :) = filt.a(~filt_ts_nan, :) - (a_psiT .* (filt.t(~filt_ts_nan) - Tref)) - (a_psiS .* filt.s(~filt_ts_nan));
    filt.c(~filt_ts_nan, :) = filt.c(~filt_ts_nan, :) - (c_psiT .* (filt.t(~filt_ts_nan) - Tref)) - (c_psiS .* filt.s(~filt_ts_nan));
    % correct T/S in tot
    tot_ts_nan = isnan(filt_interp.t) | isnan(filt_interp.s);
    tot.a(~tot_ts_nan, :) = tot.a(~tot_ts_nan, :) - (a_psiT .* (filt_interp.t(~tot_ts_nan) - Tref)) - (a_psiS .* filt_interp.s(~tot_ts_nan));
    tot.c(~tot_ts_nan, :) = tot.c(~tot_ts_nan, :) - (c_psiT .* (filt_interp.t(~tot_ts_nan) - Tref)) - (c_psiS .* filt_interp.s(~tot_ts_nan));
    % adjust T/S correction flags
    filt.flag_deltaTSnotcorrected_filt(~filt_ts_nan) = false; % change flag if corrected
    filt_interp.flag_deltaTSnotcorrected_total(~tot_ts_nan) = false; % change flag if corrected
  end

  % prepare fCDOM data if available
  if ~isempty(cdom_base)
    % interpolate spline and extrapolate nearest fdom on filt_interp table
    filt_interp = addvars(filt_interp, interp_extrap(cdom_base, filt_interp.dt, 'fdom', 30, false, 'linear'), ...
      'NewVariableNames', 'fdom', 'After', 'dt');
    % interpolate fdom on filtered data
    filt = addvars(filt, interp_extrap(cdom_base, filt.dt, 'fdom', 30, false, 'linear'), ...
      'NewVariableNames', 'fdom', 'After', 'dt');
    % add fdom variables on filt_avg table
    filt_avg = addvars(filt_avg, NaN(size(filt_avg,1),1), 'NewVariableNames', 'fdom', 'After', 'dt');
  end

  % compute filter averages
  for i=1:size(sel_start, 1)
    sel_filt = fth_interp.dt(sel_start(i)) <= filt.dt & filt.dt <= fth_interp.dt(sel_end(i));
    foo = filt(sel_filt,:);
    % filter average of a and c
    if sum(sel_filt) == 1
      filt_avg.dt(i) = foo.dt;
      filt_avg.start(i) = foo.dt;
      filt_avg.end(i) = foo.dt;
      filt_avg.a(i,:) = foo.a;
      filt_avg.c(i,:) = foo.c;
      filt_avg.a_avg_sd(i,:) = foo.a_avg_sd;
      filt_avg.c_avg_sd(i,:) = foo.c_avg_sd;
      filt_avg.a_avg_n(i) = foo.a_avg_n;
      filt_avg.c_avg_n(i) = foo.c_avg_n;
      filt_avg.flag_deltaTSnotcorrected_filt(i) = foo.flag_deltaTSnotcorrected_filt;
      % repeat for T/S if tsg is loaded
      if tsg_loaded
        filt_avg.t(i,:) = foo.t;
        filt_avg.s(i,:) = foo.s;
      end
      % repeat for fdom available
      if any(strcmp(foo.Properties.VariableNames, 'fdom'))
        filt_avg.fdom(i,:) = foo.fdom;
      end
    elseif sum(sel_filt) > 1
      a_perc25 = foo.a > prctile(foo.a, 25, 1);
      c_perc25 = foo.c > prctile(foo.c, 25, 1);
      foo.a_avg_sd(a_perc25) = NaN;
      foo.c_avg_sd(c_perc25) = NaN;
      foo.a(a_perc25) = NaN;
      foo.c(c_perc25) = NaN;
      % compute average of all values smaller than 25th percentile for each filter event
      filt_avg.dt(i) = mean(foo.dt(any(~a_perc25, 2) | any(~c_perc25, 2)), 'omitnan');
      if any(any(any(~a_perc25, 2) | any(~c_perc25, 2), 2))
        filt_avg.start(i) = min(foo.dt(any(~a_perc25, 2) | any(~c_perc25, 2)));
        filt_avg.end(i) = max(foo.dt(any(~a_perc25, 2) | any(~c_perc25, 2)));
      else
        filt_avg.start(i) = NaN;
        filt_avg.end(i) = NaN;
      end
      filt_avg.a(i,:) = mean(foo.a, 1, 'omitnan');
      filt_avg.c(i,:) = mean(foo.c, 1, 'omitnan');
      filt_avg.a_avg_sd(i,:) = mean(foo.a_avg_sd, 1, 'omitnan');
      filt_avg.c_avg_sd(i,:) = mean(foo.c_avg_sd, 1, 'omitnan');
      filt_avg.a_avg_n(i) = sum(foo.a_avg_n(any(~isnan(foo.a), 2)), 'omitnan');
      filt_avg.c_avg_n(i) = sum(foo.c_avg_n(any(~isnan(foo.c), 2)), 'omitnan');
      % change flag only if all filtered data are T/S corrected
      if ~any(foo.flag_deltaTSnotcorrected_filt) 
        filt_avg.flag_deltaTSnotcorrected_filt(i) = false;
      end
      % repeat for T/S if tsg is loaded
      if tsg_loaded
        foo.t(foo.t > prctile(foo.t, 25, 1)) = NaN;
        foo.s(foo.s > prctile(foo.s, 25, 1)) = NaN;
        % compute average of all values smaller than 25th percentile for each filter event
        filt_avg.t(i) = mean(foo.t, 1, 'omitnan');
        filt_avg.s(i) = mean(foo.s, 1, 'omitnan');
      end
      % repeat for fdom if available
      if any(strcmp(foo.Properties.VariableNames, 'fdom'))
        foo.fdom(foo.fdom > prctile(foo.fdom, 25, 1)) = NaN;
        % compute average of all values smaller than 25th percentile for each filter event
        filt_avg.fdom(i) = mean(foo.fdom, 1, 'omitnan');
      end
    end
  end
  % remove empty filter events
  rm_filter_event = all(isnan(filt_avg.a), 2) | all(isnan(filt_avg.c), 2);
  filt_avg(rm_filter_event, :) = [];

  % Find switch events from filtered to total in fth_interp table
  total_events = table();
  total_events.start = NaN(size(filt_avg, 1), 1);
  total_events.end = NaN(size(filt_avg, 1), 1);
  for i = 1:size(total_events, 1)
    % find start/end of filter events in filt_interp
    if i == 1 && ~isempty(find(filt_interp.dt <= filt_avg.dt(i), 1, 'last'))
      total_events.start(i) = 1;
      total_events.end(i) = find(filt_interp.dt <= filt_avg.dt(i), 1, 'last');
    elseif i > 1
      total_events.start(i) = find(filt_interp.dt >= filt_avg.dt(i-1), 1, 'first');
      total_events.end(i) = find(filt_interp.dt <= filt_avg.dt(i), 1, 'last');
    end
  end
  if total_events.end(i) < size(filt_interp, 1)
    total_events = [total_events; table(NaN, NaN, 'VariableNames',{'start', 'end'})];
    total_events.start(i+1) = find(filt_interp.dt >= filt_avg.dt(i), 1, 'first');
    total_events.end(i+1) = size(filt_interp, 1);
  end
  total_events(all(isnan([total_events.start total_events.end]), 2), :) = [];

  % adjust T/S correction flags only for remaining filter events
  for i=1:size(total_events, 1) % sel_start
    % adjust flag for T/S correction
    if i == 1 && total_events.start(i) > 1 && ~filt_avg.flag_deltaTSnotcorrected_filt(i)
      filt_interp.flag_deltaTSnotcorrected_filt1(total_events.start(i):total_events.end(i)) = false;
    end
    if i > 1 && ~filt_avg.flag_deltaTSnotcorrected_filt(i-1)
      filt_interp.flag_deltaTSnotcorrected_filt0(total_events.start(i):total_events.end(i)) = false;
    end
    if i < size(total_events, 1) && ~filt_avg.flag_deltaTSnotcorrected_filt(i)
      filt_interp.flag_deltaTSnotcorrected_filt1(total_events.start(i):total_events.end(i)) = false;
    end
    if i == size(total_events, 1) && total_events.end(i) < size(filt_interp, 1)
      filt_interp.flag_deltaTSnotcorrected_filt0(total_events.start(i):total_events.end(i)) = false;
    end
  end
  
  switch interpolation_method
    case 'CDOM'
      % interpolate ag and cg using fcdom
      cluster_param = table();
      cluster_param.time_weight = [0.001; 0.01; 0.1; 1; 5; 10];
      cluster_param.perc_negative_slope = NaN(size(cluster_param.time_weight));
      cluster_param.perc_not_clustered = NaN(size(cluster_param.time_weight));
      cluster_param.filt_interp = cell(size(cluster_param.time_weight));
      cluster_param.filt = cell(size(cluster_param.time_weight));
      cluster_param.regress_stats = cell(size(cluster_param.time_weight));
      % Try multiple weights for time variable used for clustering and keep the one with the least perc_negative_slope,
      % the least perc_not_clustered, and the smallest time_weight
      for i = 1:size(cluster_param,1)
        [cluster_param.filt_interp{i}, cluster_param.filt{i}, cluster_param.regress_stats{i}] = agcg_fdom_interpolation(...
          tot, filt_interp, filt, lambda, filt_avg, param.min_nb_pts_per_cluster, cluster_param.time_weight(i));
        cluster_param.perc_negative_slope(i) = (sum(cluster_param.filt_interp{i}.flag_a_negative_slope(:)) + ...
          sum(cluster_param.filt_interp{i}.flag_c_negative_slope(:))) / (size(filt_interp, 1)*size(lambda.a, 2)*2);
        cluster_param.perc_not_clustered(i) = (sum(cluster_param.filt_interp{i}.flag_a_filt0_not_clustered) + ...
          sum(cluster_param.filt_interp{i}.flag_c_filt0_not_clustered)) / (size(filt_interp, 1)*2);
      end
      cluster_param = sortrows(cluster_param, {'perc_negative_slope','perc_not_clustered','time_weight'});
      fprintf('%.3f time weight selected for clustering\n', cluster_param.time_weight(1))
      filt_interp = cluster_param.filt_interp{1};
      filt = cluster_param.filt{1};
      regression_stats = cluster_param.regress_stats{1};
      
      % id 450 nm wl to plot
      id450_a = abs(lambda.a - 450) == min(abs(lambda.a - 450));
      id450_c = abs(lambda.c - 450) == min(abs(lambda.c - 450));

      % plot clusters
      fprintf('Plotting clustering results ... ')
      figure(33); subplot(2, 2, 1);
      gscatter(datetime(filt.dt, 'ConvertFrom', 'datenum'), filt.a(:, id450_a) ./ filt.fdom, ...
        filt.a_clusters, [], 'vo<s>pdh+*x', 8, 'on', 'Time', 'a/fdom');
      title('a cluster time series'); set(gca, 'Fontsize', 12)
      subplot(2, 2, 2); gsc = gscatter(filt.fdom, filt.a(:,id450_a), filt.a_clusters, [], 'vo<s>pdh+*x');
      xlimit = get(gca, 'XLim');
      ylimit = get(gca, 'YLim');
      for c = 1:size(regression_stats.a.slope, 2)
        rl = refline(regression_stats.a.slope(id450_a, c), regression_stats.a.intercept(id450_a, c));
        set(rl, 'Color', gsc(c).Color, 'LineWidth', 1)
      end
      set(gca, 'XLim', xlimit)
      set(gca, 'YLim', ylimit)
      xlabel('fdom'); ylabel('a'); title('a cluster: a/fdom'); set(gca, 'Fontsize', 12)
      leg = findobj(gcf, 'Type', 'Legend');
      title(leg,'Clusters')
    
      subplot(2, 2, 3);
      gscatter(datetime(filt.dt, 'ConvertFrom', 'datenum'), filt.c(:, id450_c) ./ filt.fdom, ...
        filt.c_clusters, [], 'vo<s>pdh+*x', 8, 'on', 'Time', 'c/fdom');
      title('c cluster time series'); set(gca, 'Fontsize', 12)
      subplot(2, 2, 4); gsc = gscatter(filt.fdom, filt.c(:,id450_c), filt.c_clusters, [], 'vo<s>pdh+*x');
      xlimit = get(gca, 'XLim');
      ylimit = get(gca, 'YLim');
      for c = 1:size(regression_stats.c.slope, 2)
        rl = refline(regression_stats.c.slope(id450_c, c), regression_stats.c.intercept(id450_c, c));
        set(rl, 'Color', gsc(c).Color, 'LineWidth', 1)
      end
      set(gca, 'XLim', xlimit)
      set(gca, 'YLim', ylimit)
      xlabel('fdom'); ylabel('c'); title('c cluster: c/fdom'); set(gca, 'Fontsize', 12)
      leg = findobj(gcf, 'Type', 'Legend');
      leg(1).String = strrep(leg(1).String, 'data', 'regression cluster ');
      leg(3).String = strrep(leg(3).String, 'data', 'regression cluster ');
      title(leg,'Clusters')
      fprintf('done\n')
      drawnow

      % TODO: propagate error from regression and filter event STD
      filt_interp.a_avg_sd = interp1(filt.dt, filt.a_avg_sd, filt_interp.dt, 'linear');
      filt_interp.c_avg_sd = interp1(filt.dt, filt.c_avg_sd, filt_interp.dt, 'linear');
      
      % remove interpolation when there is no data
      filt_interp.a_avg_sd(all(isnan(tot.a),2), :) = NaN;
      filt_interp.a(all(isnan(tot.a),2), :) = NaN;
      filt_interp.c_avg_sd(all(isnan(tot.c),2), :) = NaN;
      filt_interp.c(all(isnan(tot.c),2), :) = NaN;
      
    case 'linear'
      % Interpolate filtered on total linearly
      filt_interp.a = interp1(filt_avg.dt, filt_avg.a, filt_interp.dt, 'linear');
      filt_interp.a = fillmissing(filt_interp.a, 'nearest', 'SamplePoints', filt_interp.dt);
      filt_interp.c = interp1(filt_avg.dt, filt_avg.c, filt_interp.dt, 'linear');
      filt_interp.c = fillmissing(filt_interp.c, 'nearest', 'SamplePoints', filt_interp.dt);
      filt_interp.a_avg_sd = interp1(filt_avg.dt, filt_avg.a_avg_sd, filt_interp.dt, 'linear');
      filt_interp.a_avg_sd = fillmissing(filt_interp.a_avg_sd, 'nearest', 'SamplePoints', filt_interp.dt);
      filt_interp.c_avg_sd = interp1(filt_avg.dt, filt_avg.c_avg_sd, filt_interp.dt, 'linear');
      filt_interp.c_avg_sd = fillmissing(filt_interp.c_avg_sd, 'nearest', 'SamplePoints', filt_interp.dt);
      regression_stats = struct();
    otherwise
    error('Method not supported.');
  end

  % Remove lines full of NaNs or with inf data
  sel2rm = any(~isfinite(tot.a),2) | any(~isfinite(tot.c),2)| all(isnan(tot.a),2) | ...
           all(isnan(tot.c),2) | all(isnan(filt_interp.a),2) | all(isnan(filt_interp.c),2);
  tot(sel2rm,:) = [];
  filt_interp(sel2rm,:) = [];
  
  if exist('visFlag', 'file')
    % id only day to run in all tables to plot
    filt_interp_id = filt_interp.dt >= min(days2run) & filt_interp.dt < max(days2run)+1;
    tot_id = tot.dt >= min(days2run) & tot.dt < max(days2run)+1;
    filt_avg_id = filt_avg.dt >= min(days2run) & filt_avg.dt < max(days2run)+1;
    % plot
    fh = figure(52);
    tiledlayout(2, 1, "TileSpacing", "compact", "Padding", "compact")
    nxt = nexttile;
    visFlag([], filt_interp(filt_interp_id, :), tot(tot_id, :), [], filt_avg(filt_avg_id, :), [], 'a', round(size(tot.a, 2)/2), [], [], 'spd', false, fh);
    if any(strcmp(filt_avg.Properties.VariableNames, 'fdom'))
      yyaxis('left')
      ax1 = gca;
      ax1.YColor = [0	205	205]/255;
      scatter(datetime(filt_interp.dt(filt_interp_id),'ConvertFrom','datenum'), filt_interp.fdom(filt_interp_id), '.', 'MarkerEdgeColor', [0	205	205]/255, 'MarkerEdgeAlpha', 0.5)
      plot(datetime(filt_avg.dt(filt_avg_id), 'ConvertFrom', 'datenum'), filt_avg.fdom(filt_avg_id), 'o', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', [0	205	205]/255)
      ylabel('FDOM (volts)')
      legend('fCDOM (raw)', 'fCDOM (avg filter event)', 'Filtered a interpolated', 'Total a', 'Filtered a (avg filter event)', 'AutoUpdate','off', 'FontSize', 12)
    else
      legend('Filtered a interpolated', 'Total a', 'Filtered a median',...
        'AutoUpdate','off', 'FontSize', 12)
    end
    title('Check filter event interpolation, press q to continue', 'FontSize', 14)
    nxt = nexttile;
    visFlag([], filt_interp(filt_interp_id, :), tot(tot_id, :), [], filt_avg(filt_avg_id, :), [], 'c', round(size(tot.a, 2)/2), [], [], 'spd', false, nxt);
    if any(strcmp(filt_avg.Properties.VariableNames, 'fdom'))
      yyaxis('left')
      ax1 = gca;
      ax1.YColor = [0	205	205]/255;
      scatter(datetime(filt_interp.dt(filt_interp_id),'ConvertFrom','datenum'), filt_interp.fdom(filt_interp_id), '.', 'MarkerEdgeColor', [0	205	205]/255, 'MarkerEdgeAlpha', 0.5)
      plot(datetime(filt_avg.dt(filt_avg_id), 'ConvertFrom', 'datenum'), filt_avg.fdom(filt_avg_id), 'o', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', [0	205	205]/255)
      ylabel('FDOM (volts)')
      legend('fCDOM (raw)', 'fCDOM (avg filter event)', 'Filtered c interpolated', 'Total c', 'Filtered c (avg filter event)', 'AutoUpdate','off', 'FontSize', 12)
    else
      legend('Filtered c interpolated', 'Total c', 'Filtered c median',...
        'AutoUpdate','off', 'FontSize', 12)
    end
    % guiSelectOnTimeSeries(fh);
  end
  
  % figure();
  % subplot(1, 3, 1)
  % scatter(filt_interp.a(filt_interp.flag_linear_interp,40), filt_interp.s(filt_interp.flag_linear_interp), 20, 'k+')
  % scatter(filt_interp.a(~filt_interp.flag_linear_interp,40), filt_interp.s(~filt_interp.flag_linear_interp), 20, 'bo', 'filled')
  % ylabel('Salinity')
  % xlabel('ag')
  % subplot(1, 3, 2)
  % scatter(filt_interp.c(filt_interp.flag_linear_interp,40), filt_interp.s(filt_interp.flag_linear_interp), 20, 'k+')
  % scatter(filt_interp.c(~filt_interp.flag_linear_interp,40), filt_interp.s(~filt_interp.flag_linear_interp), 20, 'bo', 'filled')
  % ylabel('Salinity')
  % xlabel('cg')
  % subplot(1, 3, 3)
  % scatter(filt_interp.fdom(filt_interp.flag_linear_interp), filt_interp.s(filt_interp.flag_linear_interp), 20, 'k+')
  % scatter(filt_interp.fdom(~filt_interp.flag_linear_interp), filt_interp.s(~filt_interp.flag_linear_interp), 20, 'bo', 'filled')
  % ylabel('Salinity')
  % xlabel('fdom')

  % Particulate = Total - FSW
  p = table(tot.dt, 'VariableNames', {'dt'});
  p.ap = tot.a - filt_interp.a;
  p.cp = tot.c - filt_interp.c;
  
  if size(lambda.a, 2) > 50 % perform two separate corrections for ap and cp only for ACS data, not AC9
    % Interpolate wavelengths for Scattering & Residual temperature correction
    ap_for_cpresiduals_corr = interp1(lambda.a', p.ap', lambda.c', 'linear', 'extrap')';
    cp_for_apresiduals_corr = interp1(lambda.c', p.cp', lambda.a', 'linear', 'extrap')';
    % ap Scattering & Residual temperature correction
    % cp Residual correction (for efficiency use the one computed from ap as it should be the same)
    fprintf('ap %s residual temperature and scattering correction ... ', strrep(scattering_correction, '_', ' '))
    if strcmp(scattering_correction, 'Rottgers2013_semiempirical')
      [p.ap, ~] = ResidualTempScatterCorrRottgers_semiempirical(p.ap, cp_for_apresiduals_corr, lambda.a, psi);
      nap_offset = true;
    elseif strcmp(scattering_correction, 'Zaneveld1994_proportional')
      [p.ap, ~] = ResidualTempScatterCorrZaneveld_proportional(p.ap, cp_for_apresiduals_corr, lambda.a, psi);
      nap_offset = false;
    elseif strcmp(scattering_correction, 'ZaneveldRottgers_blended')
      [p.ap, ~] = ResidualTempScatterCorrZaneveldRottgers_blended(p.ap, cp_for_apresiduals_corr, lambda.a, psi);
      nap_offset = true;
    end
    fprintf('Done\n')
    fprintf('cp %s residual temperature and scattering correction ... ', strrep(scattering_correction, '_', ' '))
    if strcmp(scattering_correction, 'Rottgers2013_semiempirical')
      [~, p.cp] = ResidualTempScatterCorrRottgers_semiempirical(ap_for_cpresiduals_corr, p.cp, lambda.c, psi);
      nap_offset = true;
    elseif strcmp(scattering_correction, 'Zaneveld1994_proportional')
      [~, p.cp] = ResidualTempScatterCorrZaneveld_proportional(ap_for_cpresiduals_corr, p.cp, lambda.c, psi);
      nap_offset = false;
    elseif strcmp(scattering_correction, 'ZaneveldRottgers_blended')
      [~, p.cp] = ResidualTempScatterCorrZaneveldRottgers_blended(ap_for_cpresiduals_corr, p.cp, lambda.c, psi);
      nap_offset = true;
    end
    fprintf('Done\n')
  else
    fprintf('ap and cp %s residual temperature and scattering correction ... ', strrep(scattering_correction, '_', ' '))
    if strcmp(scattering_correction, 'Rottgers2013_semiempirical')
      [p.ap, p.cp] = ResidualTempScatterCorrRottgers_semiempirical(p.ap, p.cp, lambda.ref, psi);
      nap_offset = true;
    elseif strcmp(scattering_correction, 'Zaneveld1994_proportional')
      [p.ap, p.cp] = ResidualTempScatterCorrZaneveld_proportional(p.ap, p.cp, lambda.ref, psi);
      nap_offset = false;
    elseif strcmp(scattering_correction, 'ZaneveldRottgers_blended')
      [p.ap, p.cp] = ResidualTempScatterCorrZaneveldRottgers_blended(p.ap, p.cp, lambda.ref, psi);
      nap_offset = true;
    end
  end
  
  % Remove lines full of NaNs (Rottgers2013_semiempirical potentially fail when temperature correct is too large)
  sel2rm = all(isnan(p.ap), 2) | all(isnan(p.cp), 2);
  p(sel2rm, :) = [];
  tot(sel2rm, :) = [];
  filt_interp(sel2rm, :) = [];

  % Propagate error (using geometric mean of measurement errors)
  %   Note: Error is not propagated through Scattering & Residual temperature
  %         correction as required by SeaBASS
  p.ap_sd = sqrt(tot.a_avg_sd.^2 + filt_interp.a_avg_sd.^2);
  p.cp_sd = sqrt(tot.c_avg_sd.^2 + filt_interp.c_avg_sd.^2);
  p.ap_n = tot.a_avg_n;
  p.cp_n = tot.c_avg_n;
  if any(strcmp(filt_interp.Properties.VariableNames, 'fdom'))
    p.fdom = filt_interp.fdom;
  end
  
  % % delete ap spectra full of NaNs
  % p(all(isnan(p.ap),2),:) = [];
  % p(all(isnan(p.cp),2),:) = [];
  
  % Unsmoothing ACS spectra
  % Ron Zaneveld, WET Labs, Inc., 2005
  if size(lambda.a, 2) > 50 % unsmooth only ACS, not AC9
    p = unsmoothACS(p, lambda);
  end
  
  % Auto QC on product spectra
  wla_430 = lambda.a(find(lambda.a <= 430, 1,'last')); % find lower and closest to 430nm wavelength
  wla_700 = lambda.a(find(lambda.a >= 700, 1,'first')); % find higher and closest to 700nm wavelength
  
  % replace values at each end of the spectra when < -0.0015
  p.ap_sd(p.ap < -0.0015 & lambda.a < wla_430) = NaN;
  % p.ap_sd(p.ap < -0.0015 & lambda.a >= wla_700) = NaN;
  p.ap(p.ap < -0.0015 & lambda.a < wla_430) = NaN;
  % p.ap(p.ap < -0.0015 & lambda.a >= wla_700) = NaN;

  % set flag matrix
  flag_varname = {...
    'deltaTS_not_corrected_total','deltaTS_not_corrected_filt0','deltaTS_not_corrected_filt1',...
    'fCDOM_mix_a_cluster0', 'fCDOM_mix_a_cluster1','fCDOM_mix_c_cluster0', 'fCDOM_mix_c_cluster1',...
    'fCDOM_a_cluster_chg', 'fCDOM_c_cluster_chg', 'flag_linear_interp_a', 'flag_linear_interp_c',...
    'a_filt0_not_clustered', 'a_filt1_not_clustered', 'c_filt0_not_clustered', 'c_filt1_not_clustered',...
    'flag_a_negative_slope', 'flag_c_negative_slope',...
    'cp_neg','ap_neg','ap_shape','ap_bubbles','cp_bubbles','ap430_700_neg','cp_over10',...
    'noisy600_650','ap460_640_04_450','positive_ap450_570','poc_flag','chl_ap676lh_flag',...
    'gamma_flag','chl_Halh_flag','HH_mphi_flag','HH_G50_flag','chlratio_flag',...
    'gamma_suspicious','poc_suspicious','chl_ap676lh_suspicious',...
    'chl_Halh_suspicious','HH_G50_mphi_suspicious'};
  flag = array2table(false(size(p,1), length(flag_varname)), 'VariableNames', flag_varname);
  if strcmp(interpolation_method, 'CDOM')
    flag.fCDOM_mix_a_cluster0 = filt_interp.fCDOM_mix_a_cluster0;
    flag.fCDOM_mix_a_cluster1 = filt_interp.fCDOM_mix_a_cluster1;
    flag.fCDOM_mix_c_cluster0 = filt_interp.fCDOM_mix_c_cluster0;
    flag.fCDOM_mix_c_cluster1 = filt_interp.fCDOM_mix_c_cluster1;
    flag.fCDOM_a_cluster_chg = filt_interp.fCDOM_a_cluster_chg;
    flag.fCDOM_c_cluster_chg = filt_interp.fCDOM_c_cluster_chg;
    flag.flag_linear_interp_a = filt_interp.flag_linear_interp_a;
    flag.flag_linear_interp_c = filt_interp.flag_linear_interp_c;
    flag.a_filt0_not_clustered = filt_interp.flag_a_filt0_not_clustered;
    flag.a_filt1_not_clustered = filt_interp.flag_a_filt1_not_clustered;
    flag.c_filt0_not_clustered = filt_interp.flag_c_filt0_not_clustered;
    flag.c_filt1_not_clustered = filt_interp.flag_c_filt1_not_clustered;
    flag.cluster_a_negative_slope = filt_interp.flag_a_negative_slope;
    flag.cluster_c_negative_slope = filt_interp.flag_c_negative_slope;
  else
    flag.flag_linear_interp = true(size(p,1), 1);
  end
  flag.deltaTS_not_corrected_total = filt_interp.flag_deltaTSnotcorrected_total;
  flag.deltaTS_not_corrected_filt0 = filt_interp.flag_deltaTSnotcorrected_filt0;
  flag.deltaTS_not_corrected_filt1 = filt_interp.flag_deltaTSnotcorrected_filt1;

  % delete cp spectra when any cp < -0.0015
  todelete = any(p.cp < -0.0015, 2);
  if sum(todelete) > 0
    fprintf('%.2f%% (%i) spectra failed auto-QC: cp < -0.0015\n', ...
      sum(todelete) / size(p, 1) * 100, sum(todelete))
  end
  flag.cp_neg(todelete) = true;
  bad = [p(todelete, :) table(repmat({'cp < -0.0015'}, ...
    sum(todelete), 1), 'VariableNames', {'QC_failed'})];
  p.cp(todelete, :) = NaN;

  % flag ap spectra when ap430-700 < -0.0015
  toflag = any(p.ap < -0.0015 & lambda.a >= wla_430 & lambda.a <= wla_700, 2);
  if sum(toflag) > 0
    fprintf('%.2f%% (%i) spectra flagged: ap 430-700 < -0.0015\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.ap430_700_neg(toflag) = true;
  % bad = [bad; p(toflag, :) table(repmat({'ap 430-700 < -0.0015'}, ...
  %   sum(toflag), 1), 'VariableNames', {'QC_failed'})];
  % p.ap(toflag, :) = NaN;
  
  % flag ap spectra when any ap < -0.01
  toflag = any(p.ap < -0.01 & lambda.a >= wla_430 & lambda.a <= wla_700, 2);
  if sum(toflag) > 0
    fprintf('%.2f%% (%i) spectra flagged: ap 430-700 < -0.01\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.ap_neg(toflag) = true;
  % bad = [bad; p(todelete, :) table(repmat({'ap 430-700 < -0.01'}, ...
  %   sum(todelete), 1), 'VariableNames', {'QC_failed'})];
  % p.ap(todelete, :) = NaN;
  
  % flag when ap650 > ap676: high NAP
  toflag = any(mean(p.ap(:, lambda.a >= 640 & lambda.a <= 655), 2, 'omitnan') > ...
    mean(p.ap(:, lambda.a >= 670 & lambda.a <= 680), 2, 'omitnan'), 2);
  if sum(toflag) > 0
    fprintf('%.2f%% (%i) spectra flagged: ap650 > ap676 high NAP\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.ap_shape(toflag) = true;
  % bad = [bad; p(toflag, :) table(repmat({'ap650 > ap676: high NAP'}, ...
  %   sum(toflag), 1), 'VariableNames', {'QC_failed'})];
  % p.cp(toflag, :) = NaN;
  
  % flag attenuation spectra when cp > 10
  toflag = any(p.cp > 10, 2);
  if sum(toflag) > 0
    fprintf('%.2f%% (%i) spectra flagged: p.cp > 10\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.cp_over10(toflag) = true;
  % bad = [bad; p(toflag, :) table(repmat({'p.cp > 10'}, ...
  %   sum(toflag), 1), 'VariableNames', {'QC_failed'})];
  % p.cp(toflag, :) = NaN;
  
  % find wavelength below and above which ap and cp are unrealistic and replace by NaNs
  if size(lambda.a, 2) > 50 % clean only ACS data, not AC9
    % replace unrealistic ap and cp red wavelenght by NaN
    [~, foo] = min(p.cp(:, lambda.c > 715),[],2);
    minred_wlc = sum(lambda.c <= 715) + foo;
    [~, foo] = min(abs(p.ap(:, lambda.a > 715)),[],2);
    minred_wla = sum(lambda.a <= 715) + foo;
    p.cp(lambda.c > lambda.c(minred_wlc)' & lambda.c > 720) = NaN;
    p.cp_sd(lambda.c > lambda.c(minred_wlc)' & lambda.c > 720) = NaN;
    p.ap(lambda.a > lambda.a(minred_wla)' & lambda.a > 715) = NaN;
    p.ap_sd(lambda.a > lambda.a(minred_wla)' & lambda.a > 715) = NaN;
    
    % replace unrealistic ap in blue wavelength by NaN
    blue_wl = p.ap;
    blue_wl(:,lambda.a > 550) = NaN;
    blue_wl_var = [zeros(size(blue_wl, 1), 1) abs(diff(blue_wl,[],2))]; % get absolute derivative over wavelengths  
    cutoffblue_wla = NaN(size(blue_wl_var, 1), 1);
    for i = 1:size(cutoffblue_wla,1)
      foo = lambda.a(find(blue_wl_var(i,:) > 6 * mean(blue_wl_var(i, lambda.a > 450), 2, 'omitnan'), 1, 'last'));
      if ~isempty(foo)
        cutoffblue_wla(i) = foo;
      else
        cutoffblue_wla(i) = min(lambda.a);
      end
    end
    if ~all(isnan(cutoffblue_wla))
      p.ap_sd(lambda.a < cutoffblue_wla & lambda.a <= 450) = NaN;
      p.ap(lambda.a < cutoffblue_wla & lambda.a <= 450) = NaN;
    end
    
    % Auto QC from ratio std 600-650 / 640-676
    ap_450 = p.ap(:, find(lambda.a >= 450, 1,'first'));
    ratiod600_ap450 = sum(abs(diff(p.ap(:, lambda.a > 600 & lambda.a <= 650), 2, 2)),2) ./ ap_450;
    
    % get automatic threshold for ratio std 600-650 / 640-676
    fudge_list = (0.1:0.1:10)';
    ndel_spec = NaN(size(fudge_list));
    for i=1:size(fudge_list,1)
      toflag = ratiod600_ap450 > fudge_list(i) * median(ratiod600_ap450, 'omitnan');
      ndel_spec(i) = sum(toflag);
    end
    fudge_factor = fudge_list(find(abs(diff(ndel_spec)) == min(abs(diff(ndel_spec))) & ...
      ndel_spec(2:end) < 0.001 * max(ndel_spec), 1,  'first')); % 0.05 threshold on first derivative of number of spectra deleted
    if isempty(fudge_factor)
      fudge_factor = fudge_list(find(abs(diff(ndel_spec)) == min(abs(diff(ndel_spec))) & ...
        ndel_spec(2:end) < 0.01 * max(ndel_spec), 1,  'first')); % 0.05 threshold on first derivative of number of spectra deleted
    end
    if ~isempty(fudge_factor)
      toflag = ratiod600_ap450 > fudge_factor * median(ratiod600_ap450);
      if sum(toflag) > 0
        fprintf('%.2f%% (%i) spectra flagged: sum(abs(d(ap)/d(lambda(600-650)))) / ap450nm\n', ...
          sum(toflag) / size(p, 1) * 100, sum(toflag))
      end
  
    %   delete bad spectra
      flag.noisy600_650(toflag) = true;
      bad = [bad; p(toflag, :) ...
        table(repmat({'sum(abs(d(ap)/d(lambda(600-650)))) / ap_{450nm}'}, sum(toflag), 1), ...
        'VariableNames', {'QC_failed'})];
  %     p.ap(toflag, :) = NaN;
    end
    % flag potential bubbles using the step in the center of the ap spectra
    foo_difstep = diff(p.ap(:, lambda.a > 560 & lambda.a <= 600),1, 2) ./ diff(lambda.a(lambda.a > 560 & lambda.a <= 600));
    foo_difref = diff(p.ap(:, lambda.a > 500 & lambda.a <= 560),1, 2) ./ diff(lambda.a(lambda.a > 500 & lambda.a <= 560));
    toflag = any(abs(foo_difstep) > 3 * prctile(abs(foo_difref), 95, 2), 2);
    if sum(toflag)
      fprintf('Signal of bubbles/large particles detected in %.2f%% (%i) of the ap spectra, flagged: step in ap between 550 and 600 nm\n', ...
        sum(toflag) / size(p, 1) * 100, sum(toflag))
      flag.ap_bubbles(toflag) = true;
    end
    % flag potential bubbles using the step in the center of the cp spectra
    foo_difstep = diff(p.cp(:, lambda.c > 560 & lambda.c <= 600),1, 2) ./ diff(lambda.c(lambda.c > 560 & lambda.c <= 600));
    foo_difref = diff(p.cp(:, lambda.c > 500 & lambda.c <= 560),1, 2) ./ diff(lambda.c(lambda.c > 500 & lambda.c <= 560));
    toflag = any(abs(foo_difstep) > 3 * prctile(abs(foo_difref), 95, 2), 2);
    if sum(toflag)
      fprintf('Signal of bubbles/large particles detected in %.2f%% (%i) of the cp spectra, flagged: step in cp between 550 and 600 nm\n', ...
        sum(toflag) / size(p, 1) * 100, sum(toflag))
      flag.cp_bubbles(toflag) = true;
    end
  end
  
  % Auto QC when a positive first derivatives of ap over
  % wavelentght between 460 and 640 nm is larger than 0.4 times ap at 450nm
  % or second derivatives of ap over wavelentght between 460 and 640 nm is
  % larger than 0.006
  ap_450 = p.ap(:, find(lambda.a >= 450, 1,'first'));
  d460_640 = diff(p.ap(:, lambda.a > 460 & lambda.a <= 640),[],2);
  % delete bad spectra
  toflag = any(d460_640 > 0.4 * ap_450,2) | any(abs(diff(d460_640,[],2)) > 0.05, 2);
  if sum(toflag)
    fprintf('%.2f%% (%i) spectra flagged: d(ap)/d(lambda460-640) > 0.4 * ap_{450nm} | abs(d"(ap)/d(lambda460-640)) > 0.05)\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.ap460_640_04_450(toflag) = true;
  bad = [bad; p(toflag, :) table(repmat({'d(ap)/d(lambda460-640) > 0.4 * ap_{450nm} | abs(d"(ap)/d(lambda460-640)) > 0.05)'}, ...
    sum(toflag), 1), 'VariableNames', {'QC_failed'})];
  % p.ap(toflag, :) = NaN;
  
  % Auto QC when ap spectra contains 4 consecutive positive first derivatives of ap over
  % wavelentght between 485 and 570 nm
  d485_570 = diff(p.ap(:, lambda.a > 485 & lambda.a <= 570),[],2);
  pos_d485_570 = d485_570 > 0;
  toflag = false(size(pos_d485_570, 1), 1);
  N = 4; % Required number of consecutive numbers following a first one
  for i = 1:size(toflag,1)
    t = [false pos_d485_570(i,:) false];
    if any(find(diff(t)==-1)-find(diff(t)==1)>=N) % First t followed by >=N consecutive numbers
      toflag(i) = true;
    end
  end
  if sum(toflag)
    fprintf('%.2f%% (%i) spectra flagged: 4 consecutive d(ap)/d(lambda485-570) > 0\n', ...
      sum(toflag) / size(p, 1) * 100, sum(toflag))
  end
  flag.positive_ap450_570(toflag) = true;
  bad = [bad; p(toflag, :) table(repmat({'4 consecutive d(ap)/d(lambda485-570) > 0'}, ...
    sum(toflag), 1), 'VariableNames', {'QC_failed'})];
  bad = sortrows(bad, 'dt');
  % p.ap(toflag, :) = NaN;
  
  % % Auto QC when ap spectra contains 5 consecutive positive first derivatives of cp over
  % % wavelentght between 485 and 570 nm
  % d485_570 = diff(p.cp(:, lambda.c > 485 & lambda.c <= 570),[],2);
  % pos_d485_570 = d485_570 > 0;
  % todelete = false(size(pos_d485_570, 1), 1);
  % N = 5; % Required number of consecutive numbers following a first one
  % for i = 1:size(todelete,1)
  %   t = [false pos_d485_570(i,:) false];
  %   if any(find(diff(t)==-1)-find(diff(t)==1)>=N) % First t followed by >=N consecutive numbers
  %     todelete(i) = true;
  %   end
  % end
  % if sum(todelete)
  %   fprintf('%.2f%% (%i) spectra failed auto-QC: 5 consecutive d(cp)/d(lambda485-570) > 0\n', ...
  %     sum(todelete) / size(p, 1) * 100, sum(todelete))
  % end
  % bad = [bad; p(todelete, :) table(repmat({'5 consecutive d(cp)/d(lambda485-570) > 0'}, ...
  %   sum(todelete), 1), 'VariableNames', {'QC_failed'})];
  % bad = sortrows(bad, 'dt');
  % p.cp(todelete, :) = NaN;
  
  % run gaussian decomposition
  agaus = GaussDecomp(p, lambda.a, compute_ad_aphi, nap_offset);
  p = [p agaus];
  
  fprintf('Calculating Chl line height, POC & gamma ... ')
  % Derive standard products from ap and cp
  % Derive POC (Specific to region)
  % 	The particulate organic carbon (POC) is computed using the particulate attenuation at 660 nm Using the global relationship from Gardner et al. (2006)
  % Gardner, W.D., Mishonov, A., Richardson, M.J., 2006. Global POC concentrations from in-situ and satellite data. Deep Sea Res. II 53, 718?740.
  cp660 = interp1(lambda.c, p.cp',660,'linear')';
  p.poc = cp660.*380;
  flag.poc_flag(p.poc < 0) = true;
  % p.poc(p.poc < 0) = NaN;
  
  % Derive Chl (Line heigh at 676 compared to 650 and 715)
  % 	Chlorophyll a (chl) is computed using the particulate absorption line height at 676 nm and the global relationship from Tara Ocean (Boss et al. 2013)
  % REFERENCES:
  % Emmanuel Boss, Marc Picheral, Thomas Leeuw, Alison Chase, Eric Karsenti, Gabriel Gorsky, Lisa Taylor, Wayne Slade, Josephine Ras, Herve Claustre, 2013.The characteristics of particulate absorption, scattering and attenuation coefficients in the surface ocean; Contribution of the Tara Oceans expedition, Methods in Oceanography.
  ap_a = interp1(lambda.a, p.ap',[650 676 715],'linear')';
  p.ap(ap_a(:,1) > ap_a(:,2), :) = NaN; % deleted unrealistic spectra
  ap_a(ap_a(:,1) > ap_a(:,2), :) = NaN; % deleted unrealistic spectra
  p.ap676_lh = ap_a(:,2)-(39/65*ap_a(:,1)+26/65*ap_a(:,3));  % ap_a650-(39/65*ap_a(:,1)+26/65*ap_a(:,3));
  flag.chl_ap676lh_flag(real(p.ap676_lh) ~= p.ap676_lh) = true;
  p.ap676_lh(real(p.ap676_lh) ~= p.ap676_lh) = NaN;
  flag.ap676_lh(p.ap676_lh < 0) = true;
  p.ap676_lh(p.ap676_lh < 0) = NaN;
  p.chl_ap676lh = 157*p.ap676_lh.^1.22;
  
  % 3.3 Derive Gamma (does not support NaN values) (Boss et al. 2001)
  % REFERENCES:
  % Boss, E., W.S. Pegau, W.D. Gardner, J.R.V. Zaneveld, A.H. Barnard., M.S. Twardowski, G.C. Chang, and T.D. Dickey, 2001. Spectral particulate attenuation and particle size distribution in the bottom boundary layer of a continental shelf. Journal of Geophysical Research, 106, 9509-9516.
  % [~,p.gamma] = FitSpectra_HM2(lambda.ref(:,1:end-2),p.cp(:,1:end-2));
  % Correct bu on March 5, 2018, FitSpectra_HM2 does not accept NaNs in cp
  p.gamma = NaN(size(p, 1), 1);
  cp_temp_gam = p.cp(~all(isnan(p.cp),2), :);
  sel = ~any(isnan(cp_temp_gam));
  [~,temp_gam] = FitSpectra_HM2(lambda.c(:,sel), cp_temp_gam(:,sel));
  p.gamma(~all(isnan(p.cp),2)) = temp_gam;
  flag.gamma_flag(p.gamma < 0) = true;
  % p.gamma(p.gamma < 0) = NaN;
  fprintf('Done\n')
  
  fprintf('Calculating chlorophyll from cp (H_alh) ... ')
  % Chlorophyll absorption (alh) and phytoplankton size eigenvectors (P) inferred from cp
  % Houskeeper, H.F., Draper, D., Kudela, R.M., Boss, E., 2020. Chlorophyll absorption and phytoplankton size information inferred from hyperspectral particulate beam attenuation. Appl. Opt. 59, 6765. https://doi.org/10.1364/AO.396832
  % First compute hskpr P parameters (put link to github of hkpr do not include in github).
  [p.Halh, Pr] = houskeeperetal2020(lambda.c, p.cp, false);
  flag.chl_Halh_flag(p.Halh < 0) = true;
  p.Halh(p.Halh < 0) = NaN;
  p.chl_Halh = 157*p.Halh.^1.22;
  fprintf('Done\n')
  
  fprintf('Estimating G50 and mphi (slope of PSD) ... ')
  % Use fit from Haëntjens et al. 2021v22 to get median average
  % cross-sectional area (G50) and slope of phytoplankton size distribution
  % (in abundance) (mphi):
  p.HH_G50 = modelG50.predictFcn(Pr');
  p.HH_mphi = modelmphi.predictFcn(Pr');
  flag.HH_mphi_flag(p.HH_mphi > 0 | p.HH_mphi < -8) = true;
  % p.HH_mphi(p.HH_mphi > 0 | p.HH_mphi < -8) = NaN;
  flag.HH_G50_flag(p.HH_G50 < 0 | p.HH_G50 > 500) = true;
  % p.HH_G50(p.HH_G50 < 0 | p.HH_G50 > 50) = NaN;
  fprintf('Done\n')
  
  % Extra flags for suspicious data
  % gamma_suspicious
  flag.gamma_suspicious(p.gamma < 0.4) = true;
  % poc_suspicious
  flag.poc_suspicious(p.poc > 1000) = true;
  % chl_ap676lh_suspicious
  flag.chl_ap676lh_suspicious(p.chl_ap676lh > 30) = true;
  % chl_Halh_suspicious if ap676_lh is one order of magnitude different from Halh
  flag.chlratio_flag(p.ap676_lh./p.Halh <= 0.1 | p.ap676_lh./p.Halh >= 10 | p.chl_Halh > 30) = true;
  % HH_G50_mphi_suspicious
  flag.HH_G50_mphi_suspicious(p.HH_G50 < 0.3 | p.HH_mphi < -5 | p.HH_G50 > 50) = true;
  
  % set flag column
  p.flag_bit = set_flagbit(flag);
  % flag_info = read_flagbit(p.flag_bit, 'ACS');
  
  %% ag & cg
  if ~isempty(di)
    di_dt = datetime(di.dt, 'ConvertFrom', 'datenum');
    if strcmp(di_method, 'best_di')
      % select DIW with lowest a or c values between 550-650nm
      di_orig = di;
      best_di = table();
      best_di.a = NaN(size(di_orig,1), 1);
      best_di.c = NaN(size(di_orig,1), 1);
      best_di.assign_loc = cell(size(di_orig,1), 1);
      for i = 1:size(di_orig,1)
        if i == 1 || i == size(di_orig,1)
          best_di.assign_loc{i} = abs(di_dt(i) - di_dt) < hours(96); % 72 96 48
        else
          best_di.assign_loc{i} = abs(di_dt(i) - di_dt) < hours(48); % 36 48 24
        end
        % lowest_di_a = di_orig.a(:, lambda.a >= 550 & lambda.a <= 650) == ...
        %   min(di_orig.a(iddi, lambda.a >= 550 & lambda.a <= 650), [], 1);
        % lowest_di_c = di_orig.c(:, lambda.c >= 550 & lambda.c <= 650) == ...
        %   min(di_orig.c(iddi, lambda.c >= 550 & lambda.c <= 650), [], 1);
        % foo_a = find(sum(lowest_di_a, 2) == max(sum(lowest_di_a, 2)));
        % foo_c = find(sum(lowest_di_c, 2) == max(sum(lowest_di_c, 2)));
        best_di.a(i) = find(mean(di_orig.a(:, lambda.a >= 550 & lambda.a <= 650), 2) == ...
          min(mean(di_orig.a(best_di.assign_loc{i}, lambda.a >= 550 & lambda.a <= 650), 2), [], 1), 1, 'first');
        best_di.c(i) = find(mean(di_orig.c(:, lambda.c >= 550 & lambda.c <= 650), 2) == ...
          min(mean(di_orig.c(best_di.assign_loc{i}, lambda.c >= 550 & lambda.c <= 650), 2), [], 1), 1, 'first');
      end
      for i = 1:size(di_orig,1)
        di.a(best_di.assign_loc{i}, :) = repmat(di_orig.a(best_di.a(i), :), sum(best_di.assign_loc{i}), 1);
        di.c(best_di.assign_loc{i}, :) = repmat(di_orig.c(best_di.c(i), :), sum(best_di.assign_loc{i}), 1);
        di.a_avg_sd(best_di.assign_loc{i}, :) = repmat(di_orig.a_avg_sd(best_di.a(i), :), sum(best_di.assign_loc{i}), 1);
        di.c_avg_sd(best_di.assign_loc{i}, :) = repmat(di_orig.c_avg_sd(best_di.c(i), :), sum(best_di.assign_loc{i}), 1);
      end
    
      % visProd3D(lambda.a, di.dt, di.a, false); zlabel('a_diw (m^{-1})');
      % visProd3D(lambda.c, di.dt, di.c, false); zlabel('c_diw (m^{-1})');

    else
      % duplicate DI from the start of each section to the end of each DI section
      di_end = di(1:end-1, :);
      di_end.dt = datenum(di_dt(2:end) - minutes(1));
      di = [di; di_end];
      di = sortrows(di, 'dt');
      if max(di.dt) < max(filt_avg.dt)
        di.dt(end) = datenum(datetime(max(filt_avg.dt), 'ConvertFrom', 'datenum') + minutes(1));
      end
    end
    
% visProd3D(lambda.a, di.dt, di.a, false); zlabel('a_f (m^{-1})');
% visProd3D(lambda.a, di_interp.dt, di_interp.a, false); zlabel('a_f (m^{-1})');

    % remove when a and c are full of NaNs
    filt_avg(all(isnan(filt_avg.a), 2) & all(isnan(filt_avg.c), 2),:) = [];
    
    % Interpolate DI on filtered
    di_interp = table(filt_avg.dt, 'VariableNames', {'dt'});
    di_interp.a = interp1(di.dt, di.a, di_interp.dt, 'linear', 'extrap');
    di_interp.c = interp1(di.dt, di.c, di_interp.dt, 'linear', 'extrap');
    di_interp.a_avg_sd = interp1(di.dt, di.a_avg_sd, di_interp.dt, 'linear', 'extrap');
    di_interp.c_avg_sd = interp1(di.dt, di.c_avg_sd, di_interp.dt, 'linear', 'extrap');
  
    % Dissolved = Filtered - DI
    g = table(filt_avg.dt, 'VariableNames', {'dt'});
    g.ag = filt_avg.a - di_interp.a;
    g.cg = filt_avg.c - di_interp.c;
    % set flag matrix
    flag_g_varname = {'deltaTSnotcorrected_filt'};
    flag_g = array2table(false(size(g,1), length(flag_g_varname)), 'VariableNames', flag_g_varname);
    flag_g.deltaTSnotcorrected_filt = filt_avg.flag_deltaTSnotcorrected_filt;

    % Interpolate wavelength of c on a 
    % g.cg = cell2mat(arrayfun(@(i) interp1(c_wl, g.cg(i,:), a_wl, 'linear', 'extrap'), 1:size(g,1), 'UniformOutput', false)');
    g.ag = interp1(lambda.a', g.ag', lambda.ref', 'linear', 'extrap')';
    g.cg = interp1(lambda.c', g.cg', lambda.ref', 'linear', 'extrap')';
    fprintf('Correcting for temperature & salinity dependence ... ')
    % Temperature & Salinity Correction (No Scattering correction needed)
    [g.ag, g.cg] = TemperatureAndSalinityDependence(g.ag, g.cg, lambda.a, lambda.c, psi, flag_g.deltaTSnotcorrected_filt);
    fprintf('Done\n')
  
    % correct dissolved data for biofouling
    ida450 = abs(lambda.a - 450) == min(abs(lambda.a - 450));
    idc450 = abs(lambda.c - 450) == min(abs(lambda.c - 450));
    [data_corrected, DIW_biofouling_correction] = BiofoulingCorrection(g, lambda, ...
      {'ag', 'cg'}, ida450, flow_data, flow_data.(FTH.view.spd_variable), cdom_base);
    
% visProd3D(lambda.a, data_corrected.dt, data_corrected.ag, false, 'Wavelength', false, 2); zlabel('a_g (m^{-1})');
% visProd3D(lambda.c, data_corrected.dt, data_corrected.cg, false, 'Wavelength', false, 3); zlabel('c_g (m^{-1})');
% visProd3D(lambda.c, data_corrected.dt, data_corrected.ag - data_corrected.cg, false, 'Wavelength', false, 4); zlabel('a_g - c_g (m^{-1})');
% 
% idwl = abs(lambda.a - 450) == min(abs(lambda.a - 450));
% figure; subplot(1,3,1)
% plot_linreg(filt_avg.fdom, data_corrected.ag(:, idwl), 'robust', 'linear', false)
% xlabel('fCDOM (v)')
% ylabel('a_g 450 nm [m^-1]')
% subplot(1,3,2)
% plot_linreg(filt_avg.fdom, data_corrected.cg(:, idwl), 'robust', 'linear', false)
% xlabel('fCDOM (v)')
% ylabel('c_g 450 nm [m^-1]')
% subplot(1,3,3)
% scatter(datetime(filt_avg.dt, 'convertfrom', 'datenum'), data_corrected.ag(:, idwl) - data_corrected.cg(:, idwl), 50, 'filled')
% ylabel('a_g450 - c_g450 [m^-1]')
% set(gca, 'FontSize', 14)

    % force ag to zero at 650 nm and apply this offset to the entire spectra
    forced_NiRto0 = data_corrected;
    ida650 = abs(lambda.a - 650) == min(abs(lambda.a - 650));
    forced_NiRto0.ag = forced_NiRto0.ag - forced_NiRto0.ag(:, ida650);
    % force cg to zero at 650 nm and apply this offset to the entire spectra
    idc650 = abs(lambda.c - 650) == min(abs(lambda.c - 650));
    forced_NiRto0.cg = forced_NiRto0.cg - forced_NiRto0.cg(:, idc650);

% visProd3D(lambda.a, forced_NiRto0.dt, forced_NiRto0.ag, false, 'Wavelength', false, 5); zlabel('a_g forced 0 at 650 nm (m^{-1})');
% visProd3D(lambda.c, forced_NiRto0.dt, forced_NiRto0.cg, false, 'Wavelength', false, 6); zlabel('c_g forced 0 at 650 nm (m^{-1})');
% visProd3D(lambda.a, forced_NiRto0.dt, forced_NiRto0.ag-forced_NiRto0.cg, false, 'Wavelength', false, 7); zlabel('a_g - c_g forced 0 at 650 nm (m^{-1})');
% 
% idwl = abs(lambda.a - 450) == min(abs(lambda.a - 450));
% figure; subplot(1,3,1)
% plot_linreg(filt_avg.fdom, forced_NiRto0.ag(:, idwl), 'robust', 'linear', false)
% xlabel('fCDOM (v)')
% ylabel('a_g450 forced 0 at 650 nm[m^-1]')
% subplot(1,3,2)
% plot_linreg(filt_avg.fdom, forced_NiRto0.cg(:, idwl), 'robust', 'linear', false)
% xlabel('fCDOM (v)')
% ylabel('c_g450 forced 0 at 650 nm [m^-1]')
% subplot(1,3,3)
% scatter(datetime(filt_avg.dt, 'convertfrom', 'datenum'), forced_NiRto0.ag(:, idwl) - forced_NiRto0.cg(:, idwl), 50, 'filled')
% ylabel('a_g450 - c_g450 forced 0 at 650 nm [m^-1]')
% set(gca, 'FontSize', 14)
% 
% figure; hold on
% scatter(datetime(g.dt, 'convertfrom', 'datenum'), g.base_fit_ag, 50, 'filled')
% scatter(datetime(g.dt, 'convertfrom', 'datenum'), g.base_fit_cg, 50, 'filled')
% figure; hold on
% scatter(datetime(g.dt, 'convertfrom', 'datenum'), g.y_intercp_fit_ag, 50, 'filled')
% scatter(datetime(g.dt, 'convertfrom', 'datenum'), g.y_intercp_fit_cg, 50, 'filled')
    
    % if absolute mean difference between ag and cg is lower after forcing use forced data otherwise original data
    if mean(abs(forced_NiRto0.ag(:, ida450) - forced_NiRto0.cg(:, idc450)), 'omitmissing') < ...
        mean(data_corrected.ag(:, ida450) - data_corrected.cg(:, idc450), 'omitmissing')
      g = forced_NiRto0;
      g.ag_650offset = data_corrected.ag(:, ida650);
      g.cg_650offset = data_corrected.ag(:, idc650);
    else
      g = data_corrected;
      g.ag_650offset = zeros(size(g.dt));
      g.cg_650offset = zeros(size(g.dt));
    end

    % Propagate error
    %   Note: Error is not propagated through temprature and salinity dependance
    %         correction as required by SeaBASS
    g.ag_sd = sqrt(filt_avg.a_avg_sd.^2 + di_interp.a_avg_sd.^2);
    g.cg_sd = sqrt(filt_avg.c_avg_sd.^2 + di_interp.c_avg_sd.^2);
    g.ag_n = filt_avg.a_avg_n;
    g.cg_n = filt_avg.c_avg_n;
    g.fdom = filt_avg.fdom;
    
  %   visProd3D(lambda.a, g.dt, g.ag, false, 'Wavelength', false, 70);
  %   title('ag with auto best DI 72h')
  %   saveGraph('ag_with_auto_best_DI_72h', 'fig')
  %   
  %   visProd3D(lambda.c, g.dt, g.cg, false, 'Wavelength', false, 71);
  %   title('cg with auto best DI 72h')
  %   saveGraph('cg_with_auto_best_DI_72h', 'fig')
  
  %   visProd3D(lambda.a, di_interp.dt, di_interp.a, false, 'Wavelength', false, 72);
  %   visProd3D(lambda.c, di_interp.dt, di_interp.c, false, 'Wavelength', false, 73);
  
    fprintf('Exponential fit to ag and cg ... ')
    sel_a = lambda.a >= 410 & lambda.a < 580 & ~any(isnan(g.ag(~all(isnan(g.ag),2), :)));
    sel_c = lambda.c >= 410 & lambda.c < 580 & ~any(isnan(g.cg(~all(isnan(g.cg),2), :)));
    % fit exponential on ag
    [g.y_intercp_fit_ag, g.base_fit_ag, ~, ~, g.RMSE_fit_ag] = FitExp(lambda.a(sel_a), ...
      g.ag(:, sel_a), g.ag_sd(:, sel_a));
    % add fit flag
    g.ag_fitflag = false(size(g, 1), 1);
    g.ag_fitflag(g.RMSE_fit_ag > 0.0025) = true;
    % fit exponential on cg
    [g.y_intercp_fit_cg, g.base_fit_cg, ~, ~, g.RMSE_fit_cg] = FitExp(lambda.c(sel_c), ...
      g.cg(:, sel_c), g.cg_sd(:, sel_c));
    % add fit flag
    g.cg_fitflag = false(size(g, 1), 1);
    g.cg_fitflag(g.RMSE_fit_cg > 0.0025) = true;
    fprintf('Done\n')

    % keep only data 410 < data < 580 nm (blue wl too noisy and temperature impact > 580 nm)
    g.ag(:, ~sel_a) = NaN;
    g.ag_sd(:, ~sel_a) = NaN;
    g.cg(:, ~sel_c) = NaN;
    g.cg_sd(:, ~sel_c) = NaN;

  %   % QC with ag and cg spectra (limited testing on the QC)
  %   g(g.ag(:,1) < 0 & g.cg(:,end-3) < -0.005, :) = [];
  else
    g = table();
  end
end


%% Temperature And Salinity Dependence correction
function [a_corr, c_corr, a_slope, c_slope] = TemperatureAndSalinityDependence(a, c, wl_a, wl_c, psi, salinity_correction)
  % Note that a and c does not have to be on the same wavelength, however the function would need to be edited for that

  % Interpolate Sullivan values on the current ACS
  a_psiT=interp1(psi.wl, psi.psiT, wl_a,'spline'); % PCHIP or SPLINE -> better than linear
  c_psiT=interp1(psi.wl, psi.psiT, wl_c,'spline'); % PCHIP or SPLINE -> better than linear
  a_psiS=interp1(psi.wl, psi.a_psiS, wl_a,'spline');
  c_psiS=interp1(psi.wl, psi.c_psiS, wl_c,'spline');
  % Center psiS on 0 instead of +/- 0.001
  a_psiS = a_psiS - median(a_psiS(wl_a <= 590));
  c_psiS = c_psiS - median(c_psiS(wl_c <= 590));
  
  % Parameters of minization routine
  opts = optimset('fminsearch');      
  opts = optimset(opts,'MaxIter',20000000); 
  opts = optimset(opts,'MaxFunEvals',20000); % 20000
  opts = optimset(opts,'TolX',1e-8);
  opts = optimset(opts,'TolFun',1e-8);
  
  % Wavelength selection
  iwla = wl_a >= 450;
  iwlc = wl_c >= 450;
  % iwl = wl >= 700;
  
  % Init minimization parameters
  deltaT = 10;
  deltaS = 20;
  amp = a(find(~any(isnan(a),2), 1, 'first'),20);
  slope = 0.014;
  
  % Init loop
  n = size(a,1);
  a_deltaT = NaN(n,1);
  a_deltaS = NaN(n,1);
  a_amp = NaN(n,1);
  c_deltaT = NaN(n,1);
  c_deltaS = NaN(n,1);
  c_amp = NaN(n,1);
  a_slope = NaN(n,1);
  c_slope = NaN(n,1);
  
  % Run minimization
  for k=1:n
    if salinity_correction
      % Force Temperature, Salinity, and Slope
      x = fminsearch(@costfun_TSD, [deltaT, deltaS, amp, slope], opts, a(k,iwla), a_psiT(iwla), a_psiS(iwla), wl_a(iwla), salinity_correction);
      a_deltaT(k) = x(1); a_deltaS(k) = x(2); a_amp(k) = x(3); a_slope(k) = x(4);
      x = fminsearch(@costfun_TSD, [deltaT, deltaS, amp, slope], opts, c(k,iwlc), c_psiT(iwlc), c_psiS(iwlc), wl_c(iwlc), salinity_correction);
      c_deltaT(k) = x(1); c_deltaS(k) = x(2); c_amp(k) = x(3); c_slope(k) = x(4);
    else
      % Without Salinity forcing
      x = fminsearch(@costfun_TSD, [deltaT, amp, slope], opts, a(k,iwla), a_psiT(iwla), a_psiS(iwla), wl_a(iwla), salinity_correction);
      a_deltaT(k) = x(1); a_amp(k) = x(2); a_slope(k) = x(3);
      x = fminsearch(@costfun_TSD, [deltaT, amp, slope], opts, c(k,iwlc), c_psiT(iwlc), c_psiS(iwlc), wl_c(iwlc), salinity_correction);
      c_deltaT(k) = x(1); c_amp(k) = x(2); c_slope(k) = x(3);
    end
  end
  
  % Apply correction
  if salinity_correction
    a_corr = a - a_deltaT.*a_psiT - a_deltaS.*a_psiS;
    c_corr = c - c_deltaT.*c_psiT - c_deltaS.*c_psiS;
  else
    a_corr = a - a_deltaT.*a_psiT;
    c_corr = c - c_deltaT.*c_psiT;
  end
  
  % Display the forcing parameters
  % disp([a_deltaT, c_deltaT, a_deltaS, c_deltaS, a_slope, c_slope])
end

% cost function for Temperature And Salinity Dependence correction
function cost = costfun_TSD(x, spectra, psiT, psiS, wl, salinity_correction)
  if salinity_correction
    % Force Temperature, Salinity, and Slope
    cost = sum((spectra - psiT.*x(1) - psiS .* x(2) - x(3) .* exp(-x(4)*(wl-450))).^2);
  else
    % Without Salinity forcing
    cost = sum((spectra - psiT.*x(1) - x(2).*exp(-x(3)*(wl-450))).^2);
  end
end

%% Residual Temperature And Scattering Correction (Zaneveld 1994 proportional)
function [ap_corr, cp_corr] = ResidualTempScatterCorrZaneveld_proportional(ap, cp, wl, psi)
  % Function from Emmanuel Boss improved by Nils Haëntjens
  % Assumes negligible ap in NIR after Zaneveld 1994 method 3
  psiT = interp1(psi.wl, psi.psiT, wl);
  
  % Parameters of minization routine
  opts = optimset('fminsearch');
  % opts = optimset(opts,'NonlEqnAlgorithm', 'gn'); % Does not work on R2017a
  opts = optimset(opts,'MaxIter',20000000);
  opts = optimset(opts,'MaxFunEvals',20000);
  opts = optimset(opts,'TolX',1e-8);
  opts = optimset(opts,'TolFun',1e-8);
  
  % Find Near Infrared & references
  iNIR = 710 <= wl &  wl <= 750;  % spectral srange for optimization (710 to 750nm)
  if isempty(iNIR); error('Unable to perform correction as no wavelength available in NIR.'); end
  % Find nearest wavelength to greater than 715 nm to use as reference for correction
  iref = find(abs(wl - 715) == min(abs(wl - 715)), 1, 'first'); % 715 730
  % If ACS spectra do not go up to 715 nm take the closest wavelength to 715 nm
  if isempty(iref); [~, iref] = max(wl); end % works as there is data in iNIR so lowest wavelength is 710
  
  % Initialize output arrays
  deltaT = NaN(size(ap,1),1);
  
 % Init routine parameter scattering correction
  bp = cp - ap;
  
  % Run minimization routine on good spectra only
  sel = find(all(isfinite(ap),2));
  for k = sel'
    deltaT(k) = fminsearch(@costFun_RTSC, 0, opts, ap(k,:), bp(k,:), psiT, iNIR, iref);         
  end
  ap_corr = ap - psiT.*deltaT - ((ap(:,iref) - psiT(iref).*deltaT) ./ bp(:,iref)) .* bp; 
  cp_corr = cp - psiT.*deltaT;
end


%% Residual Temperature And Scattering Correction (Rottgers2013 semiempirical)
function [ap_corr, cp_corr] = ResidualTempScatterCorrRottgers_semiempirical(ap, cp, wl, psi)
  % Function from Emmanuel Boss after Rottgers et al. 2013
  % Allows absorption in NIR in case of high NAP
  psiT = interp1(psi.wl, psi.psiT, wl);
  
  % Parameters of minization routine
  opts = optimset('fminsearch');
  % opts = optimset(opts,'NonlEqnAlgorithm', 'gn'); % Does not work on R2017a
  opts = optimset(opts,'MaxIter',20000000); 
  opts = optimset(opts,'MaxFunEvals',20000);
  opts = optimset(opts,'TolX',1e-8);
  opts = optimset(opts,'TolFun',1e-8);
  
  % Find Near Infrared & references
  iNIR = 710 <= wl &  wl <= 750;  % spectral srange for optimization (710 to 750nm)
  if isempty(iNIR); error('Unable to perform correction as no wavelength available in NIR.'); end
  % Find nearest wavelength to greater than 715 nm to use as reference for correction
  iref = find(abs(wl - 715) == min(abs(wl - 715)), 1, 'first'); % 715 730
  % If ACS spectra do not go up to 715 nm take the closest wavelength to 715 nm
  if isempty(iref); [~, iref] = max(wl); end % works as there is data in iNIR so lowest wavelength is 710
  
  % Initialize output arrays
  deltaT = NaN(size(ap,1),1);
  
  % Init routine parameter scattering correction
  bp = cp - ap;
  
  % Run minimization routine on good spectra only
  sel = find(all(isfinite(ap),2));
  for k = sel'
    deltaT(k) = fminsearch(@costFun_RTSC, 0, opts, ap(k,:), bp(k,:), psiT, iNIR, iref);         
  end
  % apply temperature correction and replace negative values at iref by NaN to avoid complex solution leakage
  ap_Tcorr = ap - psiT.*deltaT;
  cp_corr = cp - psiT.*deltaT;
  ap_Tcorr_iref = ap_Tcorr(:,iref);
  ap_Tcorr_iref(ap_Tcorr_iref < 0) = NaN;
  % apply flat scattering correction
  ap_corr = ap_Tcorr - (ap_Tcorr_iref - 0.212*ap_Tcorr_iref.^1.135);
end

%% Residual Temperature And Scattering Correction (Zaneveld 1994 proportional + NIR offset from Rottgers2013 semiempirical)
function [ap_corr, cp_corr] = ResidualTempScatterCorrZaneveldRottgers_blended(ap, cp, wl, psi)
  % Function from Emmanuel Boss after Zaneveld 1994 method 3 and Rottgers et al. 2013
  % Scattering correction + allows absorption in NIR in case of high NAP
  psiT = interp1(psi.wl, psi.psiT, wl);
  
  % Parameters of minization routine
  opts = optimset('fminsearch');
  % opts = optimset(opts,'NonlEqnAlgorithm', 'gn'); % Does not work on R2017a
  opts = optimset(opts,'MaxIter',20000000); 
  opts = optimset(opts,'MaxFunEvals',20000);
  opts = optimset(opts,'TolX',1e-8);
  opts = optimset(opts,'TolFun',1e-8);
  
  % Find Near Infrared & references
  iNIR = 710 <= wl &  wl <= 750;  % spectral srange for optimization (710 to 750nm)
  if isempty(iNIR); error('Unable to perform correction as no wavelength available in NIR.'); end
  % Find nearest wavelength to greater than 730 nm to use as reference for correction
  iref = find(abs(wl - 715) == min(abs(wl - 715)), 1, 'first'); % 715 730
  % If ACS spectra do not go up to 715 nm take the closest wavelength to 715 nm
  if isempty(iref); [~, iref] = max(wl); end % works as there is data in iNIR so lowest wavelength is 710
  
  % Initialize output arrays
  deltaT = NaN(size(ap,1),1);
  
  % Init routine parameter scattering correction
  bp = cp - ap;
  
  % Run minimization routine on good spectra only
  sel = find(all(isfinite(ap),2));
  for k = sel'
    % deltaT(k) = fminsearch(@costFun_RTSC_RottgersBoss, 0, opts, ap(k,:), bp(k,:), psiT, iNIR, iref);         
    deltaT(k) = fminsearch(@costFun_RTSC, 0, opts, ap(k,:), bp(k,:), psiT, iNIR, iref);         
  end
  % apply temperature correction and replace negative values at iref by NaN to avoid complex solution leakage
  ap_Tcorr = ap - psiT.*deltaT;
  cp_corr = cp - psiT.*deltaT;
  ap_Tcorr_iref = ap_Tcorr(:,iref);
  ap_Tcorr_iref(ap_Tcorr_iref < 0) = NaN;
  bp_Tcorr = cp_corr - ap_Tcorr;
  % apply scattering correction
  ap_corr = ap_Tcorr - ap_Tcorr_iref ./ bp_Tcorr(:,iref) .* bp_Tcorr + 0.212*ap_Tcorr_iref.^1.135;
end

%% Residual Temperature And Scattering Correction cost function 
function cost = costFun_RTSC(deltaT, ap, bp, psiT, iNIR, iref)
  cost = sum(abs(ap(iNIR) - psiT(iNIR).*deltaT - ((ap(iref)-psiT(iref).*deltaT)./bp(iref)).*bp(iNIR)));
end

%% Temperature And Salinity Correction
function [a_ts, c_ts] = TemperatureAndSalinityCorrection(a, c, a_wl, c_wl, delta_t, delta_s, psi)

  %interpolate literature psiT values to acs wavelengths
  a_psi_t = interp1(psi.wl, psi.psiT, a_wl, 'linear', 'extrap');
  c_psi_t = interp1(psi.wl, psi.psiT, c_wl, 'linear', 'extrap');
  a_psi_s = interp1(psi.wl, psi.a_psiS, a_wl, 'linear', 'extrap');
  c_psi_s = interp1(psi.wl, psi.c_psiS, c_wl, 'linear', 'extrap');
  %a_sigma_psi_t = interp1(psi.wl, psi.a_sigma_psiS, wl_a, 'linear', 'extrap');
  %c_sigma_psi_t = interp1(psi.wl, psi.c_sigma_psiS, wl_c, 'linear', 'extrap');
    
  %correct acdata_raw for temp-dependent water absorbance, propagate error
  %into del_raw.  Output acdata_t and del_t.
   
  a_ts = a - (a_psi_t(ones(size(delta_t,1),1),:) .* delta_t(:,ones(size(a_psi_t,2),1)))...
    + (a_psi_s(ones(size(delta_s,1),1),:) .* delta_s(:,ones(size(a_psi_s,2),1)));
  c_ts = c - (c_psi_t(ones(size(delta_t,1),1),:) .* delta_t(:,ones(size(c_psi_t,2),1)))...
    + (c_psi_s(ones(size(delta_s,1),1),:) .* delta_s(:,ones(size(c_psi_s,2),1)));
  
  %c_ts = c - (c_psi_t .* delta_t + c_psi_s .* delta_s);
  %a_del_t(:,1) = ((del_rawa(:,1)).^2 + (a_sigma_psi_t(:,1).*delta_t).^2).^(1/2);
  %c_del_t(:,1) = ((del_rawc(:,1)).^2 + (c_sigma_psi_t(:,1).*delta_t).^2).^(1/2);
end


%%
function a_corr = ScatteringCorrection(a_wl, c_wl, a, c, method)
  % Interpolate wavelength of c on a
  c = interp1(c_wl', c, a_wl', 'linear', 'extrap');
  
  switch method
    case 'flat'
      % spectrally flat correction
      a_730 = interp1(wl, a, 730, 'linear', 'extrap');
      a_corr = a - a_730;
    case 'varying'
      % spectrally varying scattering correction
      b = c - a;
      a_730 = interp1(wl, a, 730, 'linear', 'extrap');
      b_730 = interp1(wl, b, 730, 'linear', 'extrap');
      a_corr = a - b .* a_730 ./ b_730;
    case 'rottgers'
      % Rottgers et al., 2013
      a_715 = interp1(wl, a, 715, 'linear', 'extrap');
      c_715 = interp1(wl, c, 715, 'linear', 'extrap');
      a_corr = a - a_715 .* (1/0.56 .* c - a)./(1/0.56 .* c_715 - a_715);
    otherwise
      error('Method not supported');
  end
end


%%
function acs_unsmoothed = unsmoothACS(acs_data, lambda)
  % AC-S "un-smoothing" and spectral decomposition method
  % Ron Zaneveld, WET Labs, Inc., 2005
  % Ali Chase, University of Maine, 2014
  %
  % Unsmoothing method from spectral decomposition in:
  % Chase, A., et al., Decomposition of in situ particulate absorption
  % spectra. Methods in Oceanography (2014), http://dx.doi.org/10.1016/j.mio.2014.02.022
  %%
  todo = acs_data.Properties.VariableNames(contains(acs_data.Properties.VariableNames, ...
    {'ap', 'cp'}) & ~contains(acs_data.Properties.VariableNames, {'_sd','_n'}));
  
  acs_unsmoothed = table();
  acs_unsmoothed.dt = acs_data.dt;
  for i = todo
    fprintf([i{:} ' unsmoothing ... '])
    % Set up filter factors at every 0.1 nm from 1 to 799 nm, with center
    % wavelength at centwavel (i.e. at the data wavelengths)
    wavelength = .1:.1:799; % Thus index of 1 nm = 10; 356 nm= 3560;
    SIG1 = (-9.845*10^-8.*lambda.(i{:}(1)).^3 + 1.639*10^-4*lambda.(i{:}(1)).^2 - 7.849*10^-2*lambda.(i{:}(1)) + 25.24)/2.3547 ;
    for j = 1:size(lambda.(i{:}(1)),2)
      for jkl = 1:size(wavelength,2)
        filtfunc(jkl,j) = (1/(sqrt(2*pi)*SIG1(j)))*exp(-0.5*((wavelength(jkl)-lambda.(i{:}(1))(j))/SIG1(j)).^2); % First term normalizes area under the curve to 1.
      end
    end
  
    % Convolve the measurement with the fiter factors add the difference to
    % the measured spectrum to get the first corrected spectrum.
    % This is the corrected absorption spectrum "ap".
    minwavel = min(lambda.(i{:}(1)));
    maxwavel = max(lambda.(i{:}(1)));
  
    centwavel = minwavel:.1:maxwavel;% The range of centwavel is 0.1 nm.
    splinap = spline(lambda.(i{:}(1)), acs_data.(i{:}), centwavel); % Spline the measured data to every 0.1 nm.
    % We need data from 0 to 799 nm to multiply by filtfac.
    absspec = zeros(size(acs_data.(i{:}),1), size(wavelength,2));
    absspec(:, minwavel*10:maxwavel*10) = splinap;
    absspec(:, 1:minwavel*10-1) = ones(1, size(1:minwavel*10-1,2)) .* absspec(:, minwavel*10);
    aspecprime = absspec';
  
    meassignal6 = NaN(size(aspecprime, 2), size(lambda.(i{:}(1)), 2));
    parfor j = 1:size(aspecprime, 2)        
      measur2 = aspecprime(:,j) .* filtfunc; % the measured signal for every filter factor.
      meassignal6(j,:) = 0.1 * sum(measur2); % The measured spectrum at a wavelength i is the sum of what a filter measured at
    end
    acs_unsmoothed.(i{:}) = acs_data.(i{:}) - meassignal6 + acs_data.(i{:});
    fprintf('Done\n')
  end
  
  for i = 1:size(acs_data,2)
    if ~any(strcmp(acs_unsmoothed.Properties.VariableNames, acs_data.Properties.VariableNames{i}))
      acs_unsmoothed = [acs_unsmoothed acs_data(:,i)];
    end
  end
end


%%
function agaus = GaussDecomp(p, lambda, compute_ad_aphi, nap_offset)
  % Gaussian decomposition
  % Ron Zaneveld, WET Labs, Inc., 2005
  % Ali Chase, University of Maine, 2014
  % Adaptation to InLineAnalysis: Guillaume Bourdin, March 2021
  %
  % Reference:
  % Chase, A., et al., Decomposition of in situ particulate absorption
  % spectra. Methods in Oceanography (2014), http://dx.doi.org/10.1016/j.mio.2014.02.022
  
  %% identify lines full of NaNs to reconstruct table of same size
  idnan = all(isnan(p.ap),2);
  % delete row full of NaN
  p(idnan, :) = [];
  
  % delete row and columns full of nans
  lambda(all(isnan(p.ap),1)) = [];
  p.ap_sd(:, all(isnan(p.ap),1)) = [];
  p.ap(:, all(isnan(p.ap),1)) = [];
  if nap_offset
    ap_715 = interp1(lambda, p.ap', 715, 'linear')';
    p.ap = p.ap - ap_715;
  end
  
  % Peak center values ("peak_loc") determined using a interative approach that allows the location to vary (uses the matlab
  % function LSQNONLIN), and are rounded to nearest integer. Sigma values ("lsqsig") are determined similarly. FWHM = sigma*2.355
  peak_loc = [406,434,453,470,492,523,550,584,617,638,660,675];
  lsqsig = [16,12,12,13,16,14,14,16,13,11,11,10];
  onenm = 400:1:720;
  
  fprintf('Gaussian decomposition ... ')
  
  % extrapolate NaN tails with constant nearest non-NaN value
  ap_filled = fillmissing(p.ap, 'nearest', 2, 'SamplePoints', lambda);
  ap_sd_filled = fillmissing(p.ap_sd, 'nearest', 2, 'SamplePoints', lambda);
  
  % interpolate the un-smoothed ap spectra to one nm resolution
  acorr2onenm = interp1(lambda, ap_filled', onenm, 'spline');
  
  % define the matrix of component Gaussian functions using the peaks
  % and widths (sigma) above
  coef2 = exp(-0.5 .* (((onenm .* ones(size(peak_loc,2),1))' - peak_loc .* ...
    ones(size(onenm,2),1)) ./ lsqsig) .^ 2);
  
  % define a function for non-algal particles and concatenate this to the Gaussian matrix
  coef2nap = exp(-0.01 * (onenm - 400));
  coef2 = [coef2nap', coef2];
  
  % normalize both the component functions and the measured ap
  % spectra by the uncertainty (std dev) in the ap spectra
  ap_sd_int = interp1(lambda, ap_sd_filled', onenm, 'linear', 'extrap');
  acorr2onenm_new = (acorr2onenm ./ ap_sd_int);
  
  amps = NaN(size(p,1), size(coef2,2));
  sumspec_temp = NaN(size(coef2,1), size(p,1));
  % compspec_temp = NaN(size(coef2,1), size(coef2,2), size(p,1));
  
  parfor i = 1:size(acorr2onenm_new,2)
    coef2_new = coef2 ./ ap_sd_int(:,i);
    % Inversion analysis
    amps(i, :) = lsqnonneg(coef2_new, acorr2onenm_new(:, i));
    % Using the inverted amplitudes, build a matrix of new component
    % spectra (compspec) and the sum of the Gaussian and nap functions (sumspec):
    sumspec_temp(:, i) = sum(amps(i, :) .* coef2, 2);
  %   compspec_temp(:, :, i) = amps(i, :) .* coef2;
  end
  
  % interpolate back to the original resolution
  % compspec = interp1(onenm', compspec_temp, lambda, 'spline');
  sumspec = interp1(onenm, sumspec_temp, lambda, 'spline')';
  
  uncertainty = sum(abs(p.ap(:, lambda > 440 & lambda < 705) - ...
    sumspec(:, lambda > 440 & lambda < 705)), 2, 'omitnan') / ...
      size(lambda(lambda > 440 & lambda < 705),2);
    
  % reconstruct matrix with NaN in the same place
  agaus = NaN(size(idnan,1), 14);
  agaus(~idnan,:) = [amps uncertainty];
  agaus = array2table(agaus, 'VariableNames', ...
      [{'ad_model400'} cellfun(@(x) ['agaus' x], cellstr(num2str(peak_loc'))', 'un', 0) {'agaus-mae'}]);
  
  if compute_ad_aphi
    % Compute a non-algal and a phytoplankton from Zheng, G., and D. Stramski (2013) model
    % Reference: 
    % - Zheng, G., and D. Stramski (2013), A model based on stacked-constraints approach for partitioning the light absorption coefficient of seawater 
    % into phytoplankton and non-phytoplankton components, J. Geophys. Res. Oceans, 118, 2155?2174, doi:10.1002/jgrc.20115.
    % - Zheng, G., and D. Stramski (2013), A model for partitioning the light absorption coefficient of suspended marine particles into phytoplankton 
    % and nonalgal components, J. Geophys. Res. Oceans, 118, 2977?2991, doi:10.1002/jgrc.20206.
    wl_ZS13 = [400 412 420 430 443 450 467 490 500 510 550 555 630 650 670 700];
    qwl = unique([wl_ZS13(:)' 442 676]);
    ap_ZS13 = interp1(lambda, ap_filled', qwl, 'linear', 'extrap'); % need to extrapolate for 400
  
  %   % vectorized partition_ap (DOESN'T WORK YET)
  %   agaus.ad_ZS13 = NaN(size(agaus, 1), size(qwl, 2));
  %   agaus.aphi_ZS13 = NaN(size(agaus, 1), size(qwl, 2));
  %   [agaus.ad_ZS13, agaus.aphi_ZS13] = partition_ap_vec(ap_ZS13, qwl', 50);  % Must be in colum direction
    
    % Run partition_ap
    ad_ZS13 = NaN(size(ap_ZS13, 2), size(qwl, 2));
    aphi_ZS13 = NaN(size(ap_ZS13, 2), size(qwl, 2));
  %   parfor i = 1:size(ap_ZS13, 2)
    for i = progress(1:size(ap_ZS13, 2))
      [ad_ZS13(i,:), aphi_ZS13(i,:)] = partition_ap(ap_ZS13(:, i), qwl', 50);
    end
  %   figure; hold on
  %   plot(qwl, ad)
  %   plot(qwl, aphi_ZS13)
  
    agaus.ad_ZS13 = NaN(size(agaus, 1), size(qwl, 2));
    agaus.aphi_ZS13 = NaN(size(agaus, 1), size(qwl, 2));
    agaus.ad_ZS13(~idnan,:) = ad_ZS13;
    agaus.aphi_ZS13(~idnan,:) = aphi_ZS13;
    agaus.Properties.VariableUnits = repmat({'1/m'}, 1, size(agaus, 2));
  end
  fprintf('Done\n')
end

function [filt_interp, filt, regress_stats] = agcg_fdom_interpolation(tot, filt_interp, filt, lambda, ...
  filt_avg, min_nb_pts_per_cluster, time_weight_for_cluster)
  % Reconstruct ag and cg from fCDOM data
  [filt_avg, filt_interp, filt, regress_stats] = fdom_agcg_model(filt, filt_interp, lambda, filt_avg, ...
    min_nb_pts_per_cluster, time_weight_for_cluster);

  n_periods = size(filt_avg,1)-1;

  % visProd3D(lambda.a, filt_interp.dt, filt_interp.intercept_interp_a, false, 'Wavelength', false, 27);
  % title('intercept a')
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.intercept_interp_c, false, 'Wavelength', false, 28);
  % title('intercept c')
  % 
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.slope_interp_a, false, 'Wavelength', false, 27);
  % title('slope a')
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.slope_interp_c, false, 'Wavelength', false, 28);
  % title('slope c')

  % % interpolate based on fCDOM: Emmanuel's method (local)
  % % allocate fdom0, a0, and c0 varialbes
  % filt_interp.fdom0 = NaN(size(filt_interp.fdom));
  % filt_interp.a0 = NaN(size(filt_interp.a));
  % filt_interp.c0 = NaN(size(filt_interp.c));
  % 
  % % interpolate based on fCDOM: Emmanuel's method
  % delta_acf = table(NaN(n_periods, 1), 'VariableNames', {'dt'});
  % delta_acf.delta_af = NaN(n_periods, size(filt_avg.a, 2));
  % delta_acf.delta_cf = NaN(n_periods, size(filt_avg.a, 2));
  % % Prepare coefficients for each total period going from t0 to t1, starting and finishing by a filtered time
  % for i=1:n_periods
  %   it0 = i; it1 = i + 1;
  %   it = filt_avg.dt(it0) <= filt_interp.dt & filt_interp.dt <= filt_avg.dt(it1);
  %   if any(it)
  %     % interpolation based on fCDOM if enough dynamic range in fCDOM
  %     if any((filt_interp.fdom(it) - cdom_dark) / cdom_dark > 0.05)
  %       % prepare delta a / delta f and delta c / delta f
  %       delta_acf.dt(i) = mean([filt_avg.dt(it0) filt_avg.dt(it1)], 'omitnan');
  %       delta_acf.delta_af(i, :) = (filt_avg.a(it1,:) - filt_avg.a(it0,:)) / (filt_avg.fdom(it1) - filt_avg.fdom(it0));
  %       delta_acf.delta_cf(i, :) = (filt_avg.c(it1,:) - filt_avg.c(it0,:)) / (filt_avg.fdom(it1) - filt_avg.fdom(it0));
  %       filt_interp.fdom0(it) = repmat(filt_avg.fdom(it0), sum(it), 1);
  %       filt_interp.a0(it,:) = repmat(filt_avg.a(it0,:), sum(it), 1);
  %       filt_interp.c0(it,:) = repmat(filt_avg.c(it0,:), sum(it), 1);
  %       filt_interp.flag_linear_interp(it) = false;
  %     else
  %       filt_interp.flag_linear_interp(it) = true;
  %     end
  %   end
  % end
  % delta_acf(all(isnan(delta_acf.delta_af), 2) | all(isnan(delta_acf.delta_cf), 2), :) = [];
  % % interpolated delta a / delta fdom over total events
  % filt_interp.delta_af = interp_extrap(delta_acf, filt_interp.dt, 'delta_af', [], true, 'linear', 'nearest');
  % filt_interp.delta_cf = interp_extrap(delta_acf, filt_interp.dt, 'delta_cf', [], true, 'linear', 'nearest');
  % % extrapolate a0, c0, fdom0
  % filt_interp.fdom0 = fillmissing(filt_interp.fdom0, 'nearest', 'SamplePoints', filt_interp.dt);
  % filt_interp.a0 = fillmissing(filt_interp.a0, 'nearest', 'SamplePoints', filt_interp.dt);
  % filt_interp.c0 = fillmissing(filt_interp.c0, 'nearest', 'SamplePoints', filt_interp.dt);
  % % compute ag and cg
  % filt_interp.a = filt_interp.al + filt_interp.delta_af .* (filt_interp.fdom - filt_interp.fdoml);
  % filt_interp.c = filt_interp.cl + filt_interp.delta_cf .* (filt_interp.fdom - filt_interp.fdoml);
  % % linearly interpolate ag and cg
  % filt_interp.lin_a = interp1(filt_avg.dt, filt_avg.a, filt_interp.dt, 'linear');
  % filt_interp.lin_c = interp1(filt_avg.dt, filt_avg.c, filt_interp.dt, 'linear');
  % % replace fdom interpolation by linear for filter events where fdom is constant
  % filt_interp.a(filt_interp.flag_linear_interp) = filt_interp.lin_a(filt_interp.flag_linear_interp);
  % filt_interp.c(filt_interp.flag_linear_interp) = filt_interp.lin_c(filt_interp.flag_linear_interp);

  % visProd3D(lambda.a, filt_interp.dt, filt_interp.a, false, 'Wavelength', false, 27);
  % title('ag fdom interpolated')
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.lin_a, false, 'Wavelength', false, 28);
  % title('ag linearly interpolated')
  % 
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.c, false, 'Wavelength', false, 29);
  % title('cg fdom interpolated')
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.lin_c, false, 'Wavelength', false, 30);
  % title('cg linearly interpolated')

  % interpolate ag and cg based on fdom: Guillaume's method (global)
  % allocate new variables
  filt_interp.fCDOM_mix_a_cluster0 = false(size(filt_interp.dt));
  filt_interp.fCDOM_mix_a_cluster1 = false(size(filt_interp.dt));
  filt_interp.fCDOM_mix_c_cluster0 = false(size(filt_interp.dt));
  filt_interp.fCDOM_mix_c_cluster1 = false(size(filt_interp.dt));
  filt_interp.fCDOM_a_cluster_chg = false(size(filt_interp.dt));
  filt_interp.fCDOM_c_cluster_chg = false(size(filt_interp.dt));
  filt_interp.flag_linear_interp_a = false(size(filt_interp.dt));
  filt_interp.flag_linear_interp_c = false(size(filt_interp.dt));
  filt_interp.flag_a_filt0_not_clustered = false(size(filt_interp.dt));
  filt_interp.flag_a_filt1_not_clustered = false(size(filt_interp.dt));
  filt_interp.flag_c_filt0_not_clustered = false(size(filt_interp.dt));
  filt_interp.flag_c_filt1_not_clustered = false(size(filt_interp.dt));

  %%%%% RE-WRITTEN ON 2023-12-29 and RE-WRITTEN AGAIN ON 2024-11-23
  % linear interpolation of a, c, and fdom
  filt_interp.lin_a = interp1(filt_avg.dt, filt_avg.a, filt_interp.dt, 'linear');
  filt_interp.lin_c = interp1(filt_avg.dt, filt_avg.c, filt_interp.dt, 'linear');
  filt_interp.lin_fdom = interp1(filt_avg.dt, filt_avg.fdom, filt_interp.dt, 'linear');
  filt_interp.lin_a = fillmissing(filt_interp.lin_a, 'nearest', 'SamplePoints', filt_interp.dt);
  filt_interp.lin_c = fillmissing(filt_interp.lin_c, 'nearest', 'SamplePoints', filt_interp.dt);
  filt_interp.lin_fdom = fillmissing(filt_interp.lin_fdom, 'nearest', 'SamplePoints', filt_interp.dt);

  for i=1:n_periods
    it0 = i; it1 = i + 1;
    it = filt_avg.dt(it0) <= filt_interp.dt & filt_interp.dt <= filt_avg.dt(it1);
    if any(it)
      % set flags for fCDOM interpolation
      if any(~isnan(filt_interp.fdom(it)))
        % flag when mixed a clusters within filter event 0 or filter event 1
        idnonnan = any(~isnan(tot.a), 2);
        if any(it & idnonnan)
          if filt_avg.mix_cluster_a_flag(it0)
            % fprintf('a filter interpolation #%i flagged: mixed cluster within start filter event\n', i)
            filt_interp.fCDOM_mix_a_cluster0(it & idnonnan) = true;
          end
          if filt_avg.mix_cluster_a_flag(it1)
            % fprintf('a filter interpolation #%i flagged: mixed cluster within end filter event\n', i)
            filt_interp.fCDOM_mix_a_cluster1(it & idnonnan) = true;
          end
          % flag when a clusters change between filter event 0 and filter event 1
          if any(~isnan(filt_avg.cluster_weighted_a(it0,:))) && any(~isnan(filt_avg.cluster_weighted_a(it1,:))) && ...
              any(filt_avg.cluster_weighted_a(it0,:) ~= filt_avg.cluster_weighted_a(it1, :))
            % fprintf('a filter interpolation #%i flagged: change of clusters between start and end filter events\n', i)
            filt_interp.fCDOM_a_cluster_chg(it & any(idnonnan)) = true;
          end
          % flag when previous filter event was not clustered
          if filt_avg.a_not_clustered_flag(it0)
            filt_interp.flag_a_filt0_not_clustered(it & idnonnan) = filt_avg.a_not_clustered_flag(it0);
          end
          % flag when next filter event was not clustered
          if filt_avg.a_not_clustered_flag(it1)
            filt_interp.flag_a_filt1_not_clustered(it & idnonnan) = filt_avg.a_not_clustered_flag(it1);
          end
          % set linear interpolation boolean for missing fdom data
          filt_interp.flag_linear_interp_a(it & isnan(filt_interp.fdom)) = true;
        end
        % flag when mixed c clusters within filter event 0 or filter event 1
        idnonnan = any(~isnan(tot.c), 2);
        if any(it & idnonnan)
          if filt_avg.mix_cluster_c_flag(it0)
            % fprintf('c filter interpolation #%i flagged: mixed cluster within start filter event\n', i)
            filt_interp.fCDOM_mix_c_cluster0(it & idnonnan) = true;
          end
          if filt_avg.mix_cluster_c_flag(it1)
            % fprintf('c filter interpolation #%i flagged: mixed cluster within end filter event\n', i)
            filt_interp.fCDOM_mix_c_cluster1(it & idnonnan) = true;
          end
          % flag when c clusters change between filter event 0 and filter event 1
          if ~all(isnan(filt_avg.cluster_weighted_c(it0,:))) && ~all(isnan(filt_avg.cluster_weighted_c(it1,:))) && ...
              any(filt_avg.cluster_weighted_c(it0,:) ~= filt_avg.cluster_weighted_c(it1, :))
            % fprintf('c filter interpolation #%i flagged: change of clusters between start and end filter events\n', i)
            filt_interp.fCDOM_c_cluster_chg(it & any(idnonnan)) = true;
          end
          % flag when previous filter event was not clustered
          if filt_avg.c_not_clustered_flag(it0)
            filt_interp.flag_c_filt0_not_clustered(it & idnonnan) = filt_avg.c_not_clustered_flag(it0);
          end
          % flag when next filter event was not clustered
          if filt_avg.c_not_clustered_flag(it1)
            filt_interp.flag_c_filt1_not_clustered(it & idnonnan) = filt_avg.c_not_clustered_flag(it1);
          end
          % set linear interpolation boolean for missing fdom data
          filt_interp.flag_linear_interp_c(it & isnan(filt_interp.fdom)) = true;
        end
      else
        % linearly interpolate ag and cg if no fdom available
        filt_interp.flag_linear_interp_a(it) = true;
        filt_interp.flag_linear_interp_c(it) = true;
      end
    end
  end
  % set linear interpolation flag true when slope is negative
  filt_interp.flag_linear_interp_a(filt_interp.flag_a_negative_slope) = true;
  filt_interp.flag_linear_interp_c(filt_interp.flag_c_negative_slope) = true;

  % compute ag and cg with fCDOM
  filt_interp.a(~filt_interp.flag_linear_interp_a, :) = filt_interp.lin_a(~filt_interp.flag_linear_interp_a, :) + ...
    (filt_interp.fdom(~filt_interp.flag_linear_interp_a) - filt_interp.lin_fdom(~filt_interp.flag_linear_interp_a)) .* ...
    filt_interp.slope_interp_a(~filt_interp.flag_linear_interp_a, :);
  filt_interp.c(~filt_interp.flag_linear_interp_c, :) = filt_interp.lin_c(~filt_interp.flag_linear_interp_c, :) + ...
    (filt_interp.fdom(~filt_interp.flag_linear_interp_c) - filt_interp.lin_fdom(~filt_interp.flag_linear_interp_c)) .* ...
    filt_interp.slope_interp_c(~filt_interp.flag_linear_interp_c, :);
  
  % compute ag and cg with linear interpolation when fCDOM does not work merge filt_avg rows into 
  % filt_interp to linearly interpolate ag and cg only for missing data (in case fCDOM interpolation 
  % was available for half of filter event)
  append_filt_interp = filt_interp(end-size(filt_avg, 1)+1:end, :);
  var = filt_interp.Properties.VariableNames';
  var2mv = {'dt', 'fdom','t','s','a','c'};
  for i = 1:size(var, 1)
    if any(strcmp(var{i}, var2mv))
      append_filt_interp.(var{i}) = filt_avg.(var{i});
    else
      if islogical(append_filt_interp.(var{i}))
        append_filt_interp.(var{i}) = false(size(append_filt_interp.(var{i})));
      elseif isnumeric(append_filt_interp.(var{i}))
        append_filt_interp.(var{i}) = NaN(size(append_filt_interp.(var{i})));
      end
    end
  end
  % remove 30 seconds to appended filt_avg.dt to make sure there is no duplicate for fillmissing
  append_filt_interp.dt = datenum(datetime(append_filt_interp.dt, 'ConvertFrom', 'datenum') - seconds(30));
  filt_interp = [filt_interp; append_filt_interp];
  filt_interp = sortrows(filt_interp, 'dt');
  % linear interpolation when fdom interpolation was not used
  filt_interp.a = fillmissing(filt_interp.a, 'linear', 'SamplePoints', filt_interp.dt, 'EndValues', 'nearest');
  filt_interp.c = fillmissing(filt_interp.c, 'linear', 'SamplePoints', filt_interp.dt, 'EndValues', 'nearest');
  % remove filt_avg rows added for filling missing linear interpolation
  filt_interp(ismember(filt_interp.dt, append_filt_interp.dt), :) = [];

  % visProd3D(lambda.a, filt_interp.dt, filt_interp.a, false, 'Wavelength', false, 2);
  % title('ag fdom interpolated')
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.c, false, 'Wavelength', false, 3);
  % title('cg fdom interpolated')
  % 
  % visProd3D(lambda.a, filt_interp.dt, tot.a - filt_interp.a, false, 'Wavelength', false, 2);
  % title('ap fdom interp')
  % visProd3D(lambda.c, filt_interp.dt, tot.c - filt_interp.c, false, 'Wavelength', false, 3);
  % title('cp fdom interp')
  % 
  % visProd3D(lambda.a, filt_interp.dt, tot.a - filt_interp.lin_a, false, 'Wavelength', false, 4);
  % title('ap linear interp')
  % visProd3D(lambda.c, filt_interp.dt, tot.c - filt_interp.lin_c, false, 'Wavelength', false, 5);
  % title('cp linear interp')
  % 
  % visProd3D(lambda.a, filt_interp.dt, (tot.a - filt_interp.lin_a) - (tot.a - filt_interp.a), false, 'Wavelength', false, 4);
  % title('ap diff')
  % visProd3D(lambda.c, filt_interp.dt, (tot.c - filt_interp.lin_c) - (tot.c - filt_interp.c), false, 'Wavelength', false, 5);
  % title('cp diff')

  % % Fill missing values (start/end of analyzed period) with nearest filt_interp data
  % filt_interp.a = fillmissing(filt_interp.a, 'nearest', 'SamplePoints', filt_interp.dt);
  % filt_interp.c = fillmissing(filt_interp.c, 'nearest', 'SamplePoints', filt_interp.dt);

  % visProd3D(lambda.a, filt_interp.dt, filt_interp.a, false, 'Wavelength', false, 4);
  % title('ag fdom interpolated and nearest')
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.c, false, 'Wavelength', false, 5);
  % title('cg fdom interpolated and nearest')



  % figure; hold on
  % scatter(datetime(filt_interp.dt, 'convertfrom', 'datenum'), ag_fact_tot(:,70), 20, 'ko')
  % scatter(datetime(filt_interp.dt(~filt_interp.flag_linear_interp), 'convertfrom', 'datenum'), filt_interp.delta_af(~filt_interp.flag_linear_interp,70), 20, 'r+')
  % 
  % figure; hold on
  % scatter(filt_interp.dt, cg_fact_tot(:,70), 20, 'ko')
  % scatter(filt_interp.dt(~filt_interp.flag_linear_interp), filt_interp.delta_cf(~filt_interp.flag_linear_interp,70), 20, 'r+')

  % visProd3D(lambda.a, filt_interp.dt(it), filt_interp.a(it,:), false, 'Wavelength', false, i+1);
  % title('ag fdom interpolated')
  % visProd3D(lambda.a, filt_interp.dt(it), filt_interp.lin_a(it,:), false, 'Wavelength', false, i+2);
  % title('ag linearly interpolated')
  % 
  % visProd3D(lambda.c, filt_interp.dt(it), filt_interp.c(it,:), false, 'Wavelength', false, i+3);
  % title('cg fdom interpolated')
  % visProd3D(lambda.c, filt_interp.dt(it), filt_interp.lin_c(it,:), false, 'Wavelength', false, i+4);
  % title('cg linearly interpolated')
  % 
  % visProd3D(lambda.a, filt_interp.dt(it), filt_interp.a(it,:) - filt_interp.lin_a(it,:), false, 'Wavelength', false, i+5);
  % title('ag difference')
  % visProd3D(lambda.a, filt_interp.dt(it), filt_interp.c(it,:) - filt_interp.lin_c(it,:), false, 'Wavelength', false, i+6);
  % title('cg difference')
  %
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.a, false, 'Wavelength', false, i+1);
  % title('ag fdom interpolated')
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.c, false, 'Wavelength', false, i+2);
  % title('cg fdom interpolated')

  % yaxr = [it_filt_interp.fdom; filt_avg.fdom(it0:it1)];
  % 
  % wl_toplot = 40;
  % figure(1)
  % clf
  % subplot(1,2,1)
  % yyaxis('left')
  % hold on
  % plot_dt = datetime(it_filt_interp.dt, 'ConvertFrom', 'datenum');
  % scatter(plot_dt, ...
  %   filt_interp.a(it,wl_toplot), 40, 'filled', 'MarkerFaceColor', 'b','MarkerFaceAlpha', 0.5)
  % scatter(datetime(filt_avg.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_avg.a(:,wl_toplot), 100, 'filled', 'MarkerFaceColor', [255	193	37]/255,'MarkerFaceAlpha', 1)
  % xlim([plot_dt(1)-hours(0.5) plot_dt(end)+hours(0.5)])
  % ylabel('ag')
  % yyaxis('right')
  % hold on
  % scatter(datetime(filt_interp.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_interp.fdom, 40, 'filled', 'MarkerFaceColor', 'k', 'MarkerFaceAlpha', 0.5, 'Marker', 'd')
  % scatter(datetime(filt_avg.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_avg.fdom, 100, 'filled', 'MarkerFaceColor', 'k', 'MarkerFaceAlpha', 0.5, 'Marker', 'd')
  % xlim([plot_dt(1)-hours(0.5) plot_dt(end)+hours(0.5)])
  % ylim([min(yaxr)-0.001 max(yaxr)+0.001])
  % ylabel('fdom')
  % legend('ag fdom interpolated', 'ag filter average', 'fdom smoothed', 'Location', 'Best')
  % subplot(1,2,2)
  % yyaxis('left')
  % hold on
  % scatter(plot_dt, ...
  %   filt_interp.c(it,wl_toplot), 40, 'filled', 'MarkerFaceColor', 'r','MarkerFaceAlpha', 0.5)
  % scatter(datetime(filt_avg.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_avg.c(:,wl_toplot), 100, 'filled', 'MarkerFaceColor', [255	193	37]/255,'MarkerFaceAlpha', 1)
  % xlim([plot_dt(1)-hours(0.5) plot_dt(end)+hours(0.5)])
  % ylabel('cg')
  % yyaxis('right')
  % hold on
  % scatter(datetime(filt_interp.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_interp.fdom, 40, 'filled', 'MarkerFaceColor', 'k', 'MarkerFaceAlpha', 0.5, 'Marker', 'd')
  % scatter(datetime(filt_avg.dt, 'ConvertFrom', 'datenum'), ...
  %   filt_avg.fdom, 100, 'filled', 'MarkerFaceColor', 'k', 'MarkerFaceAlpha', 0.5, 'Marker', 'd')
  % xlim([plot_dt(1)-hours(0.5) plot_dt(end)+hours(0.5)])
  % ylim([min(yaxr)-0.001 max(yaxr)+0.001])
  % ylabel('fdom')
  % legend('cg fdom interpolated', 'cg filter average', 'fdom smoothed', 'Location', 'Best')

  % visProd3D(lambda.a, filt_interp.dt, filt_interp.a, false, 'Wavelength', false, 23);
  % title('ag fdom interpolated')
  % visProd3D(lambda.a, filt_interp.dt, filt_interp.lin_a, false, 'Wavelength', false, 24);
  % title('ag linearly interpolated')
  % 
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.c, false, 'Wavelength', false, 25);
  % title('cg fdom interpolated')
  % visProd3D(lambda.c, filt_interp.dt, filt_interp.lin_c, false, 'Wavelength', false, 26);
  % title('cg linearly interpolated')

end

%% 
function [filt_avg, filt_interp, filt, regress_stats] = fdom_agcg_model(filt, filt_interp, ...
  lambda, filt_avg, min_nb_pts_per_cluster, time_weight_for_cluster)

  % define reference wavelength to do regression (450nm)
  id450_a = abs(lambda.a - 450) == min(abs(lambda.a - 450));
  id450_c = abs(lambda.c - 450) == min(abs(lambda.c - 450));

  % iterate dbscan with multiple epsilon until nRMSE of all cluster < 10%
  fprintf('Time variable weight %.3f: Finding best epsilon for dbscan Kernel Density clustering of a(filt)/fCDOM(filt) ... ', time_weight_for_cluster)
  k = 1;
  epsilons = 1:-0.01:0.01;
  keep_going_a = true;
  keep_going_c = true;
  max_nb_clusters = day(datetime(max(filt_interp.dt) - min(filt_interp.dt), 'ConvertFrom', 'datenum'))*2;
  eval_cluster = table();
  eval_cluster.nRMSE_a450 = NaN(length(epsilons), max_nb_clusters);
  eval_cluster.nRMSE_c450 = NaN(length(epsilons), max_nb_clusters);
  eval_cluster.slope_a450 = NaN(length(epsilons), max_nb_clusters);
  eval_cluster.slope_c450 = NaN(length(epsilons), max_nb_clusters);
  eval_cluster.a_perc_not_clustered = NaN(length(epsilons), 1);
  eval_cluster.c_perc_not_clustered = NaN(length(epsilons), 1);
  eval_cluster.a_perc_negative_slope = NaN(length(epsilons), 1);
  eval_cluster.c_perc_negative_slope = NaN(length(epsilons), 1);
  % normalize the time variable between 0 to 0.5
  time_var = (filt.dt-min(filt.dt))/max(filt.dt-min(filt.dt)) * time_weight_for_cluster;
  data_tocluster_a = [filt.a(:, id450_a)./filt.fdom time_var]; % time_var filt.t filt.s
  data_tocluster_c = [filt.c(:, id450_c)./filt.fdom time_var]; % time_var filt.t filt.s
  if any(~isnan(data_tocluster_a(:, 1))) || any(~isnan(data_tocluster_c(:, 1)))
    % find minimum cluster size depending on the number of days loaded
    % round(size(time_var, 1) / (max(day(datetime(filt.dt-min(filt.dt)+1, 'ConvertFrom', 'datenum'))) * 24))
    while k <= length(epsilons) && keep_going_a && keep_going_c
      if keep_going_a
        filt.a_clusters = dbscan(data_tocluster_a, epsilons(k), min_nb_pts_per_cluster, 'Distance', 'euclidean');
        filt.a_clusters(filt.a_clusters == -1) = NaN;
      end
      if keep_going_c
        filt.c_clusters = dbscan(data_tocluster_c, epsilons(k), min_nb_pts_per_cluster, 'Distance', 'euclidean');
        filt.c_clusters(filt.c_clusters == -1) = NaN;
      end
      % regress a&c with fdom
      [regress_stats_a, regress_stats_c] = regress_acfilt(filt.a(:, id450_a), filt.c(:, id450_c), filt.fdom, filt.a_clusters, filt.c_clusters);
      if keep_going_a
        eval_cluster.nRMSE_a450(k, 1:size(regress_stats_a.nRMSE, 2)) = regress_stats_a.nRMSE;
        eval_cluster.slope_a450(k, 1:size(regress_stats_a.slope, 2)) = regress_stats_a.slope;
        eval_cluster.a_perc_not_clustered(k) = sum(isnan(filt.a_clusters))/size(filt,1)*100;
        if any(regress_stats_a.slope <= 0)
          eval_cluster.a_perc_negative_slope(k) = sum(sum(filt.a_clusters==find(regress_stats_a.slope<=0))) / size(filt,1)*100;
        else
          eval_cluster.a_perc_negative_slope(k) = 0;
        end
      end
      if keep_going_c
        eval_cluster.slope_c450(k, 1:size(regress_stats_c.slope, 2)) = regress_stats_c.slope;
        eval_cluster.nRMSE_c450(k, 1:size(regress_stats_c.nRMSE, 2)) = regress_stats_c.nRMSE;
        eval_cluster.c_perc_not_clustered(k) = sum(isnan(filt.c_clusters))/size(filt,1)*100;
        if any(regress_stats_c.slope <= 0)
          eval_cluster.c_perc_negative_slope(k) = sum(sum(filt.c_clusters==find(regress_stats_c.slope<=0))) / size(filt,1)*100;
        else
          eval_cluster.c_perc_negative_slope(k) = 0;
        end
      end
      if k >= 2
        % get epsilon_a
        if all(eval_cluster.slope_a450(k,1:size(regress_stats_a.slope, 2)) > 0) && ...
            all(eval_cluster.nRMSE_a450(k,1:size(regress_stats_a.nRMSE, 2)) < 25, 2)
          % epsilon_a = epsilons(k);
          keep_going_a = false;
        elseif size(regress_stats_a.slope, 2) > max_nb_clusters
          keep_going_a = false;
        end
        % get epsilon_c
        if all(eval_cluster.slope_c450(k,1:size(regress_stats_c.slope, 2)) > 0) && ...
            all(eval_cluster.nRMSE_c450(k,1:size(regress_stats_c.nRMSE, 2)) < 25, 2)
          % epsilon_c = epsilons(k);
          keep_going_c = false;
        elseif size(regress_stats_c.slope, 2) > max_nb_clusters
          keep_going_c = false;
        end
      end
      k = k + 1;
    end
    % find best epsilon_a based on minimum % data clustered with slope < 0 and minimum % data not clustered
    [~, ord_a] = sortrows([eval_cluster.a_perc_negative_slope eval_cluster.a_perc_not_clustered], [1 2]);
    epsilon_a = epsilons(ord_a(1));
    % find best epsilon_c based on minimum % data clustered with slope < 0 and minimum % data not clustered
    [~, ord_c] = sortrows([eval_cluster.c_perc_negative_slope eval_cluster.c_perc_not_clustered], [1 2]);
    epsilon_c = epsilons(ord_c(1));
    fprintf('done\n')

    % run clustering with the best epsilon
    fprintf('Kernel Density clustering of a(filt)/fCDOM(filt) ... ')
    % cluster absorption
    filt.a_not_clustered_flag = false(size(filt, 1), 1);
    filt.a_clusters = dbscan(data_tocluster_a, epsilon_a, min_nb_pts_per_cluster, 'Distance', 'euclidean');
    filt.a_not_clustered_flag(filt.a_clusters == -1) = true;
    if ~all(filt.a_clusters == -1)
      filt.a_clusters(filt.a_clusters == -1) = NaN;
      filt.a_clusters = fillmissing(filt.a_clusters, 'nearest', 'SamplePoints', filt.dt);
    end
    % cluster attenuation
    filt.c_not_clustered_flag = false(size(filt, 1), 1);
    filt.c_clusters = dbscan(data_tocluster_c, epsilon_c, min_nb_pts_per_cluster, 'Distance', 'euclidean');
    filt.c_not_clustered_flag(filt.c_clusters == -1) = true;
    if ~all(filt.a_clusters == -1)
      filt.c_clusters(filt.c_clusters == -1) = NaN;
      filt.c_clusters = fillmissing(filt.c_clusters, 'nearest', 'SamplePoints', filt.dt);
    end
    fprintf('done\n')

    % regress a&c with fdom
    regress_stats = struct();
    [regress_stats.a, regress_stats.c] = regress_acfilt(filt.a, filt.c, filt.fdom, filt.a_clusters, filt.c_clusters);
    
    % % flag data part of cluster with negative slope
    % clusters.a_negative_slope = false(size(filt, 1), 1);
    % clusters.c_negative_slope = false(size(filt, 1), 1);
    % if any(regress_stats_a.slope <= 0)
    %   clusters.a_negative_slope = filt.a_clusters == find(regress_stats_a.slope<=0);
    % end
    % if any(regress_stats_c.slope <= 0)
    %   clusters.c_negative_slope = filt.c_clusters == find(regress_stats_c.slope<=0);
    % end
  
    % % find wavelength where fdom fit and temperature fit work poorly
    % bad_fit_a = all(regress_stats.a.R2 < prctile(regress_stats.a.R2, 5, 1), 2);
    % bad_fit_c = all(regress_stats.c.R2 < prctile(regress_stats.c.R2, 5, 1), 2);
    
    %%% apply the regression to get a and c dissolved from fdom
    groups_a = unique(filt.a_clusters(~isnan(filt.a_clusters)));
    groups_c = unique(filt.c_clusters(~isnan(filt.c_clusters)));
  
    % Weight cluster for each filter event in case clusters are changing
    % during a filter event (unlikely with dbscan method but let's keep it in case)
    filt_avg = addvars(filt_avg, NaN(size(filt_avg,1),size(groups_a,1)), NaN(size(filt_avg,1),size(groups_c,1)), ...
      false(size(filt_avg.dt)), false(size(filt_avg.dt)),...
      false(size(filt_avg.dt)), false(size(filt_avg.dt)),...
      false(size(filt_avg.dt)), false(size(filt_avg.dt)),...
      'NewVariableNames', {'cluster_weighted_a', 'cluster_weighted_c',...
      'mix_cluster_a_flag', 'mix_cluster_c_flag',...
      'cluster_a_slope_flag', 'cluster_c_slope_flag',...
      'a_not_clustered_flag', 'c_not_clustered_flag'}, 'after', 'end');
    for i=1:size(filt_avg, 1)
      sel_filt = filt_avg.start(i) <= filt.dt & filt.dt <= filt_avg.end(i);
      foo = filt(sel_filt,:);
      % flag if clusters are not all the same within a single filter event
      if size(unique(foo.a_clusters(~isnan(foo.a_clusters))), 1) > 1
        filt_avg.mix_cluster_a_flag(i) = true;
      end
      if size(unique(foo.c_clusters(~isnan(foo.c_clusters))), 1) > 1
        filt_avg.mix_cluster_c_flag(i) = true;
      end
      % flag if any a_filt data was not clustered in dbscan
      if any(foo.a_not_clustered_flag)
        filt_avg.a_not_clustered_flag(i) = true;
      end
      % flag if any a_filt data was not clustered in dbscan
      if any(foo.c_not_clustered_flag)
        filt_avg.c_not_clustered_flag(i) = true;
      end
      % weight clusters proportionally for each filter events
      for j = 1:size(groups_a,1)
        filt_avg.cluster_weighted_a(i,j) = sum(foo.a_clusters == groups_a(j)) ./ sum(~isnan(foo.a_clusters));
      end
      for j = 1:size(groups_c,1)
        filt_avg.cluster_weighted_c(i,j) = sum(foo.c_clusters == groups_c(j)) ./ sum(~isnan(foo.c_clusters));
      end
    end
    % weight a slope and intercept
    filt_avg.slope_a = zeros(size(filt_avg, 1), size(regress_stats.a, 1));
    filt_avg.intercept_a = zeros(size(filt_avg, 1), size(regress_stats.a, 1));
    % Define same exponential function as in FitExp to get the slopes computed from the base and intercept parameters
    expfun = @(p, xd) p(1) * exp(p(2) * (xd - 440));
    for j = 1:size(groups_a, 1)
      % Fit exponential function to slopes vs lambda in the blue to reconstruct the slope values in the red where S/N is very low
      [y_intercp_a_slopefit, base_a_slopefit, ~, ~, ~] = FitExp(lambda.a(lambda.a < 500), ...
        regress_stats.a.slope(lambda.a < 500, j)', regress_stats.a.RMSE(lambda.a < 500, j)');
      fitted_a_slope = expfun([y_intercp_a_slopefit base_a_slopefit], lambda.a);

      % using slopes directly
      filt_avg.slope_a = filt_avg.slope_a + fitted_a_slope.*filt_avg.cluster_weighted_a(:,j);
      filt_avg.intercept_a = filt_avg.intercept_a + regress_stats.a.intercept(:,j)'.*filt_avg.cluster_weighted_a(:,j);

      % figure; 
      % subplot(2,1,1); plot(lambda.a, regress_stats.a.slope)
      % subplot(2,1,2); plot(lambda.a, regress_stats.a.intercept)

      % figure(3)
      % subplot(1, 2, 1)
      % plot(lambda.a, regress_stats.a.slope, 'bo');  hold on;   %plot your raw data
      % plot(lambda.a, fitted_a_slope, 'r-');  %plot the fit data
      % ylabel('a slope of a_{filt}-fdom correlation')
      % xlabel('\lambda')

    end
    % weight c slope and intercept
    filt_avg.slope_c = zeros(size(filt_avg, 1), size(regress_stats.c, 1));
    filt_avg.intercept_c = zeros(size(filt_avg, 1), size(regress_stats.c, 1));
    for j = 1:size(groups_c, 1)
      % Fit exponential function to slopes vs lambda in the blue to reconstruct the slope values in the red where S/N is very low
      [y_intercp_c_slopefit, base_c_slopefit, ~, ~, ~] = FitExp(lambda.c(lambda.c < 500), ...
        regress_stats.c.slope(lambda.c < 500, j)', regress_stats.c.RMSE(lambda.c < 500, j)');
      fitted_c_slope = expfun([y_intercp_c_slopefit base_c_slopefit], lambda.c);

      % using slopes directly
      filt_avg.slope_c = filt_avg.slope_c + fitted_c_slope.*filt_avg.cluster_weighted_c(:,j);
      filt_avg.intercept_c = filt_avg.intercept_c + regress_stats.c.intercept(:,j)'.*filt_avg.cluster_weighted_c(:,j);

      % figure; 
      % subplot(4,2,2); plot(lambda.c, regress_stats.c.slope)
      % subplot(4,2,4); plot(lambda.c, regress_stats.c.intercept)

      % subplot(1, 2, 2)
      % plot(lambda.c, regress_stats.c.slope, 'bo');  hold on;   %plot your raw data
      % plot(lambda.c, fitted_c_slope, 'r-');  %plot the fit data
      % ylabel('c slope of c_{filt}-fdom correlation')
      % xlabel('\lambda')

    end
    % interpolate slope and intercept cluster specific onto filt_interp.dt
    filt_interp.slope_interp_a = interp_extrap(filt_avg, filt_interp.dt, 'slope_a', [], true, 'linear', 'nearest');
    filt_interp.intercept_interp_a = interp_extrap(filt_avg, filt_interp.dt, 'intercept_a', [],  true, 'linear', 'nearest');
    filt_interp.slope_interp_c = interp_extrap(filt_avg, filt_interp.dt, 'slope_c', [],  true, 'linear', 'nearest');
    filt_interp.intercept_interp_c = interp_extrap(filt_avg, filt_interp.dt, 'intercept_c', [],  true, 'linear', 'nearest');
    % flag data for which a/fdom and/or c/fdom correlation has a negative slope
    filt_interp.flag_a_negative_slope = false(size(filt_interp.dt));
    filt_interp.flag_c_negative_slope = false(size(filt_interp.dt));
    filt_interp.flag_a_negative_slope(any(filt_interp.slope_interp_a < 0, 2)) = true;
    filt_interp.flag_c_negative_slope(any(filt_interp.slope_interp_c < 0, 2)) = true;
  else
    % regress_stats = [];
    if all(isnan(filt.fdom))
      error('All fDOM filter event data loaded are NaNs')
    end
    if all(isnan(data_tocluster_a))
      error('All "a" filter event data loaded are NaNs')
    end
    if all(isnan(data_tocluster_c))
      error('All "c" filter event data loaded are NaNs')
    end
  end
end


%% regression between a/c and other variable in filter events
function [a_reg, c_reg] = regress_acfilt(a, c, ancillary, a_clusters, c_clusters)
  if nargin < 4
    a_clusters = ones(size(a, 1), 1);
    c_clusters = ones(size(c, 1), 1);
  end
  % robust linear regression between filt a&c and ancillary variable
  a_reg = array2table(NaN(size(a, 2), 5), 'VariableNames', {'slope', 'intercept', 'RMSE', 'nRMSE', 'R2'});
  c_reg = array2table(NaN(size(c, 2), 5), 'VariableNames', {'slope', 'intercept', 'RMSE', 'nRMSE', 'R2'});
  groups_a = unique(a_clusters(~isnan(a_clusters)));
  groups_c = unique(c_clusters(~isnan(c_clusters)));
  for i = 1:size(a, 2)
    for j = 1:size(groups_a, 1)
      id_grp_a = a_clusters == groups_a(j);
      % regress variable with filtered water absorption
      [stats.b, stats.stats] = robustfit(ancillary(id_grp_a), a(id_grp_a, i));
      [~, MSGID] = lastwarn();
      warning('off', MSGID)
      a_reg.slope(i, j) = stats.b(2);
      a_reg.intercept(i, j) = stats.b(1);
      a_reg.RMSE(i, j) = stats.stats.robust_s;
      % calculated normalize RMSE by average and R2
      id_nonan = all(~isnan([a(:, i) ancillary]), 2);
      a_reg.nRMSE(i, j) = a_reg.RMSE(i, j)/abs(mean(a(id_grp_a & id_nonan, i), 'omitnan'))*100;
      a_reg.R2(i, j) = corr(a(id_grp_a & id_nonan, i), stats.b(1) + stats.b(2) * ancillary(id_grp_a & id_nonan))^2;
    end
    for j = 1:size(groups_c, 1)
      id_grp_c = c_clusters == groups_c(j);
      % regress variable with filtered water attenuation
      [stats.b, stats.stats] = robustfit(ancillary(id_grp_c), c(id_grp_c, i));
      [~, MSGID] = lastwarn();
      warning('off', MSGID)
      c_reg.slope(i, j) = stats.b(2);
      c_reg.intercept(i, j) = stats.b(1);
      c_reg.RMSE(i, j) = stats.stats.robust_s;
      % calculated normalize RMSE by average and R2
      id_nonan = all(~isnan([c(:, i) ancillary]), 2);
      c_reg.nRMSE(i, j) = c_reg.RMSE(i, j)/abs(mean(c(id_grp_c & id_nonan, i), 'omitnan'))*100;
      c_reg.R2(i, j) = corr(c(id_grp_c & id_nonan, i), stats.b(1) + stats.b(2) * ancillary(id_grp_c & id_nonan))^2;
    end
  end
end

