1 #ifndef __JTOOLS__JQUADRATURE__
2 #define __JTOOLS__JQUADRATURE__
61 template<
class JFunction_t>
66 const double eps = 1.0e-4)
71 const double Vmin = integral(Xmin, Xmax) / (double) nx;
73 for (
int i = 0; i != nx; ++i) {
75 for (
double xmin = Xmin,
xmax = Xmax; ; ) {
78 const double v = integral(Xmin,
x);
80 if (fabs(Vmin -
v) < eps * Vmin ||
81 fabs(
xmax -
xmin) < eps * (Xmax - Xmin)) {
83 const double __x = 0.5 * (Xmin +
x);
84 const double __y = Vmin / integral(__x);
122 const double eps = 1.0e-12) :
127 const int M = (
n + 1) / 2;
129 double p0,
p1, p2, pp;
131 for (
int i = 0; i < M; ++i) {
133 double z = cos(
PI * (i+0.75) / (
n+0.5));
143 for (
int j = 0;
j <
n; ++
j) {
146 p2 = ((2*
j + 1) * z*
p1 -
j*p0) / (
j+1);
149 pp =
n * (z*p2 -
p1) / (z*z - 1.0);
154 }
while (fabs(z-z1) > eps);
156 const double y = 2.0 / ((1.0-z*z)*pp*pp);
185 const double eps = 1.0e-12) :
188 const int number_of_iterations = 100;
191 double p0,
p1, p2, pp;
193 double z = (1.0 + alf) * (3.0 + 0.92*alf) / (1.0 + 2.4*
n + 1.8*alf);
195 for (
int i = 0; i <
n; ++i) {
203 z += (15.0 + 6.25*alf) / (1.0 + 0.9*alf + 2.5*
n);
207 const double ai = i - 1;
208 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);
212 for (
int k = 0; k != number_of_iterations; ++k) {
219 for (
int j = 0;
j <
n; ++
j) {
222 p2 = ((2*
j + 1 + alf - z) *
p1 - (
j + alf)*p0) / (
j+1);
225 pp = (
n*p2 - (
n+alf)*
p1) / z;
230 if (fabs(z-z1) < eps)
234 const double y = -tgamma(alf+
n) / tgamma((
double)
n) / (pp*
n*
p1);
260 const double eps = 1.0e-12) :
265 const double pii = 1.0 /
pow(
PI,0.25);
267 const int number_of_iterations = 100;
269 const int M = (
n + 1) / 2;
271 double p0,
p1, p2, pp;
275 for (
int i = 0; i < M; ++i) {
280 z = sqrt((
double) (2*
n+1)) - 1.85575 *
pow((
double) (2*
n+1),-0.16667);
284 z -= 1.14 *
pow((
double)
n,0.426) / z;
288 z = 1.86*z + 0.86*at( 0 ).getX();
292 z = 1.91*z + 0.91*at( 1 ).getX();
296 z = 2.00*z + 1.00*at(i-2).getX();
300 for (
int k = 0; k != number_of_iterations; ++k) {
307 for (
int j = 0;
j <
n; ++
j) {
310 p2 = z * sqrt(2.0/(
double) (
j+1)) *
p1 - sqrt((
double)
j / (
double) (
j+1)) * p0;
313 pp = sqrt((
double) (2*
n)) *
p1;
318 if (fabs(z-z1) < eps)
322 const double y = 2.0 / (pp*pp);
352 const double b = -2*g * (
a + 1.0);
353 const double ai = 1.0 / (
a + 1.0);
355 const double ymin =
pow(1.0 + g, 2*(
a + 1.0)) / b;
356 const double ymax =
pow(1.0 - g, 2*(
a + 1.0)) / b;
358 const double dy = (ymax - ymin) / (
n + 1);
360 for (
double y = ymax - 0.5*dy;
y > ymin;
y -= dy) {
362 const double v =
y*b;
363 const double w =
pow(
v, ai);
364 const double x = (1.0 + g*g -
w) / (2*g);
365 const double dx =
pow(
v, -
a*ai)*dy;
388 const double b = -2*g * (
a + 1.0);
389 const double ai = 1.0 / (
a + 1.0);
391 const double ymin =
pow(1.0 + g*g -2*g*
xmin,
a + 1.0) / b;
392 const double ymax =
pow(1.0 + g*g -2*g*
xmax,
a + 1.0) / b;
394 const double dy = (ymax - ymin) / (
n + 1);
396 for (
double y = ymax - 0.5*dy;
y > ymin;
y -= dy) {
398 const double v =
y*b;
399 const double w =
pow(
v, ai);
400 const double x = (1.0 + g*g -
w) / (2*g);
401 const double dx =
pow(
v, -
a*ai)*dy;
418 const double dy = 1.0 / (
n + 1);
419 const double gi = log((1.0 + g*g) / (1.0 - g*g)) / (2.0*g);
421 for (
double y = 1.0 - 0.5*dy;
y > 0.0;
y -= dy) {
423 const double v = -
y*2.0*g*gi + log(1.0 + g*g);
424 const double w = exp(
v);
425 const double x = (1.0 + g*g -
w) / (2.0*g);
426 const double dx =
w*gi*dy;
453 const double dy = 1.0 / (
n + 1);
454 const double gi = 3.0/g + 1.0;
458 const double p = 1.0/g;
460 for (
double y = 0.5*dy;
y < 1.0;
y += dy) {
462 const double q = 0.5*gi - gi*
y;
464 const double b = sqrt(q*q + p*p*p);
465 const double u =
pow(-q + b, 1.0/3.0);
466 const double v =
pow(+q + b, 1.0/3.0);
468 const double x =
u -
v;
469 const double dx = (
u +
v) / (3.0*b);
494 for (
double ds = 1.0 / (
n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
496 const double cb = sqrt((1.0 + sb)*(1.0 - sb));
497 const double dc = ds*sb/cb;
527 for (ds = 1.0 / (
n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
529 cb = sqrt((1.0 + sb)*(1.0 - sb));
535 for (dc = (cb + 1.0) / (
n/2), cb -= 0.5*dc ; cb > -1.0; cb -= dc) {
General purpose class for a collection of sorted elements.
The elements in a collection are sorted according to their abscissa values and a given distance opera...
T pow(const T &x, const double y)
Power .
static const double PI
Mathematical constants.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).