function [timebin,depthbin,beamcbin,dspkbeamcbin,spkbeamcbin,smoothbeamcfinal,parkdepth,sampperiod] = ...
    beamtransOST_QC(time,depth,beamc,negjumpthres,window_h)

%%% This function preprocesses and quality controls the data collected by
%%% a beam transmissometer used as an "optical sediment trap" sensor on a 
%%% drifting platform such as a profiling float.
%%% 
%%% Constants and thresholds may need adjustment depending on platform 
%%% behavior and sensor configuration, and the user is strongly urged to 
%%% conduct sensitivity tests with real or synthetic data from their 
%%% application to optimize this code.  
%%%
%%% Despiking procedure is based on N. Briggs et al., Deep Sea Research 
%%% Part I: Oceanographic Research Papers. 58, 1031–1039 (2011).
%%%
%%% This software is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
%%% https://creativecommons.org/licenses/by-sa/4.0/
%%%
%%% If you use this software in your own work, please acknowledge the most
%%% recent version of the accompanying publication referenced alongside
%%% this code on Zenodo (doi:10.5281/zenodo.7595966)
%%% February 2023, margaret.estapa@maine.edu

% INPUTS
% 
% time = vector of times [d]
% depth = vector of depths [m]
% beamc = vector of beam attenuation [1/m].  For a transmissometer,
%    attenuance can be obtained by multiplying beamc by the transmissometer
%    pathlength.  This script works with apparent beam attenuation rather than
%    attenuance for maximum compatibility.
% 
% negjumpthres = scalar rate of change of beam attenuation [1/(m*d)] below
%    which data should be treated as a "negative jump" and corrected.
%    Suggested starting value:  -0.2 [1/(m*d)]
%
% window_h = scalar smoothing window size in units of hours. 
%    Suggested starting value:  1.5 [h]
    
% OUTPUTS
% 
% timebin = Nx1. [d]. Regular time interval 
% depthbin = Nx2. [m]. Mean platform depth (col1) and stdev during sample
%     (col2). Value = NaN if platform fails QC tests at that time point.
% beamcbin = Nx2. [1/m]. Mean beamc (col1) and stdev during sample (col2).
%     Value = NaN if beamc fails QC tests at that time point.
% dspkbeamcbin.  Nx1. [1/m]. Beam attenuation despiked with a running median
%     followed by running mean (see below; after Briggs et al. 2011)
% spkbeamcbin.  Nx2. [1/m]. Difference between beamcbin and dspkbeamcbin.
% smoothbeamcfinal. Nx1. [1/m]. Despiked beamc after removing negative jumps and
%     smoothing with a median filter. To compute attenuance flux [m2/m2/d], 
%     compute the gradient of this vector with respect to time and multiply 
%     by the sensor pathlength.
% parkdepth.  Scalar. [m].  Computed median platform depth during deployment
% sampperiod.  Scalar. [d]. Computed time interval of timebin

%%%%%%%%%%%%%%%%%

% Set up constants and thresholds

burstsdthres = 0.005;  %max threshold for sd of beam-c values in one burst [1/m]
maxbeamcthres = Inf; % maximum beam-c increase allowed during deployment [1/m] 
depththres = 25;  % [m]  outside +/- this depth window do not use data
meanwindow = window_h/24;  % Convert window to units of days
medianwindow = window_h/24;  % Convert window to units of days

%%%%%%%%%%%%%%

% Preliminary QC

% Average the sample "bursts" recorded during drift
% Binaverager function is at the end of this m-file.
% Outputs have two columns: 1st column mean, 2nd column stdev
% If platform did not sample in bursts then this does nothing
[timebin,depthbin,beamcbin,sampperiod] = binaverager(time,depth,beamc); 

medianwindowN = medianwindow./sampperiod; %number of sampling periods in half the median window
meanwindowN = meanwindow./sampperiod;  %number of sampling periods in half the mean window

% Set to NaN any bursts with too much noise 
beamcbin(find(beamcbin(:,2) > burstsdthres),:) = NaN;

% Set to NaN bursts that are too high relative to baseline 
beamcbin(find(beamcbin(:,1)>(min(beamcbin(:,1))+maxbeamcthres)),:)=NaN;

% Set to NaN if the float was too deep or shallow
parkdepth = median(depthbin(:,1),'omitnan');
beamcbin(union(find(depthbin(:,1)>depththres+parkdepth),...
    find(depthbin(:,1)<parkdepth-depththres)),:) = NaN;

%%%%%%%%%%%%%%

% Despike the data and remove negative excursions

% Compute a moving median followed by a mean to remove spikes (Briggs et al
% 2011)  Spikes are transient and thus do not contribute to attenuance rate of change,
% but they may contain useful information so save them.
dspkbeamcbin = movmedian(beamcbin(:,1),medianwindowN,'omitnan');
dspkbeamcbin = movmean(dspkbeamcbin,meanwindowN);  
spkbeamcbin = [beamcbin(:,1)-dspkbeamcbin,beamcbin(:,2)];

