function [doxy_tcorr] = timecorr_optode_B14(doxy, tc, secs, optode_type, blthickness)

% Applies time correction for optodes based on equation 37 in Bittig et al.
% 2014.
%
% doxy is an array for oxygen, tc is an array from temperature, and secs is an array for time
% elapsed in seconds. Expected order is from surface at i = 1, and bottom at i = end. So secs should
% be be declining. This assumes that this is an ascending profile. 
%
% optode_type must be 'sbe63', 'aand', or 'aandF'. 
%
% blthickness is the thickness of the boundary layer, and is in um. 
%
% for pumped SBE63 on gliders and floats: 34 um (Bittig et al. 2017)
%
% yui Takeshita
% May-3-2020
%
% Update 11/25/2020
%
% allows input for blthickness to be an array or scalar. Use scalar for SBE63 where BLT is constant
% throughout the profile because it's pumped. Array for varying BLT because of changing float ascent
% speed for non-pumped Aanderaas. All changes were made in tau_optode_bittig17

% calculate tau based on optode type and temperature
tau = tau_optode_bittig17(blthickness, tc, optode_type);

% model parameters
dt = -diff(secs); % time interval between measurements
b = (1 + 2*tau(2:end)./dt).^-1; % Eq 34 in Bittig et al. 2014
a = 1-2*b; % Eq 34 in Bittig et al. 2014

% Inputs into equation 37, Bittig et al. 2014
% cfi = cfilt at time i, and cfiplus1 = cfilt at time i+1
cfi = doxy(2:end);
cfiplus1 = doxy(1:end-1);

% this is the left side of equation 37: the averaged c_O2 between i & i+1
cin_avg = 1./(2.*b).*(cfiplus1 - a.*cfi);

% need to interpolate above value onto original time stamp

% calculate the associated time stamp for cin_avg
t0 = secs(2:end);
t1 = secs(1:end-1);
% new timestamp
tavg = (t0+t1)./2;

% for interpolation only use non-NaN data
iuse = ~isnan(tavg) & ~isnan(cin_avg);
doxy_tcorr = interp1(tavg(iuse), cin_avg(iuse), secs, 'linear', 'extrap');
% 
% % ========= Correction for Takeshita et al. 2013 ==================
% %
% % use the same tau calculated from Bittig
% %
% % calculate do2/dt
% ddoxy = diff(s.doxy(:,divenum));
% do2dt = ddoxy./dt; 
% % calculate seconds passed since start
% t = s.secs(2:end,divenum) - nanmin(s.secs(:,divenum));
% 
% doxy_ = s.doxy(2:end,divenum) - tau.*do2dt.*(1 - exp(-t./tau));
% 
% % try different running mean sizes
% do2dtrm1 = runmean(do2dt,1,1,'zero');
% do2dtrm2 = runmean(do2dt,2,1,'zero');
% 
% doxy_1 = s.doxy(2:end,divenum) - tau.*do2dtrm1.*(1 - exp(-t./tau));
% doxytcorr_t13 = s.doxy(2:end,divenum) - tau.*do2dtrm2.*(1 - exp(-t./tau));
% doxytcorr_t13 = [s.doxy(1,divenum); doxytcorr_t13];
% 


end
