#ifndef __JQUADRATURE__
#define __JQUADRATURE__

#include <vector>
#include <cmath>

#include "TMath.h"

#include "JTools/JConstants.hh"
#include "JTools/JElement.hh"


namespace JTOOLS {

  namespace {
    using JTOOLS::PI;
    using JTOOLS::JElement2D_t;
  }


  /**
   * Type definition for numerical integration.
   *
   *          
   *           _x_2                 N-1
   *          |                     --
   *          | dx W(x) f(x)   =    \   w_i f(x_i)
   *          |                     /_ 
   *         _|                      i=0
   *      x_1            
   */
  class JQuadrature :
    public std::vector<JElement2D_t> 
  {
  public:
    /**
     * Default constructor.
     */
    JQuadrature() :
      std::vector<JElement2D_t>()
    {}


    /**
     * Constructor.
     * The template argument corresponds to a 2D function.
     * This function should return the value of the integral between the two integration limits.
     *
     * \param  xmin   minimal x
     * \param  xmax   maximal x
     * \param  n      number of points
     * \param  fcn    function(xmin, xmax)
     * \param  eps    precision
     */
    template<class T>
    JQuadrature(const double xmin,
		const double xmax,
		const int    n,
		const T&     fcn,
		const double eps = 1.0e-4) :
      std::vector<JElement2D_t>()
    {
      double Xmin = xmin;
      double Xmax = xmax;

      const double Vmin = fcn(Xmin, Xmax) / (double) n;
    
      for (int i = 0; i != n; ++i) {
      
	for (double xmin = Xmin, xmax = Xmax; ; ) {
	  
	  const double x = 0.5 * (xmin + xmax);
	  const double v = fcn(Xmin, x);
	
	  if (fabs(Vmin - v)    < eps *  Vmin ||
	      fabs(xmax - xmin) < eps * (Xmax - Xmin)) {
	    
	    const double __x = 0.5 * (Xmin + x);
	    const double __y = Vmin / fcn(__x);
	  
	    push_back(JElement2D_t(__x,__y));
	  
	    Xmin = x;
	    xmax = Xmax;
	    
	    break;
	  }
	  
	  if (v < Vmin)
	    xmin = x;
	  else
	    xmax = x;
	}
      }
    }
  };


