1 #ifndef __JTOOLS__JQUADRATURE__
2 #define __JTOOLS__JQUADRATURE__
18 namespace JPP {
using namespace JTOOLS; }
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; ; ) {
77 const double x = 0.5 * (xmin + 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 for (
int i = 0; i <
n; ++i) {
198 z = (1.0 + alf) * (3.0 + 0.92*alf) / (1.0 + 2.4*n + 1.8*alf);
202 z += (15.0 + 6.25*alf) / (1.0 + 0.9*alf + 2.5*n);
206 const double ai = i - 1;
207 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);
214 for (k = 0; k != number_of_iterations; ++
k) {
221 for (
int j = 0;
j <
n; ++
j) {
224 p2 = ((2*
j + 1 + alf - z) * p1 - (
j + alf)*p0) / (
j+1);
227 pp = (n*p2 - (n+alf)*p1) / z;
232 if (fabs(z-z1) < eps)
236 const double y = -tgamma(alf+n) / tgamma((
double) n) / (pp*n*
p1);
262 const double eps = 1.0e-12) :
267 const double pii = 1.0 /
pow(
PI,0.25);
269 const int number_of_iterations = 100;
271 const int M = (n + 1) / 2;
273 double p0,
p1,
p2, pp;
277 for (
int i = 0; i <
M; ++i) {
282 z = sqrt((
double) (2*n+1)) - 1.85575 *
pow((
double) (2*n+1),-0.16667);
286 z -= 1.14 *
pow((
double) n,0.426) / z;
290 z = 1.86*z + 0.86*at( 0 ).getX();
294 z = 1.91*z + 0.91*at( 1 ).getX();
298 z = 2.00*z + 1.00*at(i-2).getX();
302 for (
int k = 0;
k != number_of_iterations; ++
k) {
309 for (
int j = 0;
j <
n; ++
j) {
312 p2 = z * sqrt(2.0/(
double) (
j+1)) * p1 - sqrt((
double)
j / (
double) (
j+1)) * p0;
315 pp = sqrt((
double) (2*n)) *
p1;
320 if (fabs(z-z1) < eps)
324 const double y = 2.0 / (pp*pp);
355 const double b = -2*g * (a + 1.0);
356 const double ai = 1.0 / (a + 1.0);
358 const double ymin =
pow(1.0 + g, 2*(a + 1.0)) / b;
359 const double ymax =
pow(1.0 - g, 2*(a + 1.0)) / b;
361 const double dy = (ymax - ymin) / (n + 1);
363 for (
double y = ymax - 0.5*dy; y > ymin; y -= dy) {
365 const double v = y*b;
366 const double w =
pow(v, ai);
367 const double x = (1.0 + g*g -
w) / (2*g);
368 const double dx =
pow(v, -a*ai)*dy;
391 const double b = -2*g * (a + 1.0);
392 const double ai = 1.0 / (a + 1.0);
394 const double ymin =
pow(1.0 + g*g -2*g*xmin, a + 1.0) / b;
395 const double ymax =
pow(1.0 + g*g -2*g*xmax, a + 1.0) / b;
397 const double dy = (ymax - ymin) / (n + 1);
399 for (
double y = ymax - 0.5*dy; y > ymin; y -= dy) {
401 const double v = y*b;
402 const double w =
pow(v, ai);
403 const double x = (1.0 + g*g -
w) / (2*g);
404 const double dx =
pow(v, -a*ai)*dy;
421 const double dy = 1.0 / (n + 1);
422 const double gi = log((1.0 + g*g) / (1.0 - g*g)) / (2.0*g);
424 for (
double y = 1.0 - 0.5*dy; y > 0.0; y -= dy) {
426 const double v = -y*2.0*g*gi + log(1.0 + g*g);
427 const double w =
exp(v);
428 const double x = (1.0 + g*g -
w) / (2.0*g);
429 const double dx = w*gi*dy;
457 const double dy = 1.0 / (n + 1);
458 const double gi = 3.0/g + 1.0;
462 const double p = 1.0/g;
464 for (
double y = 0.5*dy; y < 1.0; y += dy) {
466 const double q = 0.5*gi - gi*y;
468 const double b = sqrt(q*q + p*p*p);
469 const double u =
pow(-q + b, 1.0/3.0);
470 const double v =
pow(+q + b, 1.0/3.0);
472 const double x = u -
v;
473 const double dx = (u +
v) / (3.0*b);
498 for (
double ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
500 const double cb = sqrt((1.0 + sb)*(1.0 - sb));
501 const double dc = ds*sb/cb;
531 for (ds = 1.0 / (n/2), sb = 0.5*ds; sb < 1.0; sb += ds) {
533 cb = sqrt((1.0 + sb)*(1.0 - sb));
539 for (dc = (cb + 1.0) / (n/2), cb -= 0.5*dc ; cb > -1.0; cb -= dc) {
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
The elements in a collection are sorted according to their abscissa values and a given distance opera...
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
T pow(const T &x, const double y)
Power .
static const double PI
Mathematical constants.
General purpose class for a collection of sorted elements.
JTOOLS::JElement2D< double, double > JElement2D_t