
% 
% Computes sound velocity from in-situ temperature, salinity, and pressure
% 
% Method 1: UNESCO (Chen and Millero) equation
% Method 2: Mackenzie formula
% 
% 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 SV = soundspeed_from_temppsal(temp,psal,pres,lat,method)

if method == 1  % Chen and Millero's formula
    
    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 ;
    
    SV = cw + a.*psal + b.*(psal.^1.5) + d.*(psal.^2);
    
elseif method == 2  % Mackenzie's formula
    
    % Depth from pressure, without dynamic height
    z = -gsw_z_from_p(pres,lat);

    % Salinity from sound velocity, temperature, and depth
    SV = 1448.96 + 4.591.*temp - 5.304e-2.*(temp.^2) + 2.374e-4.*(temp.^3) + ...
         1.340*(psal-35) + 1.630e-2.*z + 1.675e-7.*(z.^2) -1.025e-2.*temp.*(psal-35) - ...
         7.139e-13.*temp.*(z.^3);
    
else
    
    error('Only two methods are available. 1: Chen-Millero (UNESCO), 2: Mackenzie');
    
end
