#ifndef __JPHYSICS__JLED__
#define __JPHYSICS__JLED__

#include "JPhysics/JAbstractPMT.hh"
#include "JPhysics/JAbstractMedium.hh"
#include "JPhysics/JDispersionInterface.hh"
#include "JPhysics/JDispersion.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JElement.hh"
#include "JTools/JQuadrature.hh"


/**
 * \author mdejong
 */

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

namespace JPHYSICS {


  using JTOOLS::JQuadrature;
  using JTOOLS::JGaussLegendre;
  using JTOOLS::JCotangent;

  typedef JTOOLS::JElement2D<double, double> JElement2D_t;
  typedef JTOOLS::JElement3D<double, double> JElement3D_t;


  /**
   * Interface for emission profile from LED.
   */
  class JAbstractLED {
  public:

    /**
     * Virtual destructor.
     */
    virtual ~JAbstractLED()
    {}


    /**
     * Light yield from LED (number of p.e. per unit solid angle per unit time).
     *
     * \param  ct         zenith  angle of emission
     * \param  phi        azimuth angle of emission
     * \param  dt         time of emission [ns]
     * \return            d^2P / dOmega dt [npe/ns/sr]
     */
    virtual double getLightFromLED(const double ct, 
				   const double phi,
				   const double dt) const = 0;
  };