%%% The section below generates a cumulative correction to eliminate rapid
%%% decreases in attenuance (more negative than "negjumpthres").  This
%%% is because there is no physical reason for a decrease in attenuance 
%%% that is only due to passively sinking particles, 
%%% except if material that has already accumulated on the window is
%%% removed (e.g., it fell off, was resuspended, or something ate it).

% Find the indices of points where negative jumps are occurring -- where
% rate of change is less than the negative threshold.
grad = gradient(dspkbeamcbin,timebin); 
njinds = find(grad<negjumpthres);

% Step through the vector, and set injump = 1 when a negative
% jump starts.  When the slope next exceeds the threshold, set injump = 0.  
% For this interval, the correction is the opposite of the mean slope
% during the interval. It is an ugly loop but the only way I could think
% about it.
cmlcorr = zeros(size(timebin));
injump = 0;
for i = 1:length(timebin)
    if ismember(i,njinds) && injump == 0  % Jump has just started
        injump = 1;
        startjump = i;
    else    % Check if jump has just ended, otherwise do nothing
        if ~ismember(i,njinds) && injump == 1 % Jump has just ended
            injump = 0;
            % Compute the necessary correction for the jump that just
            % ended, except can't correct a neg. jump in progress at start
            if i > 1 && startjump > 1  
                  cmlcorr(startjump:i) = -(gradient(dspkbeamcbin(startjump:i)));  
            end
        end
    end
end

% Take the cumulative sum of all the corrections
cmlcorr = cumsum(cmlcorr,'omitnan');  

% Add the cumulative correction to the beam attenuation timeseries
beamcfinal = dspkbeamcbin+cmlcorr;

% Smooth one more time to make the derivative nice (same running mean window 
% as was used after despiking).  movmean is centered.
smoothbeamcfinal = movmean(beamcfinal,meanwindowN); 

%%%%%%%%%%%%%%%

% Make a figure.  Commment out if not desired.

% This figure shows all the correction steps
figure;
set(gcf,'Units','Inches','Position',[1 1 5 4],'Papersize',[5 4]);
plot(timebin,beamcbin(:,1),'k');
hold on;plot(timebin,dspkbeamcbin,'r');
hold on;plot(timebin,beamcfinal,'b');
hold on;plot(timebin,smoothbeamcfinal,'m');
hold on;plot(time,beamc,'.','color',[0.5 0.5 0.5],'markersize',2);
legend({'Original after QC','Despiked','Neg. jumps removed','Final, smoothed','Original raw'},'location','northwest');
ylabel('Apparent beam attenuation (m^{-1})');
datetick('x','keeplimits','keepticks');
% % Overlay depth plot
% ah = axes('position',get(gca,'position'),'color','none','xlim',get(gca,'xlim'),...
%     'xtick',[],'yaxislocation','right','ydir','reverse','ycolor',[0.2 0.7 0]);
% hold on;
% plot(time,depth,'-','color',[0.2 0.7 0]);
% plot(timebin,depthbin(:,1),'-','color',[0.2 0.7 0],'linewidth',1.5);
% plot(timebin,ones(size(timebin)).*(parkdepth-depththres),'--','color',[0.2 0.7 0]);
% plot(timebin,ones(size(timebin)).*(parkdepth+depththres),'--','color',[0.2 0.7 0]);
% ylabel('Depth +/- threshold');

%%%%%%%%%%%%%%%%%%

% Function below averages bursts to return means on a regular time
% coordinate plus the standard deviation for each bin.  If the data were
% not sampled in "bursts" then this doesn't really have much effect.

function [timebin,depthbin,beamcbin,sampperiod] = binaverager(time,depth,beamc)

sampperiod = max(unique(diff(time)));
timebinedges = [min(time):sampperiod:max(time)]';
timebin = mean([timebinedges(1:end-1),timebinedges(2:end)],2);
depthbin = NaN(length(timebin),2);
beamcbin = NaN(length(timebin),2);
[~,~,inds] = histcounts(time,timebinedges);
for i = 1:length(timebin)
    if ~isempty(find(inds==i))
        depthbin(i,1) = mean(depth(find(inds==i)),'omitnan');
        depthbin(i,2) = std(depth(find(inds==i)),'omitnan');
        beamcbin(i,1) = mean(beamc(find(inds==i)),'omitnan');
        beamcbin(i,2) = std(beamc(find(inds==i)),'omitnan');
    end
end
% These lines added to deal with occasional APEX gaps in pressure record
gooddepth = find(~isnan(depthbin(:,1)));
depthbininterp(:,1) = interp1(timebin(gooddepth),depthbin(gooddepth,1),timebin,'linear');
depthbininterp(:,2) = interp1(timebin(gooddepth),depthbin(gooddepth,2),timebin,'nearest');
depthbin = depthbininterp;