  /**
   * Numerical integrator for W(x) = 1.
   * 
   * Gauss-Legendre integration code is taken from reference:
   * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, 
   * Cambridge University Press.
   */ 
  class JGaussLegendre : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  eps    precision
     */
    JGaussLegendre(const int    n,
		   const double eps = 1.0e-12) :
      JQuadrature()
    {
      resize(n);

      const int M = (n + 1) / 2;
      
      double p0, p1, p2, pp;
      
      for (int i = 0; i < M; ++i) {
	
	double z  = cos(PI * (i+0.75) / (n+0.5));
	double z1;
	
	do {
	  
	  p1 = 0.0;
	  p2 = 1.0;
	  
	  // recurrence relation
	  
	  for (int j = 0; j < n; ++j) {
	    p0 = p1;
	    p1 = p2;
	    p2 = ((2*j + 1) * z*p1  - j*p0) / (j+1);
	  }
	  
	  pp = n * (z*p2 - p1) / (z*z - 1.0);
	  
	  z1 = z;
	  z  = z1 - p2/pp;
	  
	} while (fabs(z-z1) > eps);

	const double y = 2.0 / ((1.0-z*z)*pp*pp);

	at(  i  ) = JElement2D_t(-z,y);
	at(n-i-1) = JElement2D_t(+z,y);
      }
    }
  };


  /**
   * Numerical integrator for W(x) = x^a e^(-x).
   *
   * Gauss-Laguerre integration code is taken from reference:
   * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, 
   * Cambridge University Press.
   */ 
  class JGaussLaguerre : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  alf    power
     * \param  eps    precision
     */
    JGaussLaguerre(const int    n,
		   const double alf,
		   const double eps = 1.0e-12) :
      JQuadrature()
    {
      const int number_of_iterations = 100;

      double z, z1;
      double p0, p1, p2, pp;
      
      for (int i = 0; i < n; ++i) {
	
	switch (i) {

	case 0:
	  z  = (1.0 + alf) * (3.0 + 0.92*alf) / (1.0 + 2.4*n + 1.8*alf);
	  break;

	case 1:
	  z += (15.0 + 6.25*alf) / (1.0 + 0.9*alf + 2.5*n); 
	  break;

	default:
	  const double ai = i - 1;
	  z += ((1.0+2.55*ai)/(1.9*ai) + (1.26*ai*alf)/(1.0+3.5*ai)) * (z - at(i-2).getX()) / (1.0 + 0.3*alf);
	  break;
	}


	int k;

	for (k = 0; k != number_of_iterations; ++k) {
	  
	  p1 = 0.0;
	  p2 = 1.0;
	  
	  // recurrence relation
	  
	  for (int j = 0; j < n; ++j) {
	    p0 = p1;
	    p1 = p2;
	    p2 = ((2*j + 1 + alf - z) * p1  - (j + alf)*p0) / (j+1);
	  }
	  
	  pp = (n*p2 - (n+alf)*p1) / z;
	  
	  z1 = z;
	  z  = z1 - p2/pp;

	  if (fabs(z-z1) < eps)
	    break;
	}

	const double y = -TMath::Gamma(alf+n) / TMath::Gamma((double) n) / (pp*n*p1);
	
	push_back(JElement2D_t(z,y));
      }
    }
  };


  /**
   * Numerical integrator for W(x) = e^-(x^2).
   *
   * Gauss-Hermite integration code is taken from reference:
   * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, 
   * Cambridge University Press.
   */ 
  class JGaussHermite : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  eps    precision
     */
    JGaussHermite(const int    n,
		  const double eps = 1.0e-12) :
      JQuadrature()
    {
      resize(n);

      const double pii = 1.0 / pow(PI,0.25);

      const int number_of_iterations = 100;

      const int M = (n + 1) / 2;
      
      double p0, p1, p2, pp;
      double z = 0.0;
      double z1;
      
      for (int i = 0; i < M; ++i) {

	switch (i) {

	case 0:
	  z  = sqrt((double) (2*n+1))  -  1.85575 * pow((double) (2*n+1),-0.16667);
	  break;

	case 1:
	  z -= 1.14 * pow((double) n,0.426) / z;
	  break;

	case 2:
	  z  = 1.86*z + 0.86*at( 0 ).getX();
	  break;

	case 3:
	  z  = 1.91*z + 0.91*at( 1 ).getX();
	  break;

	default:
	  z  = 2.00*z + 1.00*at(i-2).getX();
	  break;
	}

	for (int k = 0; k != number_of_iterations; ++k) {
	  
	  p1 = 0.0;
	  p2 = pii;
	  
	  // recurrence relation
	  
	  for (int j = 0; j < n; ++j) {
	    p0 = p1;
	    p1 = p2;
	    p2 = z * sqrt(2.0/(double) (j+1)) * p1  -  sqrt((double) j / (double) (j+1)) * p0;
	  }
	  
	  pp = sqrt((double) (2*n)) * p1;
	  
	  z1 = z;
	  z  = z1 - p2/pp;

	  if (fabs(z-z1) < eps)
	    break;
	} 

	const double y = 2.0 / (pp*pp);

	at(  i  ) = JElement2D_t(-z,y);
	at(n-i-1) = JElement2D_t(+z,y);
      }
    }
  };


  /**
   * Numerical integrator for W(x) = (1 + g^2 - 2gx)^a.
   * For this, g > 0.
   *
   * Henyey-Greenstein integration points and weights.
   */ 
  class JHenyeyGreenstein : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  g      angular dependence parameter
     * \param  a      power
     */
    JHenyeyGreenstein(const int    n,
		      const double g,
		      const double a) :
      JQuadrature()
    {
      const double b  = -2*g * (a + 1.0);
      const double ai =  1.0 / (a + 1.0);

      const double ymin = pow(1.0 + g, 2*(a + 1.0)) / b;
      const double ymax = pow(1.0 - g, 2*(a + 1.0)) / b;

      const double dy = (ymax - ymin) / (n + 1);
      
      for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
	
	const double v  = y*b;
	const double w  = pow(v, ai);
	const double x  = (1.0 + g*g - w) / (2*g);
	const double dx = pow(v, -a*ai)*dy;
	
	push_back(JElement2D_t(x,dx));
      }
    }


    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  g      angular dependence parameter
     * \param  a      power
     * \param  xmin   minimal value
     * \param  xmax   maximal value
     */
    JHenyeyGreenstein(const int    n,
		      const double g,
		      const double a,
		      const double xmin,
		      const double xmax) :
      JQuadrature()
    {
      const double b  = -2*g * (a + 1.0);
      const double ai =  1.0 / (a + 1.0);

      const double ymin = pow(1.0 + g*g -2*g*xmin, a + 1.0) / b;
      const double ymax = pow(1.0 + g*g -2*g*xmax, a + 1.0) / b;

      const double dy = (ymax - ymin) / (n + 1);
      
      for (double y = ymax - 0.5*dy; y > ymin; y -= dy) {
	
	const double v  = y*b;
	const double w  = pow(v, ai);
	const double x  = (1.0 + g*g - w) / (2*g);
	const double dx = pow(v, -a*ai)*dy;
	
	push_back(JElement2D_t(x,dx));
      }
    }


    /**
     * Constructor for special case where a = -1.
     *
     * \param  n      number of points
     * \param  g      angular dependence parameter
     */
    JHenyeyGreenstein(const int    n,
		      const double g) :
      JQuadrature()
    {
      const double dy = 1.0 / (n + 1);
      const double gi = log((1.0 + g*g) / (1.0 - g*g)) / (2.0*g);
      
      for (double y = 1.0 - 0.5*dy; y > 0.0; y -= dy) {
	
	const double v  = -y*2.0*g*gi  +  log(1.0 + g*g);
	const double w  = exp(v);
	const double x  = (1.0 + g*g - w) / (2.0*g);
	const double dx = w*gi*dy;
	
	push_back(JElement2D_t(x,dx));
      }
    }
  };


  /**
   * Numerical integrator for W(x) = (1 + g*x*x).
   * For this, g > 0.
   *
   * Rayleigh integration points and weights.
   */ 
  class JRayleigh : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     * \param  g      angular dependence parameter
     */
    JRayleigh(const int    n,
	      const double g) :
      JQuadrature()
    {
      const double dy = 1.0 / (n + 1);
      const double gi = 3.0/g  +  1.0;
      
      // t^3 + 3pt + 2q = 0
      
      const double p  = 1.0/g;
      
      for (double y = 0.5*dy; y < 1.0; y += dy) {
	
	const double q  = 0.5*gi - gi*y;
	
	const double b  = sqrt(q*q + p*p*p);
	const double u  = pow(-q + b, 1.0/3.0);
	const double v  = pow(+q + b, 1.0/3.0);
	
	const double x  =  u - v;
	const double dx = (u + v) / (3.0*b);
	
	push_back(JElement2D_t(x, dx*gi*dy)); 
      }
    }
  };


  /**                                       _______
   * Numerical integrator for W(x) = |x| / V1 - x*x
   *
   * Co-tangent integration points and weights.
   */ 
  class JCotangent : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     */
    JCotangent(const int n) :
      JQuadrature()
    {
      for (double ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {

	const double cb = sqrt((1.0 + sb)*(1.0 - sb));
	const double dc = ds*sb/cb;
	  
	insert(end(),   JElement2D_t(+cb, dc));   
	insert(begin(), JElement2D_t(-cb, dc));   
      }
    }
  };


  /**                                       _______
   * Numerical integrator for W(x) = |x| / V1 - x*x ,  x > 0
                                   =  1             ,  x < 0
   *
   * Bi-tangent integration points and weights.
   */ 
  class JBitangent : 
    public JQuadrature 
  {
  public:
    /**
     * Constructor.
     *
     * \param  n      number of points
     */
    JBitangent(const int n) :
      JQuadrature()
    {
      double sb, ds;
      double cb = 0.0;
      double dc = 0.0;
	  
      for (ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {

	cb = sqrt((1.0 + sb)*(1.0 - sb));
	dc = ds*sb/cb;
	  
	insert(begin(), JElement2D_t(cb, dc));   
      }

      for (dc = (cb + 1.0) / (n/2), cb -= 0.5*dc ; cb > -1.0; cb -= dc)
	insert(begin(), JElement2D_t(cb, dc));   
    }
  };
}

#endif