13 #ifndef SI_PIXEL_TEMPLATE_STANDALONE
23 #include<boost/bind.hpp>
26 namespace VVIObjDetails {
27 void sincosint(
double x,
double & sint,
double & cint);
32 inline double f1(
double x,
double const * h_) {
return h_[0]+h_[1]*
std::log(h_[2]*x)-h_[3]*
x;}
35 int dzero(
double a,
double b,
double& x0,
36 double& rv,
double eps,
int mxf, F func);
51 const double xp[9] = { 9.29,2.47,.89,.36,.15,.07,.03,.02,0.0 };
52 const double xq[7] = { .012,.03,.08,.26,.87,3.83,11.0 };
54 double q, u,
x,
c1, c2, c3, c4,
d1, h4, h5, h6,
q2, x1, d, ll, ul, xf1, xf2, rv;
59 if(kappa < 0.01) kappa = 0.01;
60 if(kappa > 10.) kappa = 10.;
61 if(beta2 < 0.) beta2 = 0.;
62 if(beta2 > 1.) beta2 = 1.;
64 h_[4] = 1. - beta2*0.42278433999999998 + 7.6/kappa;
67 h4 = -7.6/kappa - (beta2 * .57721566 + 1);
74 for (lp = 0; lp < 9; ++lp) {
75 if (kappa >= xp[lp])
break;
78 for (lq = 0; lq < 7; ++lq) {
79 if (kappa <= xq[lq])
break;
88 h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
89 if (kappa >= .07) {h_[0] += 6.90775527;}
90 h_[1] = beta2 * kappa;
92 h_[3] = omega_ * 1.5707963250000001;
96 d =
exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
99 a_[n - 1] = omega_ * .31830988654751274;
103 for (k = 1; k <
n; ++
k) {
111 xf1 = kappa * (beta2 * c1 - c4) - x * c2;
112 xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
114 d1 = q * d * omega_ *
exp(xf1);
115 a_[l - 1] = d1 *
cos(xf2);
116 b_[l - 1] = -d1 *
sin(xf2);
118 d1 = q * d *
exp(xf1)/
k;
119 a_[l - 1] = d1 *
sin(xf2);
120 b_[l - 1] = d1 *
cos(xf2);
121 a_[n - 1] += q2 *
a_[l - 1];
140 double f, u,
y, a0, a1;
142 double b1, b0, b2, cof;
148 }
else if (x <=
t1_) {
150 u =
omega_ * y - 3.141592653589793;
155 for (k = 2; k <= n1; ++
k) {
158 a0 =
a_[k - 1] + cof * a1 - a2;
162 for (k = 2; k <=
n; ++
k) {
165 b0 =
b_[k - 1] + cof * b1 - b2;
167 f = (a0 - a2) * .5 + b0 *
sin(u);
171 if (
mode_ != 0) {f = 1.;}
193 namespace VVIObjDetails {
197 const double zero = 0.;
198 const double one = 1.;
199 const double two = 2.;
200 const double eight = 8.;
201 const double ce = .57721566490153;
202 const double c__[14] = { 1.9405491464836,.9413409132865,
203 -.579845034293,.3091572011159,-.0916101792208,.0164437407515,
204 -.0019713091952,1.692538851e-4,-1.09393296e-5,5.522386e-7,
205 -2.23995e-8,7.465e-10,-2.08e-11,5
e-13 };
206 const double p[23] = { .96074783975204,-.0371138962124,
207 .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
208 6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
209 -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7
e-12,2.54e-12,
210 -9.5e-13,3.7e-13,-1.4e-13,6
e-14,-2
e-14,1e-14 };
211 const double q[20] = { .98604065696238,-.0134717382083,
212 4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
213 6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
214 -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6
e-14,
223 double r__,
y, b0, b1, b2,
pp, qq, alfa;
230 if (fabs(x) <= eight) {
234 h__ = two * (d__1 * d__1) - one;
238 for (i__ = 13; i__ >= 0; --i__) {
239 b0 = c__[i__] - alfa * b1 - b2;
243 b1 = ce +
log((fabs(x))) - b0 + h__ * b2;
249 h__ = two * (d__1 * d__1) - one;
253 for (i__ = 22; i__ >= 0; --i__) {
254 b0 = p[i__] - alfa * b1 - b2;
261 for (i__ = 19; i__ >= 0; --i__) {
262 b0 = q[i__] - alfa * b1 - b2;
267 b1 = r__ * (qq *
sin(x) - r__ * pp *
cos(x));
275 const double zero = 0.;
276 const double one = 1.;
277 const double two = 2.;
278 const double eight = 8.;
279 const double pih = 1.5707963267949;
280 const double s[14] = { 1.9522209759531,-.6884042321257,
281 .4551855132256,-.1804571236838,.0410422133759,-.0059586169556,
282 6.001427414e-4,-4.44708329e-5,2.5300782e-6,-1.141308e-7,4.1858e-9,
283 -1.273e-10,3.3e-12,-1
e-13 };
284 const double p[23] = { .96074783975204,-.0371138962124,
285 .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
286 6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
287 -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7
e-12,2.54e-12,
288 -9.5e-13,3.7e-13,-1.4e-13,6
e-14,-2
e-14,1e-14 };
289 const double q[20] = { .98604065696238,-.0134717382083,
290 4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
291 6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
292 -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6
e-14,
301 double r__,
y, b0, b1, b2,
pp, qq, alfa;
303 if (fabs(x) <= eight) {
306 h__ = two * (d__1 * d__1) - one;
310 for (i__ = 13; i__ >= 0; --i__) {
311 b0 = s[i__] - alfa * b1 - b2;
320 h__ = two * (d__1 * d__1) - one;
324 for (i__ = 22; i__ >= 0; --i__) {
325 b0 = p[i__] - alfa * b1 - b2;
332 for (i__ = 19; i__ >= 0; --i__) {
333 b0 = q[i__] - alfa * b1 - b2;
339 if(x < 0.) d__1 = -d__1;
340 b1 = d__1 - r__ * (r__ * pp *
sin(x) + qq *
cos(x));
349 const double zero = 0.;
350 const double one = 1.;
351 const double two = 2.;
352 const double eight = 8.;
353 const double ce = .57721566490153;
354 const double pih = 1.5707963267949;
355 const double s__[14] = { 1.9522209759531,-.6884042321257,
356 .4551855132256,-.1804571236838,.0410422133759,-.0059586169556,
357 6.001427414e-4,-4.44708329e-5,2.5300782e-6,-1.141308e-7,4.1858e-9,
358 -1.273e-10,3.3e-12,-1
e-13 };
360 const double c__[14] = { 1.9405491464836,.9413409132865,
361 -.579845034293,.3091572011159,-.0916101792208,.0164437407515,
362 -.0019713091952,1.692538851e-4,-1.09393296e-5,5.522386e-7,
363 -2.23995e-8,7.465e-10,-2.08e-11,5
e-13 };
365 const double p[23] = { .96074783975204,-.0371138962124,
366 .00194143988899,-1.7165988425e-4,2.112637753e-5,-3.27163257e-6,
367 6.0069212e-7,-1.2586794e-7,2.932563e-8,-7.45696e-9,2.04105e-9,
368 -5.9502e-10,1.8323e-10,-5.921e-11,1.997e-11,-7
e-12,2.54e-12,
369 -9.5e-13,3.7e-13,-1.4e-13,6
e-14,-2
e-14,1e-14 };
370 const double q[20] = { .98604065696238,-.0134717382083,
371 4.5329284117e-4,-3.067288652e-5,3.13199198e-6,-4.2110196e-7,
372 6.907245e-8,-1.318321e-8,2.83697e-9,-6.7329e-10,1.734e-10,
373 -4.787e-11,1.403e-11,-4.33e-12,1.4e-12,-4.7e-13,1.7e-13,-6
e-14,
382 double r__,
y, b0, b1, b2,
pp, qq, alfa;
388 if (fabs(x) <= eight) {
392 h__ = two * (d__1 * d__1) - one;
399 for (i__ = 13; i__ >= 0; --i__) {
400 b0 = c__[i__] - alfa * b1 - b2;
404 cint = ce +
log((fabs(x))) - b0 + h__ * b2;
409 for (i__ = 13; i__ >= 0; --i__) {
410 b0 = s__[i__] - alfa * b1 - b2;
414 sint = y * (b0 - b2);
421 h__ = two * (d__1 * d__1) - one;
425 for (i__ = 22; i__ >= 0; --i__) {
426 b0 = p[i__] - alfa * b1 - b2;
433 for (i__ = 19; i__ >= 0; --i__) {
434 b0 = q[i__] - alfa * b1 - b2;
440 cint = r__ * (qq *
sin(x) - r__ * pp *
cos(x));
443 if(x < 0.) d__1 = -d__1;
444 sint = d__1 - r__ * (r__ * pp *
sin(x) + qq *
cos(x));
453 const double zero = 0.;
454 const double q2[7] = { .10340013040487,3.319092135933,
455 20.449478501379,41.280784189142,32.426421069514,10.041164382905,
457 const double p3[6] = { -2.3909964453136,-147.98219500504,
458 -254.3763397689,-119.55761038372,-19.630408535939,-.9999999999036
460 const double q3[6] = { 177.60070940351,530.68509610812,
461 462.23027156148,156.81843364539,21.630408494238,1. };
462 const double p4[8] = { -8.6693733995107,-549.14226552109,
463 -4210.0161535707,-249301.39345865,-119623.66934925,
464 -22174462.775885,3892804.213112,-391546073.8091 };
465 const double q4[8] = { 34.171875,-1607.0892658722,35730.029805851,
466 -483547.43616216,4285596.2461175,-24903337.574054,89192576.757561,
468 const double a1[8] = { -2.1808638152072,-21.901023385488,
469 9.3081638566217,25.076281129356,-33.184253199722,60.121799083008,
470 -43.253113287813,1.0044310922808 };
471 const double b1[8] = { 0.,3.9370770185272,300.89264837292,
472 -6.2504116167188,1003.6743951673,14.325673812194,2736.2411988933,
474 const double a2[8] = { -3.4833465360285,-18.65454548834,
475 -8.2856199414064,-32.34673303054,17.960168876925,1.7565631546961,
476 -1.9502232128966,.99999429607471 };
477 const double b2[8] = { 0.,69.500065588743,57.283719383732,
478 25.777638423844,760.76114800773,28.951672792514,-3.4394226689987,
480 const double a3[6] = { -27.780928934438,-10.10479081576,
481 -9.1483008216736,-5.0223317461851,-3.0000077799358,
483 const double one = 1.;
484 const double b3[6] = { 0.,122.39993926823,2.7276100778779,
485 -7.1897518395045,-2.9990118065262,1.999999942826 };
486 const double two = 2.;
487 const double three = 3.;
488 const double x0 = .37250741078137;
489 const double xl[6] = { -24.,-12.,-6.,0.,1.,4. };
490 const double p1[5] = { 4.293125234321,39.894153870321,
491 292.52518866921,425.69682638592,-434.98143832952 };
492 const double q1[5] = { 1.,18.899288395003,150.95038744251,
493 568.05252718987,753.58564359843 };
494 const double p2[7] = { .43096783946939,6.9052252278444,
495 23.019255939133,24.378408879132,9.0416155694633,.99997957705159,
499 double v,
y, ap, bp, aq, dp, bq, dq;
503 for (
int i__ = 2; i__ <= 5; ++i__) {
505 ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
507 y =
exp(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
508 }
else if (x <= xl[1]) {
510 for (
int i__ = 2; i__ <= 7; ++i__) {
511 ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
513 y =
exp(-x) / x * (a2[7] + b2[7] / ap);
514 }
else if (x <= xl[2]) {
516 for (
int i__ = 2; i__ <= 7; ++i__) {
517 ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
519 y =
exp(-x) / x * (a1[7] + b1[7] / ap);
520 }
else if (x < xl[3]) {
521 v = -two * (x / three + one);
524 for (
int i__ = 2; i__ <= 8; ++i__) {
527 dp = p4[i__ - 1] - ap + v * bp;
531 for (
int i__ = 2; i__ <= 8; ++i__) {
534 dq = q4[i__ - 1] - aq + v * bq;
536 y = -
log(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
537 }
else if (x == xl[3]) {
539 }
else if (x < xl[4]) {
542 for (
int i__ = 2; i__ <= 5; ++i__) {
543 ap = p1[i__ - 1] + x * ap;
544 aq = q1[i__ - 1] + x * aq;
546 y = -
log(x) + ap / aq;
547 }
else if (x <= xl[5]) {
551 for (
int i__ = 2; i__ <= 7; ++i__) {
552 ap = p2[i__ - 1] + y * ap;
553 aq = q2[i__ - 1] + y * aq;
555 y =
exp(-x) * ap / aq;
560 for (
int i__ = 2; i__ <= 6; ++i__) {
561 ap = p3[i__ - 1] + y * ap;
562 aq = q3[i__ - 1] + y * aq;
564 y =
exp(-x) * y * (one + y * ap / aq);
572 double& rv,
double eps,
int mxf, F func) {
574 double d__1, d__2, d__3, d__4;
577 double f1,
f2,
f3, u1, u2, x1, x2, u3, u4, x3, ca, cb, cc,
fa,
fb, ee,
ff;
579 double xa, xb, fx, xx, su4;
594 ee = eps * (fabs(x0) + 1);
608 rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
624 if (u2 == 0. || u4 == 0.) {
goto L1;}
630 cb = (x1 + x2) * u2 - (x2 + x0) * u1;
631 cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
633 if (cb == 0.) {
goto L1;}
637 u4 = u3 * u3 - cc / ca;
638 if (u4 < 0.) {
goto L1;}
640 if (x0 + u3 < 0.
f) {su4 = -su4;}
643 if (x0 < xa || x0 > xb) {
goto L1;}
645 d__3 = (d__1 = x0 - x3, fabs(d__1)), d__4 = (d__2 = x0 - x2, fabs(d__2));
647 ee = eps * (fabs(x0) + 1);
685 rv = (d__1 = xb - xa, fabs(d__1)) * -.5;
void limits(double &xl, double &xu) const
density (mode=0) or distribution (mode=1) function
Sin< T >::type sin(const T &t)
Exp< T >::type exp(const T &t)
double sinint(double x)
Private version of the cosine integral.
int dzero(double a, double b, double &x0, double &rv, double eps, int mxf, F func)
VVIObj(double kappa=0.01, double beta2=1., double mode=0)
Constructor.
double cosint(double x)
Private version of the cosine and sine integral.
const T & max(const T &a, const T &b)
Cos< T >::type cos(const T &t)
double f2(double x, double const *h_)
double fcn(double x) const
void sincosint(double x, double &sint, double &cint)
Log< T >::type log(const T &t)
double f1(double x, double const *h_)
Private version of the exponential integral.
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
double expint(double x)
Private version of the sine integral.