function ox=oxyg_eq(p_equation, p_slope, sbecoef, sbecoefd, params)
%calculate CTD oxygen from T/S/P and parameters
%p_equation=
% 1: Johnson /PMEL: Soc*(Oxv+voffset+T20_DO*exp(DI*P+D2*T).*dVdt).*osat.*exp(Tcor*T).*exp(Pcor*P./(273.15+T))
%    Fits Soc, Voffset, Tau20_D0, Tcorr and E.  Here Tau20 and D0 have been set to a single parameter.
% 2: Seabird calibration sheet 1: Soc.*(Oxv+voffset).*(1.0+A*T+B*T.^2+C*T.^3).*osat.*exp(pcor*P./(T+273.15))
% Other Seabird equation: Soc.*(Oxv+voffset).*exp(Tcor*T).*btl.dosat.*exp(Pcor*P)
%
%p_slope=1 uses same but with station-dependant slope (Soc)
%A. Ganachaud, 2014
%params is a structure that must contain appropriate fields
%  sbecoef is as follows:
%  sbecoef(1)=soc (slope)
%  sbecoef(2)=voffset (bias)
%  sbecoef(3)=tau_20 (lag) (if needed) 
%  sbecoef(4)=tcor (if needed)
%  sbecoef(5)=pcor (if needed) or E (if needed)

%because we are not sure that params includes dnstat if p_slope=0
if p_slope==1
    sbecoef1=sbecoef(1)+sbecoef(6)*params.dnstat;
else
    sbecoef1=sbecoef(1)*ones(size(params.stat));
end

switch p_equation
    case 1 %PMEL equation
       ox=sbecoef1.*...
           (params.doxv+sbecoef(2)+sbecoef(3)*exp(sbecoefd.d1*params.dpres+sbecoefd.d2*params.dtemp68).*params.dsdvdt).*...
           params.dosat.*...
           exp(sbecoef(4)*params.dtemp68).*exp(sbecoef(5)*params.dpres./(273.15+params.dtemp68));
    case 2 %SEABIRD 1 (ABC coefficients)
       ox=sbecoef1.*...
           (params.doxv+sbecoef(2)).*(1.0 + sbecoefd.A*params.dtemp68 + sbecoefd.B*params.dtemp68.^2 + sbecoefd.C*params.dtemp68.^3).*...
           params.dosat.*...
           exp(sbecoef(5)*params.dpres./(params.dtemp68+273.15));
    case 3 %SEABIRD 2 (pcor equation)
       ox=sbecoef1.*...
           (params.doxv+sbecoef(2)).*...
           exp(sbecoef(3)*params.dtemp68).*...
           params.dosat.*...
           exp(sbecoef(4)*params.dpres);
    otherwise
        error('check p_equation')
end    