#ifndef __JACOUSTICS__JSOUNDVELOCITY__
#define __JACOUSTICS__JSOUNDVELOCITY__

#include <istream>
#include <ostream>
#include <cmath>

#include "JSystem/JStat.hh"
#include "JLang/JStringStream.hh"

#include "JAcoustics/JAbstractSoundVelocity.hh"


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

namespace JACOUSTICS {

  /**
   * Implementation for depth dependend velocity of sound.
   *
   * Note that
   *  - z-axis is pointing upwards; and
   *  - velocity is defined at a given depth.
   *
   * It is assumed that the dependence of the sound velocity on the depth is non-zero.\n
   * The depth is commonly referred to as the seabed.
   */
  struct JSoundVelocity :
    public JAbstractSoundVelocity
  {
    /**
     * Constructor.
     *
     * Note that the depth is relative to the sea surface and therefore should be negative.
     *
     * \param  a              velocity  [m/s]
     * \param  b              d(v)/d(z) [m/s/m]
     * \param  z0             depth     [m]
     */
    JSoundVelocity(const double a, const double b, const double z0) :
      a (a),
      b (b),
      z0(z0)
    {}


    /**
     * Set depth.
     *
     * Note that the depth is relative to the sea surface and therefore should be negative.
     *
     * \param  z0             depth     [m]
     * \return                this sound velocity
     */
    JSoundVelocity& set(const double z0)
    {
      this->a  += this->b * (z0 - this->z0);
      this->z0  = z0;
      
      return *this;
    }


    /**
     * Get sound velocity at given depth.
     *
     * Note that the depth is relative to the sea surface and therefore should be negative.
     *
     * \param  z0             depth     [m]
     * \return                sound velocity
     */
    JSoundVelocity operator[](const double z0) const
    {
      return JSoundVelocity(*this).set(z0);
    }


    /**
     * Get velocity of sound at seabed.
     *
     * \return                velocity [m/s]
     */
    virtual double operator()() const override
    {
      return a;
    }


    /**
     * Get velocity of sound at given depth relative to seabed.
     *
     * \param  z              depth [m]
     * \return                velocity [m/s]
     */
    virtual double operator()(const double z) const override
    {
      return a + b * z;
    }


    /**
     * Get distance travelled by sound.
     *
     * The distance travelled is obtained by consistency with method JSoundVelocity::getTime.\n
     * Note that the depth is relative to the seabed.
     *
     * \param  t_s            time  [s]
     * \param  z1             depth [m]
     * \param  z2             depth [m]
     * \return                distance [m]
     */
    virtual double getDistance(const double t_s,
			       const double z1,
			       const double z2) const override
    {
      const double eps = 1.0e-8;

      const double v1  = (*this)(z1);
      const double v2  = (*this)(z2);
      const double lv  = log(v2/v1);

      if (fabs(lv) > eps)
	return t_s * fabs((z2 - z1) * b / lv);
      else
	return 0.5 * t_s * (v1 + v2);
    }


    /**
     * Get propagation time of sound.
     *
     * The propagation time is obtained from the integral \f$\Delta T = \frac{D}{\Delta z} \int\frac{1}{V(z)}dz\f$.\n
     * Note that the depth is relative to the seabed.
     *
     * \param  D_m            distance [m]
     * \param  z1             depth [m]
     * \param  z2             depth [m]
     * \return                time  [s]
     */
    virtual double getTime(const double D_m,
			   const double z1,
			   const double z2) const override
    {
      const double eps = 1.0e-8;

      const double ct  = (z2 - z1) / D_m;
      const double v1  = (*this)(z1);
      const double v2  = (*this)(z2);

      if (fabs(ct) > eps)
	return fabs(log(v2/v1) / (ct*b));
      else
	return 2.0 * D_m / (v1 + v2);
    }


    /**
     * Get inverse velocity of sound.
     *
     * The inverse velocity is obtained from the propagation time given by method JSoundVelocity::getTime.\n
     * Note that the depth is relative to the seabed.
     *
     * \param  D_m            distance [m]
     * \param  z1             depth [m]
     * \param  z2             depth [m]
     * \return                inverse velocity [s/m]
     */
    virtual double getInverseVelocity(const double D_m,
				      const double z1,
				      const double z2) const override
    {
      const double eps = 1.0e-8;

      const double ct  = (z2 - z1) / D_m;
      const double v1  = (*this)(z1);
      const double v2  = (*this)(z2);

      if (fabs(ct) > eps)
	return fabs(log(v2/v1) / ((z2 - z1) * b));
      else
	return 2.0 / (v1 + v2);
    }


    /**
     * Read sound velocity from input.
     *
     * \param  in            input stream
     * \param  velocity      sound velocity
     * \return               input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JSoundVelocity& velocity)
    {
      using namespace JPP;

      JStringStream is(in);

      if (getFileStatus(is.str().c_str())) {
        is.load();
      }

      if (is >> velocity.a) {

	is >> velocity.b
	   >> velocity.z0;

	is.clear();
      }

      return in;
    }


    /**
     * Write sound velocity to output.
     *
     * \param  out           output stream
     * \param  velocity      sound velocity
     * \return               output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JSoundVelocity& velocity)
    {
      out << velocity.a << ' ' << velocity.b << ' ' << velocity.z0;

      return out;
    }

  private:
    double a;
    double b;
    double z0;
  };


  /**
   * Function object for velocity of sound.
   */
  static const JSoundVelocity  getSoundVelocity(1541.0, -17.0e-3, -2000.0);
}

#endif