JQuadrature.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
 
#include "JTools/JQuadrature.hh"
 
#include "JMath/JPolynome.hh"
#include "JMath/JGauss.hh"
 
#include "JPhysics/JPhysicsSupportkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test JTOOLS::JQuadrature.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  double    precision;
  int       debug;
 
  try {
 
    JParser<> zap("Example program to test JTOOLS::JQuadrature.");
 
    zap['e'] = make_field(precision)      = 1.0e-4;
    zap['d'] = make_field(debug)          = 3;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  {
    const int    numberOfPoints = 10;
    const double epsilon        = 1.0e-12;
 
    const double xmin = -5.0;
    const double xmax = +3.0;
 
    /**
     * Function.
     */
    const JPolynome f1(1.0, 2.0, 3.0);
 
    /**
     * Integral.
     */
    const JPolynome F1 = f1.getIntegral();
          
    double V = F1(xmax) - F1(xmin);
 
    JGaussLegendre engine(numberOfPoints, epsilon);
 
    const double x0 = 0.5 * (xmax + xmin);
    const double x1 = 0.5 * (xmax - xmin);
 
    double W = 0.0;
 
    for (JQuadrature::const_iterator i = engine.begin(); i != engine.end(); ++i) {
 
      const double x = x0 + i->getX() * x1;
     
      W += f1(x) * i->getY() * x1;
    }
 
    DEBUG("integral analytical " << FIXED(9,5) << V << endl);
    DEBUG("integral numerical  " << FIXED(9,5) << W << endl);
 
    ASSERT(fabs(V - W) <= precision, "Test Gauss-Legendre integration.");
  }
  {
    const int    numberOfPoints = 10;
    const double a              = 1.0;
    const double epsilon        = 1.0e-12;
 
    const double xmin =  0.0;
    const double xmax = +1.0e10;  // approximation of infinity
 
    /**
     * Function.
     */
    const struct {
      inline double operator()(const double x) const
      {
        return sin(x);
      }
    } f1;
 
    /**
     * Integral.
     */
    const struct {
      inline double operator()(const double x) const
      {
        // integral of sin(x) * x^1 * exp(-x)
 
        return -0.5 * exp(-x) * (x*sin(x) + (x + 1.0)*cos(x));
      }
    } F1;
 
    const double V = F1(xmax) - F1(xmin);
 
    JGaussLaguerre engine(numberOfPoints, a, epsilon);
 
    double W = 0.0;
 
    for (JQuadrature::const_iterator i = engine.begin(); i != engine.end(); ++i) {
 
      const double x = i->getX();
     
      W += f1(x) * i->getY();
    }
 
    DEBUG("integral analytical " << FIXED(9,5) << V << endl);
    DEBUG("integral numerical  " << FIXED(9,5) << W << endl);
 
    ASSERT(fabs(V - W) <= precision, "Test Gauss-Laguerre integration.");
  }
  {
    const int    numberOfPoints = 10;
    const double epsilon        = 1.0e-12;
 
    const double sigma[] = { 4.0, 0.5 };
  
    const double xmin = -3.0 * sigma[0];
    const double xmax = +3.0 * sigma[0];
 
    /**
     * Function.
     */
    const JGauss f1(0.0, sigma[0]);
 
    /**
     * Integral.
     */
    const JGauss F1(0.0, hypot(sigma[0], sigma[1]));
 
    JGaussHermite engine(numberOfPoints, epsilon);
 
    for (double x = xmin; x <= xmax; x += 0.05*(xmax - xmin)) {
 
      const double V = F1(x);
      
      double W = 0.0;
 
      for (JQuadrature::const_iterator i = engine.begin(); i != engine.end(); ++i) {
 
        const double u = i->getX() * sigma[1] * sqrt(2.0);
        const double v = i->getY() / sqrt(PI);
     
        W += f1(x + u) * v;
      }
 
      DEBUG(FIXED(7,3) << x << ' ' << FIXED(9,5) << V << ' ' << FIXED(9,5) << W << endl);
 
      ASSERT(fabs(V - W) <= precision * V, "Test Gauss-Hermite integration.");
    }
  }
  {
    const int    numberOfPoints = 10;
    
    const double a = -1.5;
    const double g = 0.924;
 
    const double V = 1.0;
 
    JHenyeyGreenstein engine(numberOfPoints, g, a);
 
    double W = 0.0;
 
    for (JQuadrature::const_iterator i = engine.begin(); i != engine.end(); ++i) {
 
      const double x = i->getX();
     
      W += henyey_greenstein(g,x) * i->getY() * 2*PI;
    }
 
    DEBUG("integral analytical " << FIXED(9,5) << V << endl);
    DEBUG("integral numerical  " << FIXED(9,5) << W << endl);
 
    ASSERT(fabs(V - W) <= precision, "Test Henyey-Greenstein integration.");
  }
  {
    const int    numberOfPoints = 10;
    
    const double a = 0.835;
 
    const double V = 1.0;
 
    JRayleigh engine(numberOfPoints, a);
 
    double W = 0.0;
 
    for (JQuadrature::const_iterator i = engine.begin(); i != engine.end(); ++i) {
 
      const double x = i->getX();
     
      W += rayleigh(a,x) * i->getY() * 2*PI;
    }
 
    DEBUG("integral analytical " << FIXED(9,5) << V << endl);
    DEBUG("integral numerical  " << FIXED(9,5) << W << endl);
 
    ASSERT(fabs(V - W) <= precision, "Test Rayleigh integration.");
  }
 
  return 0;
}