  /**
   * Probability Density Functions of the time response of a PMT.
   */
  class JLED : 
    public virtual JDispersionInterface, 
    public JAbstractPMT, 
    public JAbstractLED,
    public JAbstractMedium
  {
  public:
    /**
     * Constructor.
     *
     * \param  lambda               wavelength photon [nm]
     * \param  Tmin_ns              minimal time of emmision [ns]
     * \param  Tmax_ns              minimal time of emmision [ns]
     * \param  engine               scattering angle integrator
     * \param  numberOfPoints       number    of points for integration
     * \param  epsilon              precision of points for integration
     */
    JLED(const double lambda,
	 const double Tmin_ns,
	 const double Tmax_ns,
	 const JQuadrature& engine         = JCotangent(20),
	 const int          numberOfPoints = 20,
	 const double       epsilon        = 1e-12) :
      wavelength(lambda),
      tmin(Tmin_ns),
      tmax(Tmax_ns),
      main_engine(JGaussLegendre(numberOfPoints,epsilon)),
      beta_engine(engine)
    {
      const double dp = PI / numberOfPoints;

      for (double x = 0.5*dp; x < PI; x += dp) {
        phi_engine.push_back(JElement3D_t(cos(x),sin(x), dp));
      }
    }


    /**
     * Probability density function for direct light from LED.
     *
     * \param  D_m        distance between LED and PMT [m]
     * \param  cd         cosine angle LED orientation and LED - PMT position
     * \param  theta      zenith  angle orientation PMT
     * \param  phi        azimuth angle orientation PMT
     * \param  t_ns       time difference relative to direct light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getDirectLightFromLED(const double D_m,
				 const double cd,
				 const double theta,
				 const double phi,
				 const double t_ns) const
    {
      const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
      const double d  =  D_m;                                  // photon path [m]

      const double A  =  getPhotocathodeArea();

      const double px =  sin(theta)*cos(phi);
      //const double py =  sin(theta)*sin(phi);
      const double pz =  cos(theta);

      const double ct = sd*px + cd*pz;                         // cosine angle of incidence on PMT

      const double qe    = getQE(wavelength);
      const double l_abs = getAbsorptionLength(wavelength);
      const double ls    = getScatteringLength(wavelength);

      const double U  = getAngularAcceptance(ct);
      const double V  = exp(-d/l_abs) * exp(-d/ls);            // absorption & scattering
      const double W  = A/(d*d);                               // solid angle
  
      return qe * getLightFromLED(cd, 0.0, t_ns) * U * V * W;
    }

				 
    /**
     * Probability density function for scattered light from LED.
     *
     * \param  D_m        distance between LED and PMT [m]
     * \param  cd         cosine angle LED orientation and LED - PMT position
     * \param  theta      zenith  angle orientation PMT
     * \param  phi        azimuth angle orientation PMT
     * \param  t_ns       time difference relative to direct light [ns]
     * \return            dP/dt [npe/ns]
     */
    double getScatteredLightFromLED(const double D_m,
				    const double cd,
				    const double theta,
				    const double phi,
				    const double t_ns) const
    {
      using namespace std;

      double value = 0;

      if (t_ns >= tmin) {

	const double ng =  getIndexOfRefractionGroup(wavelength);

	const double sd =  sqrt((1.0 + cd)*(1.0 - cd));
	const double l  =  D_m;                                    // distance [m]

	const double A  =  getPhotocathodeArea();
	//const double sr =  sqrt(A/PI) / D_m;
	//const double cr =  sqrt((1.0 + sr)*(1.0 - sr));

	const double px =  sin(theta)*cos(phi);
	const double py =  sin(theta)*sin(phi);
	const double pz =  cos(theta);

	const double qx = cd*px +  0  - sd*pz;
	const double qy =         1*py;
	const double qz = sd*px +  0  + cd*pz;

	const double qe    = getQE(wavelength);
	const double l_abs = getAbsorptionLength(wavelength);
	const double ls    = getScatteringLength(wavelength);

	const double Jc = 1.0 / ls;                                // dN/dx

	const double xmin = tmin;
	const double xmax = std::min(tmax, t_ns);

	JQuadrature buffer;
	
	if (xmax > xmin) {

	  for (JQuadrature::const_iterator i = main_engine.begin(); i != main_engine.end(); ++i) {

	    const double t  = 0.5 * (xmax + xmin)  +  i->getX() * 0.5 * (xmax - xmin);
	    const double dt = i->getY() * 0.5 * (xmax - xmin);

	    buffer[t] = dt;
	  }

	} else {

	  buffer[tmin] = 1.0;
	}

	for (JQuadrature::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {

	  const double t  = i->getX();
	  const double dt = i->getY();

	  const double d  = l + C*(t_ns - t)/ng;                   // photon path
	  const double V  = exp(-d/l_abs);                         // absorption

	  for (JQuadrature::const_iterator j = beta_engine.begin(); j != beta_engine.end(); ++j) {

	    const double cb = j->getX();
	    const double dc = j->getY();
	    
	    const double sb = sqrt((1.0 + cb)*(1.0 - cb));
	    
	    const double v  = 0.5 * (d + l) * (d - l) / (d - l*cb);
	    const double u  = d - v;
	    
	    if (u <= 0) continue;
	    if (v <= 0) continue;
	    
	    const double cts = (l*cb - v) / u;                   // cosine scattering angle
	    
	    const double W  = min(A/(v*v), 2.0*PI);              // solid angle
	    const double Ja = getScatteringProbability(cts);     // P(cos(),phi)
	    const double Jd = ng * (1.0 - cts) / C;              // dt/du
	    
	    const double ca = (l - v*cb) / u;
	    const double sa =      v*sb  / u;
	    
	    for (std::vector<JElement3D_t>::const_iterator k = phi_engine.begin(); k != phi_engine.end(); ++k) {
	      
	      const double cp = k->getX();
	      const double sp = k->getY();
	      const double dp = k->getZ();

	      const double dom = dc * dp * v*v / (u*u);
	      
	      const double ct0 = cd*ca - sd*sa*cp;
	      const double ph0 = atan2(sa*sp, cd*sa*cp + sd*ca);
	      
	      const double vx  = -sb*cp * qx;
	      const double vy  = -sb*sp * qy;
	      const double vz  =  cb    * qz;
	      
	      const double U  =
		getAngularAcceptance(vx + vy + vz) * getLightFromLED(ct0, +ph0, t) +
		getAngularAcceptance(vx - vy + vz) * getLightFromLED(ct0, -ph0, t);
	      
	      value += qe * dt * dom * Ja * Jc * U * V * W / Jd;
	    }
	  }
	}
      }

      return value; 
    }				 

  protected:
    double wavelength;                    // wavelength photon [nm]

    double tmin;                          // minimal time of emmision [ns]
    double tmax;                          // maximal time of emmision [ns]

    JQuadrature main_engine;              // numerical integration
    JQuadrature beta_engine;              // numerical integration of scattering angle
    std::vector<JElement3D_t> phi_engine; // numerical integration of phi angle
  };


  /**
   * Probability Density Functions of the time response of a PMT (C-like interface)
   */
  class JLED_C : 
    public JLED,
    public JDispersion 
  {
  public:
    /**
     * Constructor.
     *
     * \param  Area            photo-cathode area [m^2]
     * \param  LED             pointer to interface for emission profile from LED
     * \param  QE              pointer to function for quantum efficiency of PMT
     * \param  Pmt             pointer to function for angular acceptance of PMT
     * \param  L_abs           absorption length [m]
     * \param  L_s             scattering length [m]
     * \param  Km3             pointer to model specific function to describe light scattering in water
     * \param  P_atm           ambient pressure [atm]
     * \param  wavelength      wavelength of light [nm]
     * \param  Tmin_ns         minimal time of emmision [ns]
     * \param  Tmax_ns         minimal time of emmision [ns]
     * \param  engine          scattering angle integrator
     * \param  numberOfPoints  number    of points for integration
     * \param  epsilon         precision of points for integration
     */
    JLED_C(const double Area,
           //
           // pointers to static functions
           //
	   const JAbstractLED* LED,
           double (*QE)  (const double),
           double (*Pmt) (const double),
           const double L_abs,
           const double L_s,
           double (*Km3) (const double),
           //
           // parameters for base class
           //
	   const double P_atm,
	   const double wavelength,
	   const double Tmin_ns,
	   const double Tmax_ns,
	   const JQuadrature& engine   = JCotangent(20),
           const int    numberOfPoints = 20,
           const double epsilon        = 1e-12) :
      JLED(wavelength, Tmin_ns, Tmax_ns, engine, numberOfPoints, epsilon),
      JDispersion(P_atm),
      A    (Area ),
      led  (LED  ),
      qe   (QE   ),
      l_abs(L_abs),
      ls   (L_s  ),
      pmt  (Pmt  ),
      km3  (Km3  )
    {}


    /**
     * Light yield from LED (number of p.e. per unit solid angle per unit time).
     *
     * \param  ct         zenith  angle of emission
     * \param  phi        azimuth angle of emission
     * \param  dt         time of emission [ns]
     * \return            d^2P / dOmega dt [npe/ns/sr]
     */
    virtual double getLightFromLED(const double ct, 
				   const double phi, 
				   const double dt) const
    { 
      return led->getLightFromLED(ct, phi, dt);
    }


    /**
     * Photo-cathode area of PMT.
     *
     * \return            photo-cathode area [m^2]
     */
    virtual double getPhotocathodeArea() const override 
    { 
      return A;
    }


    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.).
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE
     */
    virtual double getQE(const double lambda) const override 
    { 
      return qe(lambda);
    }

 
    /**
     * Angular acceptance of PMT.
     *
     * \param  ct         cosine angle of incidence
     * \return            relative efficiency
     */
    virtual double getAngularAcceptance(const double ct) const override 
    { 
      return pmt(ct);
    }


    /**
     * Absorption length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            absorption length [m]
     */
    virtual double getAbsorptionLength(const double lambda) const override 
    { 
      return l_abs;
    }


    /**
     * Scattering length.
     *
     * \param  lambda     wavelenth [nm]
     * \return            scattering length [m]
     */
    virtual double getScatteringLength(const double lambda) const override 
    { 
      return ls;
    }


    /**
     * Model specific function to describe light scattering in water
     * (integral over full solid angle normalised to one).
     *
     * \param  ct         cosine scattering angle
     * \return            probability
     */
    virtual double getScatteringProbability(const double ct) const override 
    { 
      return km3(ct);
    }

  protected:
    /**
     * photo-cathode area [m2]
     */
    const double A;

    /**
     * Pointer to interface for emission profile from LED
     */
    const JAbstractLED* led;

    /**
     * Quantum efficiency of PMT (incl. absorption in glass, gel, etc.)
     *
     * \param  lambda     wavelenth [nm]
     * \return            QE
     */
    double (*qe)(const double lambda);

    /**
     * Absorption length
     */
    double l_abs;

    /**
     * Scattering length
     */
    double ls;

    /**
     * Angular acceptance of PMT
     *
     * \param  ct     cosine angle of incidence
     * \return        relative efficiency
     */
    double (*pmt)(const double ct);

    /**
     * Model specific function to describe light scattering in water
     *
     * \param  ct     cosine scattering angle
     * \return        probability
     */
    double (*km3)(const double ct);
   }; 
}

#endif