marley::Integrator Class Reference

Numerical integrator that uses Clenshaw-Curtis quadrature More...

#include <Integrator.hh>

Public Member Functions
	Integrator (size_t num=N_DEFAULT_)
	Create a Clenshaw-Curtis quadrature integrator that uses 2*num sampling points.
double	num_integrate (const std::function< double(double)> &f, double a, double b) const
	Numerically integrate an arbitrary 1D function.

Detailed Description

Numerical integrator that uses Clenshaw-Curtis quadrature

Definition at line 28 of file Integrator.hh.

Constructor & Destructor Documentation

Integrator()

marley::Integrator::Integrator

(

size_t

num = N_DEFAULT_

)

Create a Clenshaw-Curtis quadrature integrator that uses 2*num sampling points.

Parameters

num	half the number of sampling points to use

Todo

add error check for when num == 0 or too small

Todo

Todo

add error check for when num is ridiculously large

Definition at line 22 of file Integrator.cc.

                                       : N_(num), weights_(num + 1, 0.),
  offsets_(num - 1)
{
 
  // Precompute the N_ - 1 offsets for speed
  double arg = 0.;
  for (size_t n = 0; n < N_ - 1; ++n) {
    arg += marley_utils::half_pi / N_;
    offsets_[n] = std::cos(arg);
  }
 
  // Also precompute the N_ + 1 weights
  for (size_t n = 0; n <= N_; ++n) {
    for (size_t k = 0; k <= N_; ++k) {
      double weight_piece = std::cos(n * k * marley_utils::pi / N_)
        / (1. - std::pow(2*k, 2));
      if (k != 0 && k != N_) weight_piece *= 2.;
      weights_[n] += weight_piece;
    }
    weights_[n] /= N_;
  }
}

Member Function Documentation

num_integrate()

double marley::Integrator::num_integrate	(	const std::function< double(double)> &	f,
		double	a,
		double	b ) const

Numerically integrate an arbitrary 1D function.

Numerically integrate a std::function<double(double)> over the interval [a,b] using Clenshaw-Curtis quadrature at 2N_ sampling points.

Definition at line 47 of file Integrator.cc.

{
  double A = (b - a) / 2.;
  double B = (b + a) / 2.;
  double C = (f(a) + f(b)) / 2.;
 
  double integral = weights_[0] * C; // n = 0 term
  integral += weights_[N_] * f(B); // n = N_ term
 
  // n = 1 to n = N_ - 1 terms
  for (size_t n = 1; n < N_; ++n) {
    double epoint = A * offsets_[n - 1];
    integral += weights_[n] * (f(B + epoint) + f(B - epoint));
  }
 
  return A * integral;
}

Referenced by marley::DiscreteNuclearReaction::total_xs().

---

The documentation for this class was generated from the following files:

• /home/sjg/Desktop/marley/include/marley/Integrator.hh
• /home/sjg/Desktop/marley/src/Integrator.cc