#ifndef __JACOUSTICS__JACOUSTICSSUPPORTKIT__
#define __JACOUSTICS__JACOUSTICSSUPPORTKIT__

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


/**
 * \file
 *
 * Acoustics support kit.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  static const double  TOAMIN_S                  =    0.0000000;   //!< Minimal allowed time-of-arrival [s]
  static const double  TOAMAX_S                  =    0.6310912;   //!< Maximal allowed time-of-arrival [s]


  /**
   * Auxiliary data structure for calculation of attenuation length.
   *
   * Attenuation according to Ainslie and McColm, J. Acoust. Soc. AM. 103 (3) 1671 1998.\n
   * See also http://resource.npl.co.uk/acoustics/techguides/seaabsorption/physics.html for numerical values.
   */
  struct JAttenuationLength {
    /**
     * Constructor.
     *
     * \param  T                    temperature [C]
     * \param  S                    salinity [pm]
     * \param  pH                   pH
     */
    JAttenuationLength(const double T,
		       const double S,
		       const double pH) :
      T (T),
      S (S),
      pH(pH)
    {}


    /**
     * Get attentuation length for given frequency at given depth.
     *
     * \param  D_m        depth [m]
     * \param  f_kHz      frequency [kHz]
     * \return            attenuation length [m]
     */
    double operator()(const double D_m,
		      const double f_kHz) const
    {
      const double z  = D_m * 1.0e-3;
      const double f1 = 0.78    * sqrt(S/35.0) * exp(T/26.0);
      const double f2 = 42.0    * exp(T/17.0);
      const double a1 = 0.106   * f1*(f_kHz*f_kHz) / (f_kHz*f_kHz + f1*f1) * exp((pH-8.0)/0.56);
      const double a2 = 0.52    * (1.0 + T/43.0) * (S/35.0) * f2 * f_kHz*f_kHz / (f_kHz*f_kHz + f2*f2) * exp(-z/6.0);
      const double a0 = 0.00049 * exp(-T/27.0 + z/17.0) * f_kHz*f_kHz;
      const double A  = a0 + a1 + a2;  // dB/km
      
      return 10.0e3 / (A*log(10.0));
    }

    double T;
    double S;
    double pH;
  };


  /**
   * Function object to calculate attenutation length.
   */
  static const JAttenuationLength getAttenuationLength(13.2,      // temperature [deg]
						       38.0,      // salinity [pm]
						       8.0);      // pH




  /**
   * Get relative quality for given frequency at given distance.
   *
   * \param  D_m        depth     [m]
   * \param  f_kHz      frequency [kHz]
   * \param  d_m        distance  [m]
   * \return            quality
   */
  inline double getQ(const double D_m, 
		     const double f_kHz,
		     const double d_m)
  {
    return exp(-0.5 * d_m / getAttenuationLength(D_m, f_kHz)) / d_m;
  }


  /**
   * Utility class for emitter power and frequency.
   */
  struct JWaveform {
    /**
     * Default constructor.
     */
    JWaveform() :
      Q0   (0.0),
      f_kHz(0.0)
    {}


    /**
     * Constructor.
     *
     * Note that the power corresponds to the minimal distance (see method JWaveform::getDmin()).
     *
     * \param  Q0         power [quality]
     * \param  f_kHz      frequency [kHz]
     */
    JWaveform(const double Q0,
	      const double f_kHz) :
      Q0   (Q0),
      f_kHz(f_kHz)
    {}


    /**
     * Get minimal distance.
     *
     * \return            distance [m]
     */
    static double getDmin()
    {
      return 1.0;
    }


    /**
     * Get quality at given distance.
     *
     * \param  D_m        depth     [m]
     * \param  d_m        distance  [m]
     * \return            quality
     */
    double getQ(const double D_m,
		const double d_m) const
    {
      if (d_m > getDmin())
	return Q0 * JPP::getQ(D_m, this->f_kHz, d_m) / JPP::getQ(D_m, this->f_kHz, getDmin());
      else
	return Q0;
    }


    /**
     * Read waveform from input stream.
     *
     * \param  in         input stream
     * \param  waveform   waveform
     * \return            input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JWaveform& waveform)
    {
      return in >> waveform.Q0 >> waveform.f_kHz;
    }


    /**
     * Write waveform to output stream.
     *
     * \param  out        output stream
     * \param  waveform   waveform
     * \return            output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JWaveform& waveform)
    {
      return out << waveform.Q0 << ' ' << waveform.f_kHz;
    }

    double Q0;
    double f_kHz;
  };
}

#endif