Lead Ball in Mariana Trench

#!/usr/bin/env python
# encoding: utf-8
"""
lead_ball.py

Created by hariya on 2011-05-11.
Copyright (c) 2011 hariya. All rights reserved.
"""

import numpy as np
from scipy.special import cbrt

def temperature(s,t,p):
  # Sea water potential temperature computation
  # Ref: HL Bryden - Deep Sea Research and Oceanographic Abstracts, 1973

  p2 = p*p
  p3 = p2*p

  t2 = t*t
  t3 = t2*t

  ss = s - 35.0

  temp = t - p*(3.65e-4 + 8.32e-5*t - 5.41e-7*t2 + 4.03e-9*t3) \
           - p*ss*(1.74e-5 - 2.98e-7*t) - p2*(8.93e-7 - 3.16e-8*t + 2.20e-10*t2) \
           + 4.11e-9*p2*ss - p3*(-1.61e-10 + 5.05e-12*t)
  return temp

def density(s,t,p):
  # Sea water density computation
  # Ref: GL Mellor - Journal of Atmospheric and Oceanic Technology, 1991
  t2 = t*t
  t3 = t*t2
  t4 = t*t3
  t5 = t*t4

  c = 1449.1 + 0.0821*p + 4.55*t -0.045*t2 + 1.34*(s - 35.0)
  pcc = p/(c*c)

  rho = -0.157406 + 6.793952e-2*t - 9.095290e-3*t2 + 1.001685e-4*t3 - 1.120083e-6*t4 + 6.536332e-9*t5 \
      + (0.824493 - 4.0899e-3*t + 7.6438e-5*t2 - 8.2467e-7*t3 + 5.3875e-9*t4)*s \
      + (-5.72466e-3 + 1.0227e-4*t - 1.6546e-6*t2)*abs(s)**1.5 + 4.8314e-4*s*s

  rho = 1000.0 + rho + 1.0e5*pcc*(1.0 - 2.0*pcc)

  return rho

def main():
  # Lead ball radius computation
  pb_m = 1.0
  pb_v = pb_m/11340.0
  pb_r = cbrt(3*pb_v/(4*np.pi))

  print "Radius of lead ball is ",pb_r," m."

  # Guesses for salinity, and surface temperature
  sal = 35.0 # Salinity in parts per thousand
  sst = 25.0 # Sea surface temperature in C

  # Initialize some values
  dt = 0.01
  time = 0.0
  depth = 0.0
  vel = 0.0
  rho = 1025.0

  # Loop till we hit the bottom
  while depth > -10911.0:
    p = -9.8*rho*depth*1.0e-5
    tt = temperature(sal,sst,p)
    rho = density(sal,tt,p)

    buoy = 9.8*(11340 - rho)*pb_v
    drag = 6*np.pi*1.08e-3*vel*pb_r

    acc = (buoy - drag)

    time = time + dt
    vel = vel + acc*dt
    depth = depth - vel*dt

    print time,depth,vel,p,tt,rho,buoy,drag,acc


if __name__ == '__main__':
      main()