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.
%% 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
% Tps=[repmat(NaN,pshift,1); T([1:length(T)-pshift])];
% Sps=[repmat(NaN,pshift,1); S([1:length(S)-pshift])];
Tfs=[repmat(NaN,pshift,1); T([1:length(T)-pshift])];
Sfs=[repmat(NaN,pshift,1); T([1:length(T)-pshift])];
% Tfs = interp1(time,T,time-Olag(1)); Sfs = interp1(time,S,time-Olag(1));

%% Advance oxygen voltages and their appropriate T/S relative to pressure
% ocV=[repmat(NaN,tshift,1); ocV([1:length(ocV)-tshift])];
% otV=[repmat(NaN,tshift,1); otV([1:length(otV)-tshift])];
% Tfs=[repmat(NaN,tshift,1); Tps([1:length(Tps)-tshift])];
% Sfs=[repmat(NaN,tshift,1); Sps([1:length(Sps)-tshift])];
% Tfs=T; Sfs=S;

%%----------------------------------
%% 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);

%%------------------------------------------------
%% Oxygen concentration assuming min T = min To
%%------------------------------------------------
% if strcmp(pltit,'y')==1;
%  Tdif = min(Tfs)-min(To); To=To+Tdif;
%  O2c = [SOc * (oc+tau*docdt) + BOc] .* OXSAT .* exp(tcor*[Tfs+wt*(To-Tfs)] + pcor*P);
%  plot(time,To,'g'); subplot(212)
%  plot(O2,(O2-O2c)./O2,'b.'); xlabel('O2'); ylabel('O2 % difference assuming min To = min T');
%  disp('Pause'); pause
% end

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
