14 #ifndef SI_PIXEL_TEMPLATE_STANDALONE 24 #include "vdt/vdtMath.h" 27 void sincosint(
float x,
float & sint,
float & cint);
31 int dzero(
float a,
float b,
float& x0,
32 float& rv,
float eps,
int mxf,
F func);
47 const float xp[9] = { 9.29,2.47,.89,.36,.15,.07,.03,.02,0.0 };
48 const float xq[7] = { .012,.03,.08,.26,.87,3.83,11.0 };
50 float q, u,
x,
c1, c2, c3, c4, d1, h4, h5, h6,
q2, x1,
d, ll, ul, xf1, xf2, rv;
55 if(kappa < 0.01
f) kappa = 0.01f;
56 if(kappa > 10.
f) kappa = 10.f;
57 if(beta2 < 0.
f) beta2 = 0.f;
58 if(beta2 > 1.
f) beta2 = 1.f;
60 float invKappa = 1.f/
kappa;
61 h_[4] = 1.f - beta2*0.42278433999999998f + (7.6f*invKappa);
64 h4 = - (7.6f*invKappa) - (beta2 * .57721566
f + 1.
f);
71 for (lp = 0; lp < 9; ++lp) {
72 if (kappa >= xp[lp])
break;
74 ll = -
float(lp) - 1.5f;
75 for (lq = 0; lq < 7; ++lq) {
76 if (kappa <= xq[lq])
break;
85 h_[0] = kappa * (beta2 * .57721566f + 2.f) + 9.9166128600000008
f;
86 if (kappa >= .07) {h_[0] += 6.90775527f;}
87 h_[1] = beta2 *
kappa;
89 h_[3] = omega_ * 1.5707963250000001f;
93 d =
vdt::fast_expf(kappa * (beta2 * (.57721566
f - h5) + 1.
f)) * .31830988654751274f;
96 a_[n - 1] = omega_ * .31830988654751274f;
100 for (k = 1; k <
n; ++
k) {
106 vdt::fast_sincosf(x1,c3,c4);
107 xf1 = kappa * (beta2 * c1 - c4) - x * c2;
108 xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
109 float s,
c; vdt::fast_sincosf(xf2,s,c);
118 a_[n - 1] += q2 *
a_[l - 1];
137 float f, u,
y, a0, a1;
139 float b1, b0, b2, cof;
145 }
else if (x <=
t1_) {
147 u =
omega_ * y - 3.141592653589793f;
148 float su,cu; vdt::fast_sincosf(u,su,cu);
153 for (k = 2; k <= n1; ++
k) {
156 a0 =
a_[k - 1] + cof * a1 - a2;
160 for (k = 2; k <=
n; ++
k) {
163 b0 =
b_[k - 1] + cof * b1 - b2;
165 f = (a0 - a2) * .5f + b0 * su;
169 if (
mode_ != 0) {f = 1.f;}
202 const float zero = 0.;
203 const float q2[7] = { .10340013040487,3.319092135933,
204 20.449478501379,41.280784189142,32.426421069514,10.041164382905,
206 const float p3[6] = { -2.3909964453136,-147.98219500504,
207 -254.3763397689,-119.55761038372,-19.630408535939,-.9999999999036
209 const float q3[6] = { 177.60070940351,530.68509610812,
210 462.23027156148,156.81843364539,21.630408494238,1. };
211 const float p4[8] = { -8.6693733995107,-549.14226552109,
212 -4210.0161535707,-249301.39345865,-119623.66934925,
213 -22174462.775885,3892804.213112,-391546073.8091 };
214 const float q4[8] = { 34.171875,-1607.0892658722,35730.029805851,
215 -483547.43616216,4285596.2461175,-24903337.574054,89192576.757561,
217 const float a1[8] = { -2.1808638152072,-21.901023385488,
218 9.3081638566217,25.076281129356,-33.184253199722,60.121799083008,
219 -43.253113287813,1.0044310922808 };
220 const float b1[8] = { 0.,3.9370770185272,300.89264837292,
221 -6.2504116167188,1003.6743951673,14.325673812194,2736.2411988933,
223 const float a2[8] = { -3.4833465360285,-18.65454548834,
224 -8.2856199414064,-32.34673303054,17.960168876925,1.7565631546961,
225 -1.9502232128966,.99999429607471 };
226 const float b2[8] = { 0.,69.500065588743,57.283719383732,
227 25.777638423844,760.76114800773,28.951672792514,-3.4394226689987,
229 const float a3[6] = { -27.780928934438,-10.10479081576,
230 -9.1483008216736,-5.0223317461851,-3.0000077799358,
232 const float one = 1.;
233 const float b3[6] = { 0.,122.39993926823,2.7276100778779,
234 -7.1897518395045,-2.9990118065262,1.999999942826 };
235 const float two = 2.;
236 const float three = 3.;
237 const float x0 = .37250741078137;
238 const float xl[6] = { -24.,-12.,-6.,0.,1.,4. };
239 const float p1[5] = { 4.293125234321,39.894153870321,
240 292.52518866921,425.69682638592,-434.98143832952 };
241 const float q1[5] = { 1.,18.899288395003,150.95038744251,
242 568.05252718987,753.58564359843 };
243 const float p2[7] = { .43096783946939,6.9052252278444,
244 23.019255939133,24.378408879132,9.0416155694633,.99997957705159,
248 float v,
y, ap, bp, aq,
dp, bq, dq;
252 for (
int i__ = 2; i__ <= 5; ++i__) {
253 ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
256 }
else if (x <= xl[1]) {
258 for (
int i__ = 2; i__ <= 7; ++i__) {
259 ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
262 }
else if (x <= xl[2]) {
264 for (
int i__ = 2; i__ <= 7; ++i__) {
265 ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
268 }
else if (x < xl[3]) {
269 v = -two * (x / three + one);
272 for (
int i__ = 2; i__ <= 8; ++i__) {
275 dp = p4[i__ - 1] - ap + v * bp;
279 for (
int i__ = 2; i__ <= 8; ++i__) {
282 dq = q4[i__ - 1] - aq + v * bq;
285 }
else if (x == xl[3]) {
287 }
else if (x < xl[4]) {
290 for (
int i__ = 2; i__ <= 5; ++i__) {
291 ap = p1[i__ - 1] + x * ap;
292 aq = q1[i__ - 1] + x * aq;
295 }
else if (x <= xl[5]) {
299 for (
int i__ = 2; i__ <= 7; ++i__) {
300 ap = p2[i__ - 1] + y * ap;
301 aq = q2[i__ - 1] + y * aq;
308 for (
int i__ = 2; i__ <= 6; ++i__) {
309 ap = p3[i__ - 1] + y * ap;
310 aq = q3[i__ - 1] + y * aq;
322 float& rv,
float eps,
int mxf,
F func) {
324 float d__1, d__2, d__3, d__4;
327 float f1,
f2,
f3,
u1,
u2, x1, x2, u3, u4, x3, ca, cb, cc,
fa,
fb, ee,
ff;
329 float xa, xb, fx,
xx, su4;
336 rv = (xb - xa) * -2.
f;
342 x0 = (xa + xb) * 0.5
f;
358 rv = (d__1 = xb - xa, fabs(d__1)) * -0.5
f;
374 if (u2 == 0.
f || u4 == 0.
f) {
goto L1;}
380 cb = (x1 + x2) * u2 - (x2 + x0) *
u1;
381 cc = (x1 - x0) * f1 - x1 * (ca * x1 + cb);
383 if (cb == 0.
f) {
goto L1;}
386 u3 = cb / (ca * 2.f);
387 u4 = u3 * u3 - cc / ca;
388 if (u4 < 0.
f) {
goto L1;}
390 if (x0 + u3 < 0.
f) {su4 = -su4;}
393 if (x0 < xa || x0 > xb) {
goto L1;}
395 d__3 = (d__1 = x0 - x3,
std::abs(d__1));
396 d__4 = (d__2 = x0 - x2,
std::abs(d__2));
436 rv = (d__1 = xb - xa,
std::abs(d__1)) * -0.5
f;
const int mode_
returns the limits on the non-zero (mode=0) or normalized region (mode=1)
void sincosint(float x, float &sint, float &cint)
Abs< T >::type abs(const T &t)
int dzero(float a, float b, float &x0, float &rv, float eps, int mxf, F func)
Private version of the exponential integral.
VVIObjF(float kappa=0.01, float beta2=1., int mode=0)
Constructor.
float expint(float x)
Private version of the cosine and sine integral.
static const G4double kappa
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
void limits(float &xl, float &xu) const
density (mode=0) or distribution (mode=1) function
int sicif(float xx, float &si, float &ci)