function O2 = calc_O2(ocV,otV,T,S,P,m,b,k,c,SOc,BOc,time,Olag,units);
%% To calculate Oxygen from sbe voltages for oxygen current 
%%  and membrane temperature 
%% Reference: W.B. Owens and R.C. Millard, 1985. A new algorithm for CTD oxygen
%%	calibration, J. Phys. Oceanogr. 15:621-631. 
%%	Also, described in SeaBird application note.
pltit = 'n'; disp(['Pressure should be in dB']);

%% Frame shifting for oxygen calculations, plumbing lag and sensor response
%% SeaBird advances oxygen 1 to 5 seconds relative to pressure.
%% Oxygen calculations account for both lags, 
%% but pressure is not directly used in the oxygen calculation. 
%% Due to plumbing lag, use T/S (t-plumbing lag = shift forward) ; 
%% Shift oxygen voltages forward (t-total lag = shift forward);
%% Define lags
Ofreq = mean(diff(time)); 	%  oxygen sampling frequency 
plag = Olag(1);	pshift = round(plag/Ofreq); 	% plumbing lag and shift
rlag = Olag(2); rshift = round(rlag/Ofreq); 	% response lag and shift
tlag = plag+rlag; tshift = round(tlag/Ofreq);	% total lag and shift
%% T and S only need pluming lag shift (+ chamber mixing lag)
%% Estimate is 3 seconds
Tfs=[repmat(NaN,pshift,1); T([1:length(T)-pshift])];
Sfs=[repmat(NaN,pshift,1); T([1:length(T)-pshift])];

%%----------------------------------
%% Oxygen current
%%----------------------------------
oc = 1e6*(m*ocV+b);			% oxygen current (uA);
disp(['Minimum of oxygen current ' num2str(min(oc)) ' uA']);
disp(['Maximum of oxygen current ' num2str(max(oc)) ' uA']);
disp(['Mean of oxygen current ' num2str(mean(oc)) ' uA']);
docdt = [NaN;[diff(oc)./diff(time)]];	% slope of oxygen current (uA/sec);

%%----------------------------------
%% Oxygen membrane temperature
%%----------------------------------
To = k*otV+c;

if strcmp(pltit,'y')==1;
 figure(1); clf; giffig
 subplot(211); plot(time,To,'b'); hold on; plot(time,Tfs,'r--');
 legend('Oxygen sensor internal temperature','water temperature');
 disp('Pause'); pause
end

%%----------------------------------
%% Oxygen saturation: f(T,S); ml/l
%%----------------------------------
TK = 273.15+Tfs;	% degrees Kelvin
A1 = -173.4292; A2 = 249.6339; A3 = 143.3483; A4 = -21.8492;
B1 = -0.033096; B2 = 0.014259; B3 =  -0.00170;
OXSAT = exp(A1 + A2*(100./TK) + A3*log(TK/100) + A4*(TK/100) + [Sfs .* (B1 + B2*(TK/100) + (B3*(TK/100).*(TK/100)))] );

%%----------------------------------
%% Oxygen concentration
%%----------------------------------
%% Constants
tcor = -.033; pcor = 1.5e-4; tau=2; wt = .67;
%% Frame-shifted temperature
O2 = [SOc * (oc+tau*docdt) + BOc] .* OXSAT .* exp(tcor*[Tfs+wt*(To-Tfs)] + pcor*P);

%%----------------------------------
%% Frame Shifted Oxygen concentration
%% to account for total lag
%% (T and S for O2 calculation only need pluming lag shift)
%%----------------------------------
Ofs=[O2([tshift+1:length(O2)]); repmat(NaN,tshift,1)];
O2=Ofs;

if strcmp(units,'umolkg')==1
 %%----------------------------------
 %% Convert to umol/kg  
 %%----------------------------------
 %% SeaBird equations are for ml/l computations
 %%  Can convert OXSAT at atmospheric pressure to mg/l by 1.4276
 %%  Convert dissolved O2 to mg/l using density of oxygen = 1.4276 kg/m^3
 dens=sw_dens(Sfs,Tfs,P);
 O2 = (O2 * 1.4276) .* (1e6./(dens*32));
end
