
% 
% Computes salinity from in-situ sound velocity, temperature, and pressure or depth
% 
% 
% INPUTS:
% - SV:   sound velocity (m/s)
% - temp: temperature (degC)
% - pres: pressure (dbar)
% - lat:  latitude (degN)
% - method: 1: invert UNESCO (Chen and Millero) equation (nonlinear in S)
%           2: invert Mackenzie formula (linear in S)
% 
% 
% References:
% C-T. Chen and F.J. Millero, "Speed of sound in seawater at high
% pressures", J. Acoust. Soc. Am., vol. 62, No. 5,1977, pp. 1129-1135
% N.P. Fofonoff and R.C. Millard Jr, "Algorithms for computation of
% fundamental properties of seawater", UNESCO technical papers in marine
% science, No. 44, Division of Marine Sciences, 1983
% K.V. Mackenzie, "Nine-term Equation for Sound Speed in the Oceans", 
% JASA, vol. 70, No. 3, Sept. 1981, pp. 807 - 812
% K.V. Mackenzie, "Discussion of Sea Water Sound-speed Determinations",
% JASA, vol. 70, No. 3, Sept. 1981  pp. 801 - 806
% 
% 


function psal = MVP_salinity_from_soundspeed(SV,temp,pres,lat,method)

%% Compute the linear solution with Mackenzie's formula
% Depth from pressure, without dynamic height
z = -gsw_z_from_p(pres,lat);

% Salinity from sound velocity, temperature, and depth
psalMack = (SV - 1448.96 - 4.591.*temp + 5.304e-2.*(temp.^2) - 2.374e-4.*(temp.^3) - ...
            1.630e-2.*z - 1.675e-7.*(z.^2) + 7.139e-13.*temp.*(z.^3) ) ./ ...
           ( 1.34 - 1.025e-2.*temp ) + 35;
   

%% If the more accurate Chen and Millero's formula is asked, use Mackenzie
%%  to choose among the roots
if method == 1
    
    pres = pres/10; % dbar -> bar
    
    c00 =  1402.388  ; c01 =  5.03830    ; c02 = -5.81090e-2 ; c03 =  3.3432e-4  ; c04 = -1.47797e-6 ; c05 = 3.1419e-9 ; 
    c10 =  0.153563  ; c11 =  6.8999e-4  ; c12 = -8.1829e-6  ; c13 =  1.3632e-7  ; c14 = -6.1260e-10 ; 
    c20 =  3.1260e-5 ; c21 = -1.7111e-6  ; c22 =  2.5986e-8  ; c23 = -2.5353e-10 ; c24 =  1.0415e-12 ; 
    c30 = -9.7729e-9 ; c31 =  3.8513e-10 ; c32 = -2.3654e-12 ; 
    a00 =  1.389     ; a01 = -1.262e-2   ; a02 =  7.166e-5   ; a03 =  2.008e-6   ; a04 = -3.21e-8    ; 
    a10 =  9.4742e-5 ; a11 = -1.2583e-5  ; a12 = -6.4928e-8  ; a13 = 1.0515e-8   ; a14 = -2.0142e-10 ; 
    a20 = -3.9064e-7 ; a21 =  9.1061e-9  ; a22 = -1.6009e-10 ; a23 = 7.994e-12   ; 
    a30 =  1.100e-10 ; a31 =  6.651e-12  ; a32 = -3.391e-13  ; 
    b00 = -1.922e-2  ; b01 = -4.42e-5    ; b10 =  7.3637e-5  ; b11 = 1.7950e-7 ; 
    d00 =  1.727e-3  ; d10 = -7.9836e-6  ; 

    cw = (c00 + c01*temp + c02*temp.^2 + c03*temp.^3 + c04*temp.^4 + c05*temp.^5) + ...
         (c10 + c11*temp + c12*temp.^2 + c13*temp.^3 + c14*temp.^4).*pres         + ...
         (c20 + c21*temp + c22*temp.^2 + c23*temp.^3 + c24*temp.^4).*( pres.^2 )  + ...
         (c30 + c31*temp + c32*temp.^2).*( pres.^3 ) ;
    a  = (a00 + a01*temp + a02*temp.^2 + a03*temp.^3 + a04*temp.^4)    + ...
         (a10 + a11*temp + a12*temp.^2 + a13*temp.^3 + a14*temp.^4).*pres + ...
         (a20 + a21*temp + a22*temp.^2 + a23*temp.^3).*( pres.^2 )     + ...
         (a30 + a31*temp + a32*temp.^2).*( pres.^3 ) ;
    b  =  b00 + b01*temp + (b10 + b11*temp).*pres ;
    d  =  d00 + d10*pres ;
    
    coefs = [d(:).^2, ...
             2*a(:).*d(:)-b(:).^2, ...
             a(:).^2+2*(cw(:)-SV(:)).*d(:), ...
             2*(cw(:)-SV(:)).*a(:), ...
             (cw(:)-SV(:)).^2];
    
    psal = nan(size(SV));
    for iPt = 1 : length(SV)
        if all(~isnan(coefs(iPt,:)))
            psalTmp = roots(coefs(iPt,:));  % all 4 roots
            for iSol = 1 : 4
                if ~isreal(psalTmp(iSol))
                    psalTmp(iSol) = 999;    % set complex roots at a very large value
                end
            end
            [~, ind] = min(abs(psalTmp-psalMack(iPt))); % pick the nearest root to Mackenzie's one
            psal(iPt) = psalTmp(ind);
        end
    end
    
elseif method == 2
    
    psal = psalMack;
    
else
    
    error('Only two methods are available. 1: Chen-Millero (UNESCO), 2: Mackenzie');
    
end
