# Salinity, temperature, pressure -> depth equations in Python

def depth_StandardOcean(mbar, lat=23):
    # Following Algorithms for computation of fundamental properties
    # of seawater: Unesco technical papers in marine science 44
    p = mbar / 100.0
    x = sin(lat/180.0*pi)**2
    C = (((-1.82e-15*p + 2.279e-10)*p - 2.2512e-5)*p + 9.72659)*p
    gr = 9.780318*(1.0 + (5.2788e-3 + 2.36e-5 * x)*x) + 1.092e-6*p
    z = C / gr
    return z

def depth_Freshwater(mbar):
    return mbar / 100.0 * 1.019716

def depth(S,t,cb,lat=23):
    # Following Algorithms for computation of fundamental properties
    # of seawater: Unesco technical papers in marine science 44
    # EOS-80, valid for S=0...42 ppt, t=-2...40 degC, p=0...100000 millibars
    b = cb/1000.0 # Convert to bars
    db = cb / 100.0 # Convert to decibars
    x = sin(lat/180.0*pi)**2

    # Pure water terms of the secant bulk modulus
    Kw = 19652.21 + 148.4206*t - 2.327105*t**2 + 1.360477e-2*t**3 - 5.155288e-5*t**4
    Aw = 3.239908 + 1.43713e-3*t + 1.16092e-4*t**2 - 5.77905e-7*t**3
    Bw = 8.50935e-5 - 6.12293e-6*t + 5.2787e-8*t**2
    
    K0 = Kw + (54.6746 - 0.603459*t + 1.09987e-2*t**2 - 6.1670e-5*t**3)*S \
         + (7.944e-2 + 1.6483e-2*t - 5.3009e-4*t**2)*S**(1.5)
    A = Aw + (2.2838e-3 - 1.0981e-5*t - 1.6078e-6*t**2)*S + 1.91075e-4*S**(1.5)
    B = Bw + (-9.9348e-7 + 2.0816e-8*t + 9.1697e-10*t**2)*S
    
    # Density of Standard Mean Oceam Water (SMOW)
    rw = 999.842594 + 6.793952e-2*t - 9.095290e-3*t**2 + 1.001685e-4*t**3 \
         - 1.120083e-6*t**4 + 6.536332e-9*t**5
    r0 = rw + (8.24493e-1 - 4.0899e-3*t + 7.6438e-5*t**2 - 8.2467e-7*t**3 \
               + 5.3875e-9*t**4)*S \
               + (-5.72466e-3 + 1.0227e-4*t - 1.6546e-6*t**2)*S**(1.5) \
               + 4.8314e-4*S**2
    #K = K0 + A*b + B*b**2
    #r = r0 / (1 - b / K)
    
    # Specific volume
    V0 = 1.0 / r0
    
    R = sqrt(A**2 - 4*B*K0)
    C = V0 * (b - 1/(2*B) * log((K0+A*b+B*b**2)/K0) + A/(2*B*R) * log((1+(2*B/(A-R))*b)/(1+(2*B/(A+R))*b))) * 100000.0
    gr = 9.780318*(1.0 + (5.2788e-3 + 2.36e-5 * x)*x) + 1.092e-6*db
    z = C/gr
    return z