#ifndef __JPHYSICS__JDISPERSION__
#define __JPHYSICS__JDISPERSION__

#include "JPhysics/JDispersionInterface.hh"


/**
 * \author mdejong
 */

namespace JPHYSICS {}
namespace JPP { using namespace JPHYSICS; }

namespace JPHYSICS {

  /**
   * Implementation of dispersion for water in deep sea.
   * This class implements the JDispersionInterface interface.
   *
   * Light dispersion data are taken from reference:
   * David J.L. Bailey, 
   * "Monte Carlo tools and analysis methods for understanding the ANTARES experiment and 
   * predicting its sensitivity to Dark Matter", 
   * PhD thesis,  University of Oxford, United Kingdom, 2002.
   */
  class JDispersion : 
    public virtual JDispersionInterface
  {
  public:
    /**
     * Constructor.
     *
     * \param  P_atm           ambient pressure [atm]
     */
    JDispersion(const double P_atm) :
      P (P_atm),              // ambient pressure [atm]
      a0(  1.3201),           // offset
      a1(  1.4e-5),           // dn/dP
      a2( 16.2566),           // d^1n/(dx)^1
      a3(-4383.0),            // d^2n/(dx)^2
      a4(  1.1455e6)          // d^3n/(dx)^3
    {}


    /**
     * Index of refraction (phase velocity).
     *
     * \param  lambda     wavelenth [nm]
     * \return            index of refraction
     */
    virtual double getIndexOfRefractionPhase(const double lambda) const
    {
      const double x = 1.0 / lambda;

      return a0  +  a1*P  +  x*(a2 + x*(a3 + x*a4));
    }

   
    /**
     * Dispersion of light for phase velocity.
     *
     * \param  lambda     wavelength of light [nm]
     * \return            dn/dlambda
     */
    virtual double getDispersionPhase(const double lambda) const
    {
      const double x = 1.0 / lambda;
      
      return -x*x*(a2 + x*(2.0*a3 + x*3.0*a4));
    }


    /**
     * Dispersion of light for group velocity.
     *
     * \param  lambda     wavelength of light [nm]
     * \return            dn/dlambda
     */
    virtual double getDispersionGroup(const double lambda) const
    {
      const double x = 1.0 / lambda;

      const double n   = getIndexOfRefractionPhase(lambda);
      const double np  = getDispersionPhase(lambda);
      const double npp = x*x*x*(2.0*a2 + x*(6.0*a3 + x*12.0*a4));
      const double ng  = n / (1.0 + np*lambda/n);
      
      return ng*ng * (2*np*np - n*npp) * lambda / (n*n*n);
    }


    /**
     * Dispersion parameters (x = 1/lambda)
     */
    const double P;            //!<  ambient pressure [atm] 
    const double a0;           //!<  offset
    const double a1;           //!<  dn/dP
    const double a2;           //!<  d^1n/(dx)^1
    const double a3;           //!<  d^2n/(dx)^2
    const double a4;           //!<  d^3n/(dx)^3    
  };
}

#endif