JQuadrature.hh

Go to the documentation of this file.

#ifndef __JTOOLS__JQUADRATURE__
#define __JTOOLS__JQUADRATURE__

#include <cmath>

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


/**
 * \file
 *
 * Auxiliary classes for numerical integration.
 * \author mdejong
 */
namespace JTOOLS {}
namespace JPP { using namespace JTOOLS; }

namespace JTOOLS {

  /**
   * Type definition of basic element for quadratures.
   */
  typedef JElement2D<double, double>               JElement2D_t;
  

  /**
   * Type definition for numerical integration.
   *
   *          \f$\displaystyle \int_{x_1}^{x_2} f(x) dx = \sum_{i=1}^{N} w_i f(x_i) \f$
   *
   * The abscissa and ordinate values of the collection can be used
   * as abscissa and weight values of the summation to approximately
   * determine the integral of the function.
   */
  class JQuadrature :
    public JCollection<JElement2D_t> 
  {
  public:
    /**
     * Default constructor.
     */
    JQuadrature()
    {}


    /**
     * General purpose constructor.
     *
     * The template argument should correspond to a function requiring two arguments.
     * These two arguments should correspond to the lower and upper integration limit, respectively.
     * The given function should return the value of the integral between the two integration limits. 
     *
     * \param  xmin       minimal x
     * \param  xmax       maximal x
     * \param  nx         number of points
     * \param  integral   integral
     * \param  eps        precision
     */
    template<class JFunction_t>
    JQuadrature(const double xmin,
                const double xmax,
                const int    nx,
                JFunction_t  integral,
                const double eps = 1.0e-4)
    {
      double Xmin = xmin;
      double Xmax = xmax;

      const double Vmin = integral(Xmin, Xmax) / (double) nx;
    
      for (int i = 0; i != nx; ++i) {
      
        for (double xmin = Xmin, xmax = Xmax; ; ) {
          
          const double x = 0.5 * (xmin + xmax);
          const double v = integral(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 / integral(__x);
          
            insert(JElement2D_t(__x,__y));
          
            Xmin = x;
            xmax = Xmax;
            
            break;
          }
          
          if (v < Vmin)
            xmin = x;
          else
            xmax = x;
        }
      }
    }
  };


  /**
   * Numerical integrator for <tt>W(x) = 1</tt>.
   * 
   * 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 <tt>W(x) = x^a e^(-x)</tt>.
   *
   * 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 = -tgamma(alf+n) / tgamma((double) n) / (pp*n*p1);
        
        insert(JElement2D_t(z,y));
      }
    }
  };


  /**
   * Numerical integrator for <tt>W(x) = e^-(x^2)</tt>.
   *
   * 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 <tt>W(x) = (1 + g^2 - 2gx)^a</tt>.
   * For this, <tt>g > 0</tt>.
   *
   * 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;
        
        insert(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;
        
        insert(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;
        
        insert(JElement2D_t(x,dx));
      }
    }
  };


  /**
   * Numerical integrator for <tt>W(x) = (1 + g*x*x)</tt>.
   * For this, <tt>g > 0</tt>.
   *
   * 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);
        
        insert(JElement2D_t(x, dx*gi*dy)); 
      }
    }
  };


  /**
   * Numerical integrator for <tt>W(x) = |x| / sqrt(1 - x*x)</tt>
   *
   * 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(JElement2D_t(+cb, dc));   
        insert(JElement2D_t(-cb, dc));   
      }
    }
  };


  /**
   * Numerical integrator for <tt>W(x) = |x| / sqrt(1 - x*x), x > 0 </tt> and <tt>W(x) =  1, x <= 0</tt>.
   *
   * 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(JElement2D_t(cb, dc));   
      }

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

#endif