68 double BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
69 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1,
70 BKMXX2 = 0.2 , BKMXY2 = 1 , BKMXX3 = 0.3 , BKMXY3 = 1;
72 double FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
73 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
74 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3);
76 double FNINV[] = {0, 1, 0.5, 0.33333333, 0.25, 0.2};
78 double EDGEC[]= {0, 0, 0.16666667e+0, 0.41666667e-1, 0.83333333e-2,
79 0.13888889e-1, 0.69444444e-2, 0.77160493e-3};
81 double U1[] = {0, 0.25850868e+0, 0.32477982e-1, -0.59020496e-2,
82 0. , 0.24880692e-1, 0.47404356e-2,
83 -0.74445130e-3, 0.73225731e-2, 0. ,
84 0.11668284e-2, 0. , -0.15727318e-2,-0.11210142e-2};
86 double U2[] = {0, 0.43142611e+0, 0.40797543e-1, -0.91490215e-2,
87 0. , 0.42127077e-1, 0.73167928e-2,
88 -0.14026047e-2, 0.16195241e-1, 0.24714789e-2,
89 0.20751278e-2, 0. , -0.25141668e-2,-0.14064022e-2};
91 double U3[] = {0, 0.25225955e+0, 0.64820468e-1, -0.23615759e-1,
92 0. , 0.23834176e-1, 0.21624675e-2,
93 -0.26865597e-2, -0.54891384e-2, 0.39800522e-2,
94 0.48447456e-2, -0.89439554e-2, -0.62756944e-2,-0.24655436e-2};
96 double U4[] = {0, 0.12593231e+1, -0.20374501e+0, 0.95055662e-1,
97 -0.20771531e-1, -0.46865180e-1, -0.77222986e-2,
98 0.32241039e-2, 0.89882920e-2, -0.67167236e-2,
99 -0.13049241e-1, 0.18786468e-1, 0.14484097e-1};
101 double U5[] = {0, -0.24864376e-1, -0.10368495e-2, 0.14330117e-2,
102 0.20052730e-3, 0.18751903e-2, 0.12668869e-2,
103 0.48736023e-3, 0.34850854e-2, 0. ,
104 -0.36597173e-3, 0.19372124e-2, 0.70761825e-3, 0.46898375e-3};
106 double U6[] = {0, 0.35855696e-1, -0.27542114e-1, 0.12631023e-1,
107 -0.30188807e-2, -0.84479939e-3, 0. ,
108 0.45675843e-3, -0.69836141e-2, 0.39876546e-2,
109 -0.36055679e-2, 0. , 0.15298434e-2, 0.19247256e-2};
111 double U7[] = {0, 0.10234691e+2, -0.35619655e+1, 0.69387764e+0,
112 -0.14047599e+0, -0.19952390e+1, -0.45679694e+0,
114 double U8[] = {0, 0.21487518e+2, -0.11825253e+2, 0.43133087e+1,
115 -0.14500543e+1, -0.34343169e+1, -0.11063164e+1,
116 -0.21000819e+0, 0.17891643e+1, -0.89601916e+0,
117 0.39120793e+0, 0.73410606e+0, 0. ,-0.32454506e+0};
119 double V1[] = {0, 0.27827257e+0, -0.14227603e-2, 0.24848327e-2,
120 0. , 0.45091424e-1, 0.80559636e-2,
121 -0.38974523e-2, 0. , -0.30634124e-2,
122 0.75633702e-3, 0.54730726e-2, 0.19792507e-2};
124 double V2[] = {0, 0.41421789e+0, -0.30061649e-1, 0.52249697e-2,
125 0. , 0.12693873e+0, 0.22999801e-1,
126 -0.86792801e-2, 0.31875584e-1, -0.61757928e-2,
127 0. , 0.19716857e-1, 0.32596742e-2};
129 double V3[] = {0, 0.20191056e+0, -0.46831422e-1, 0.96777473e-2,
130 -0.17995317e-2, 0.53921588e-1, 0.35068740e-2,
131 -0.12621494e-1, -0.54996531e-2, -0.90029985e-2,
132 0.34958743e-2, 0.18513506e-1, 0.68332334e-2,-0.12940502e-2};
134 double V4[] = {0, 0.13206081e+1, 0.10036618e+0, -0.22015201e-1,
135 0.61667091e-2, -0.14986093e+0, -0.12720568e-1,
136 0.24972042e-1, -0.97751962e-2, 0.26087455e-1,
137 -0.11399062e-1, -0.48282515e-1, -0.98552378e-2};
139 double V5[] = {0, 0.16435243e-1, 0.36051400e-1, 0.23036520e-2,
140 -0.61666343e-3, -0.10775802e-1, 0.51476061e-2,
141 0.56856517e-2, -0.13438433e-1, 0. ,
142 0. , -0.25421507e-2, 0.20169108e-2,-0.15144931e-2};
144 double V6[] = {0, 0.33432405e-1, 0.60583916e-2, -0.23381379e-2,
145 0.83846081e-3, -0.13346861e-1, -0.17402116e-2,
146 0.21052496e-2, 0.15528195e-2, 0.21900670e-2,
147 -0.13202847e-2, -0.45124157e-2, -0.15629454e-2, 0.22499176e-3};
149 double V7[] = {0, 0.54529572e+1, -0.90906096e+0, 0.86122438e-1,
150 0. , -0.12218009e+1, -0.32324120e+0,
151 -0.27373591e-1, 0.12173464e+0, 0. ,
154 double V8[] = {0, 0.93841352e+1, -0.16276904e+1, 0.16571423e+0,
155 0. , -0.18160479e+1, -0.50919193e+0,
156 -0.51384654e-1, 0.21413992e+0, 0. ,
159 double W1[] = {0, 0.29712951e+0, 0.97572934e-2, 0. ,
160 -0.15291686e-2, 0.35707399e-1, 0.96221631e-2,
161 -0.18402821e-2, -0.49821585e-2, 0.18831112e-2,
162 0.43541673e-2, 0.20301312e-2, -0.18723311e-2,-0.73403108e-3};
164 double W2[] = {0, 0.40882635e+0, 0.14474912e-1, 0.25023704e-2,
165 -0.37707379e-2, 0.18719727e+0, 0.56954987e-1,
166 0. , 0.23020158e-1, 0.50574313e-2,
167 0.94550140e-2, 0.19300232e-1};
169 double W3[] = {0, 0.16861629e+0, 0. , 0.36317285e-2,
170 -0.43657818e-2, 0.30144338e-1, 0.13891826e-1,
171 -0.58030495e-2, -0.38717547e-2, 0.85359607e-2,
172 0.14507659e-1, 0.82387775e-2, -0.10116105e-1,-0.55135670e-2};
174 double W4[] = {0, 0.13493891e+1, -0.26863185e-2, -0.35216040e-2,
175 0.24434909e-1, -0.83447911e-1, -0.48061360e-1,
176 0.76473951e-2, 0.24494430e-1, -0.16209200e-1,
177 -0.37768479e-1, -0.47890063e-1, 0.17778596e-1, 0.13179324e-1};
179 double W5[] = {0, 0.10264945e+0, 0.32738857e-1, 0. ,
180 0.43608779e-2, -0.43097757e-1, -0.22647176e-2,
181 0.94531290e-2, -0.12442571e-1, -0.32283517e-2,
182 -0.75640352e-2, -0.88293329e-2, 0.52537299e-2, 0.13340546e-2};
184 double W6[] = {0, 0.29568177e-1, -0.16300060e-2, -0.21119745e-3,
185 0.23599053e-2, -0.48515387e-2, -0.40797531e-2,
186 0.40403265e-3, 0.18200105e-2, -0.14346306e-2,
187 -0.39165276e-2, -0.37432073e-2, 0.19950380e-2, 0.12222675e-2};
189 double W8[] = {0, 0.66184645e+1, -0.73866379e+0, 0.44693973e-1,
190 0. , -0.14540925e+1, -0.39529833e+0,
191 -0.44293243e-1, 0.88741049e-1};
194 if (fKappa <0.01 || fKappa >12) {
195 std::cerr <<
"VavilovFast::set: illegal value of kappa=" << kappa << std::endl;
204 double x,
y, xx, yy,
x2,
x3, y2, y3,
xy, p2, p3, q2, q3, pq;
214 for (j=1; j<=4; j++) {
215 DRK[j+1] = DRK[1]*DRK[j];
216 DSIGM[j+1] = DSIGM[1]*DSIGM[j];
217 ALFA[j+1] = (FNINV[j]-
fBeta2*FNINV[j+1])*DRK[j];
221 fHC[2]=ALFA[3]*DSIGM[3];
222 fHC[3]=(3*ALFA[2]*ALFA[2] + ALFA[4])*DSIGM[4]-3;
223 fHC[4]=(10*ALFA[2]*ALFA[3]+ALFA[5])*DSIGM[5]-10*
fHC[2];
248 fAC[1] = W1[1] + W1[2]*
x + W1[4]*
x3 + W1[5]*
y + W1[6]*y2 + W1[7]*y3 +
249 W1[8]*
xy + W1[9]*p2 + W1[10]*p3 + W1[11]*q2 + W1[12]*q3 + W1[13]*pq;
250 fAC[2] = W2[1] + W2[2]*
x + W2[3]*
x2 + W2[4]*
x3 + W2[5]*
y + W2[6]*y2 +
251 W2[8]*
xy + W2[9]*p2 + W2[10]*p3 + W2[11]*q2;
252 fAC[3] = W3[1] + W3[3]*
x2 + W3[4]*
x3 + W3[5]*
y + W3[6]*y2 + W3[7]*y3 +
253 W3[8]*
xy + W3[9]*p2 + W3[10]*p3 + W3[11]*q2 + W3[12]*q3 + W3[13]*pq;
254 fAC[4] = W4[1] + W4[2]*
x + W4[3]*
x2 + W4[4]*
x3 + W4[5]*
y + W4[6]*y2 + W4[7]*y3 +
255 W4[8]*
xy + W4[9]*p2 + W4[10]*p3 + W4[11]*q2 + W4[12]*q3 + W4[13]*pq;
256 fAC[5] = W5[1] + W5[2]*
x + W5[4]*
x3 + W5[5]*
y + W5[6]*y2 + W5[7]*y3 +
257 W5[8]*
xy + W5[9]*p2 + W5[10]*p3 + W5[11]*q2 + W5[12]*q3 + W5[13]*pq;
258 fAC[6] = W6[1] + W6[2]*
x + W6[3]*
x2 + W6[4]*
x3 + W6[5]*
y + W6[6]*y2 + W6[7]*y3 +
259 W6[8]*
xy + W6[9]*p2 + W6[10]*p3 + W6[11]*q2 + W6[12]*q3 + W6[13]*pq;
260 fAC[8] = W8[1] + W8[2]*
x + W8[3]*
x2 + W8[5]*
y + W8[6]*y2 + W8[7]*y3 + W8[8]*
xy;
262 }
else if (
fKappa >= 0.12) {
279 fAC[1] = V1[1] + V1[2]*
x + V1[3]*
x2 + V1[5]*
y + V1[6]*y2 + V1[7]*y3 +
280 V1[9]*p2 + V1[10]*p3 + V1[11]*q2 + V1[12]*q3;
281 fAC[2] = V2[1] + V2[2]*
x + V2[3]*
x2 + V2[5]*
y + V2[6]*y2 + V2[7]*y3 +
282 V2[8]*
xy + V2[9]*p2 + V2[11]*q2 + V2[12]*q3;
283 fAC[3] = V3[1] + V3[2]*
x + V3[3]*
x2 + V3[4]*
x3 + V3[5]*
y + V3[6]*y2 + V3[7]*y3 +
284 V3[8]*
xy + V3[9]*p2 + V3[10]*p3 + V3[11]*q2 + V3[12]*q3 + V3[13]*pq;
285 fAC[4] = V4[1] + V4[2]*
x + V4[3]*
x2 + V4[4]*
x3 + V4[5]*
y + V4[6]*y2 + V4[7]*y3 +
286 V4[8]*
xy + V4[9]*p2 + V4[10]*p3 + V4[11]*q2 + V4[12]*q3;
287 fAC[5] = V5[1] + V5[2]*
x + V5[3]*
x2 + V5[4]*
x3 + V5[5]*
y + V5[6]*y2 + V5[7]*y3 +
288 V5[8]*
xy + V5[11]*q2 + V5[12]*q3 + V5[13]*pq;
289 fAC[6] = V6[1] + V6[2]*
x + V6[3]*
x2 + V6[4]*
x3 + V6[5]*
y + V6[6]*y2 + V6[7]*y3 +
290 V6[8]*
xy + V6[9]*p2 + V6[10]*p3 + V6[11]*q2 + V6[12]*q3 + V6[13]*pq;
291 fAC[7] = V7[1] + V7[2]*
x + V7[3]*
x2 + V7[5]*
y + V7[6]*y2 + V7[7]*y3 +
292 V7[8]*
xy + V7[11]*q2;
293 fAC[8] = V8[1] + V8[2]*
x + V8[3]*
x2 + V8[5]*
y + V8[6]*y2 + V8[7]*y3 +
294 V8[8]*
xy + V8[11]*q2;
315 fAC[1] = U1[1] + U1[2]*
x + U1[3]*
x2 + U1[5]*
y + U1[6]*y2 + U1[7]*y3 +
316 U1[8]*
xy + U1[10]*p3 + U1[12]*q3 + U1[13]*pq;
317 fAC[2] = U2[1] + U2[2]*
x + U2[3]*
x2 + U2[5]*
y + U2[6]*y2 + U2[7]*y3 +
318 U2[8]*
xy + U2[9]*p2 + U2[10]*p3 + U2[12]*q3 + U2[13]*pq;
319 fAC[3] = U3[1] + U3[2]*
x + U3[3]*
x2 + U3[5]*
y + U3[6]*y2 + U3[7]*y3 +
320 U3[8]*
xy + U3[9]*p2 + U3[10]*p3 + U3[11]*q2 + U3[12]*q3 + U3[13]*pq;
321 fAC[4] = U4[1] + U4[2]*
x + U4[3]*
x2 + U4[4]*
x3 + U4[5]*
y + U4[6]*y2 + U4[7]*y3 +
322 U4[8]*
xy + U4[9]*p2 + U4[10]*p3 + U4[11]*q2 + U4[12]*q3;
323 fAC[5] = U5[1] + U5[2]*
x + U5[3]*
x2 + U5[4]*
x3 + U5[5]*
y + U5[6]*y2 + U5[7]*y3 +
324 U5[8]*
xy + U5[10]*p3 + U5[11]*q2 + U5[12]*q3 + U5[13]*pq;
325 fAC[6] = U6[1] + U6[2]*
x + U6[3]*
x2 + U6[4]*
x3 + U6[5]*
y + U6[7]*y3 +
326 U6[8]*
xy + U6[9]*p2 + U6[10]*p3 + U6[12]*q3 + U6[13]*pq;
327 fAC[7] = U7[1] + U7[2]*
x + U7[3]*
x2 + U7[4]*
x3 + U7[5]*
y + U7[6]*y2 + U7[8]*
xy;
329 fAC[8] = U8[1] + U8[2]*
x + U8[3]*
x2 + U8[4]*
x3 + U8[5]*
y + U8[6]*y2 + U8[7]*y3 +
330 U8[8]*
xy + U8[9]*p2 + U8[10]*p3 + U8[11]*q2 + U8[13]*pq;
352 for (k=1; k<=
fNpt; k++) {
359 for (k=1; k<=
fNpt; k++)
375 double xx = (
x +
fHC[0])*
fHC[1];
378 for (k=2; k<=8; k++) {
380 h[k+1] = xx*
h[k]-fn*
h[k-1];
382 double s = 1 +
fHC[7]*
h[9];
428 if (
x <
fAC[0])
v = 0;
429 else if (
x >=
fAC[8])
v = 1;
480 if (z < 0 || z > 1)
return std::numeric_limits<double>::signaling_NaN();
484 double t = 2*z/
fAC[9];
485 double rlam =
fAC[0];
490 for (
int n = 1;
n <=
fNpt; ++
n) {
494 double x = (rlam+
fHC[0])*
fHC[1];
497 for (
int k = 2; k <= 8; ++k) {
499 h[k+1] =
x*
h[k]-fn*
h[k-1];
501 double y = 1+
fHC[7]*
h[9];
502 for (
int k = 2; k <= 6; ++k) {
508 double x = rlam*rlam;
514 double x = rlam*rlam;
530 double v = rlam-
fAC[9];
541 if (z < 0 || z > 1)
return std::numeric_limits<double>::signaling_NaN();
579 return vavilov->
Pdf (
x);
584 return vavilov->
Cdf (
x);
589 return vavilov->
Cdf_c (
x);