#ifndef __JACOUSTICS__JUNESCO__
#define __JACOUSTICS__JUNESCO__

#include "JMath/JConstants.hh"

/**
 * \file
 *
 * Sound velocity according UNESCO equation.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  using JMATH::PI;

  /**
   * Get gravitational constant.
   *
   * \param  phi          latitude [deg]
   * \return              g        [m/s^2]
   */
  inline double getGravity(const double phi)
  {
    const double sd  = sin(phi * PI / 180.0);
    const double s2  = sd*sd;
    const double s4  = s2*s2;

    return 9.780318 * (1.0 + 5.2788e-3*s2 + 2.36e-5*s4);
  }


  /**
   * Get sound velocity.
   *
   * K.V. Mackenzie, Nine-term equation for the sound speed in the oceans (1981)\n
   * J. Acoust. Soc. Am. 70(3), pp 807-812
   *
   * \param  D            depth [m]
   * \param  S            salinity [ppk]
   * \param  T            temperature [C]
   * \return              sound velocity [m/s]
   */
  inline double getVelocityMackenzie(const double D, const double S, const double T)
  {
    return (1448.96
	    + 4.591 * T
	    - 5.304e-2  * T*T
	    + 2.374e-4  * T*T*T
	    + 1.340 * (S - 35)
	    + 1.630e-2  * D
	    + 1.675e-7  * D*D
	    - 1.025e-2  * T * (S - 35)
	    - 7.139e-13 * T*D*D*D);
  }


  /**
   * Get sound velocity.
   *
   * A.B. Coppens, Simple equations for the speed of sound in Neptunian waters (1981)\n
   * J. Acoust. Soc. Am. 69(3), pp 862-863
   *
   * \param  D            depth [m]
   * \param  S            salinity [ppk]
   * \param  T            temperature [C]
   * \return              sound velocity [m/s]
   */
  inline double getVelocityCoppens(const double D, const double S, const double T)
  {
    const double t = T * 0.1;
    const double d = D * 1.0e-3;
    const double c = 1449.05 + 45.7 * t - 5.21 * t*t + 0.23 * t*t*t + (1.333 - 0.126*t + 0.009 * t*t)*(S - 35);

    return (c + 
	    (16.23 + 0.253*t) * d    + 
	    (0.213 - 0.1*t)   * d*d  + 
	    (0.016 + 0.0002*(S-35))*(S - 35)*t*d);
  }


  /**
   * Get sound velocity.
   *
   * C-T. Chen and F.J. Millero, Speed of sound in seawater at high pressures (1977)\n
   * J. Acoust. Soc. Am. 62(5) pp 1129-1135
   *
   * \param  P            pressure [bar]
   * \param  S            salinity [ppk]
   * \param  T            temperature [C]
   * \return              sound velocity [m/s]
   */
  inline double getVelocityUNESCO(const double P, const double S, const double T)
  {
    const double A00  =         1.389;
    const double A01  =     -1.262E-2;
    const double A02  =      7.166E-5;
    const double A03  =      2.008E-6;
    const double A04  =      -3.21E-8;
    const double A10  =     9.4742E-5;
    const double A11  =    -1.2583E-5;
    const double A12  =    -6.4928E-8;
    const double A13  =     1.0515E-8;
    const double A14  =   -2.0142E-10;
    const double A20  =    -3.9064E-7;
    const double A21  =     9.1061E-9;
    const double A22  =   -1.6009E-10;
    const double A23  =     7.994E-12;
    const double A30  =     1.100E-10;
    const double A31  =     6.651E-12;
    const double A32  =    -3.391E-13;
    const double B00  =     -1.922E-2;
    const double B01  =      -4.42E-5;
    const double B10  =     7.3637E-5;
    const double B11  =     1.7950E-7;
    const double C00  =      1402.388;
    const double C01  =       5.03830;
    const double C02  =   -5.81090E-2;
    const double C03  =     3.3432E-4;
    const double C04  =   -1.47797E-6;
    const double C05  =     3.1419E-9;
    const double C10  =      0.153563;
    const double C11  =     6.8999E-4;
    const double C12  =    -8.1829E-6;
    const double C13  =     1.3632E-7;
    const double C14  =   -6.1260E-10;
    const double C20  =     3.1260E-5;
    const double C21  =    -1.7111E-6;
    const double C22  =     2.5986E-8;
    const double C23  =   -2.5353E-10;
    const double C24  =    1.0415E-12;
    const double C30  =    -9.7729E-9;
    const double C31  =    3.8513E-10;
    const double C32  =   -2.3654E-12;
    const double D00  =      1.727E-3;
    const double D10  =    -7.9836E-6;

    const double T2 = T  * T;
    const double T3 = T2 * T;
    const double T4 = T3 * T;
    const double T5 = T4 * T;

    const double P2 = P  * P;
    const double P3 = P2 * P;

    const double Cw =
      (C00 + C01 * T + C02 * T2 + C03 * T3 + C04 * T4 + C05 * T5) +
      (C10 + C11 * T + C12 * T2 + C13 * T3 + C14 * T4) * P        +
      (C20 + C21 * T + C22 * T2 + C23 * T3 + C24 * T4) * P2       +
      (C30 + C31 * T + C32 * T2) * P3;

    const double A =
      (A00 + A01 * T + A02 * T2 + A03 * T3 + A04 * T4)     +
      (A10 + A11 * T + A12 * T2 + A13 * T3 + A14 * T4) * P +
      (A20 + A21 * T + A22 * T2 + A23 * T3) * P2           +
      (A30 + A31 * T + A32 * T2) * P3;

    const double B = B00 + B01 * T + (B10 + B11 * T) * P;

    const double D = D00 + D10 * P;

    return Cw + A*S + B*S*sqrt(S) + D*S*S;
  }


  /**
   * Get sound velocity.
   *
   * V.A. Del Grosso, New equation for the speed of sound in natural waters (with comparisons to other equations) (1974)\n
   * J. Acoust. Soc. Am 56(4) pp 1084-1091
   *
   * G.S.K. Wong and S Zhu, Speed of sound in seawater as a function of salinity, temperature and pressure (1995)\n
   * J. Acoust. Soc. Am. 97(3) pp 1732-1736
   *
   * \param  p            pressure [bar]
   * \param  S            salinity [ppk]
   * \param  T            temperature [C]
   * \return              sound velocity [m/s]
   */
  inline double getVelocityDelGrosso(const double p, const double S, const double T)
  {
    const double C000   =        1402.392;
    const double CT1    =     0.5012285E1;
    const double CT2    =    -0.551184E-1;
    const double CT3    =     0.221649E-3;
    const double CS1    =     0.1329530E1;
    const double CS2    =    0.1288598E-3;
    const double CP1    =       0.1560592;
    const double CP2    =    0.2449993E-4;
    const double CP3    =   -0.8833959E-8;
    const double CST    =   -0.1275936E-1;
    const double CTP    =    0.6353509E-2;
    const double CT2P2  =    0.2656174E-7;
    const double CTP2   =   -0.1593895E-5;
    const double CTP3   =    0.5222483E-9;
    const double CT3P   =   -0.4383615E-6;
    const double CS2P2  =   -0.1616745E-8;
    const double CST2   =    0.9688441E-4;
    const double CS2TP  =    0.4857614E-5;
    const double CSTP   =   -0.3406824E-3;

    // 1 bar = 100 kPa = 1.019716 kg/cm^2.

    const double f  = 1.019716;
    const double P  = p * f;

    const double T2 = T  * T;
    const double T3 = T2 * T;

    const double S2 = S  * S;

    const double P2 = P  * P;
    const double P3 = P2 * P;

    const double DCT   = CT1 * T + CT2 * T2 + CT3 * T3;
    const double DCS   = CS1 * S + CS2 * S2;
    const double DCP   = CP1 * P + CP2 * P2 + CP3 * P3;
    const double DCSTP = CTP * T*P + CT3P * T3*P + CTP2 * T*P2 + CT2P2 * T2*P2 + CTP3 * T*P3 + CST * S*T + CST2 * S*T2 + CSTP * S*T*P + CS2TP * S2*T*P + CS2P2 * S2*P2;
 
    return C000 + DCT + DCS + DCP + DCSTP;
  }


  /**
   * Get depth.
   *
   * Claude C. Leroy and Francois Parthiot, Depth-pressure relationships in the oceans and seas\n
   * J. Acoust. Soc. Am. 103(3), pp 1346-1352
   *
   * \param  P            pressure [MPa]
   * \param  phi          latitude [deg]
   * \param  option       option for Mediterranean
   * \return              depth    [m]
   */
  inline double getDepth(const double P, const double phi, const bool option = true)
  { 
    const double sd  = sin(phi * PI / 180.0);
    const double g   = 9.780318 * (1.0 + 5.2788e-3 * sd*sd + 2.363e-5 * sd*sd*sd*sd);
    const double z   = (9.72659e2 * P - 2.2512e-1 * P*P + 2.279e-4 * P*P*P - 1.82e-7 * P*P*P*P) / (g + 1.092e-4 * P);
    const double dz  = -7e-2 * P + 2e-3* (P*P);

    return z + (option ? dz : 0.0);
  }


  /**
   * Get pressure.
   *
   * Claude C. Leroy and Francois Parthiot, Depth-pressure relationships in the oceans and seas\n
   * J. Acoust. Soc. Am. 103(3), pp 1346-1352
   *
   * \param  z            depth    [m]
   * \param  phi          latitude [deg]
   * \param  option       option for Mediterranean
   * \return              pressure [MPa]
   */
  inline double getPressure(const double z, const double phi, const bool option = true)
  { 
    const double sd  = sin(phi * PI / 180.0);
    const double g   = 9.78031 * (1.0 + 5.3e-3 * sd*sd);
    const double h45 = 1.00818e-2 * z + 2.465e-8 * z*z - 1.25e-13 * z*z*z + 2.8e-19 * z*z*z*z;
    const double k   = (g - 2e-5*z) / (9.80612 - 2e-5 * z);
    const double dp  = -8.5e-6 * z + 1.4e-9 * z*z;

    return h45 * k - (option ? dp : 0.0);
  }
}

#endif