An open source method of characteristics neutron transport code.
exponentials.h
Go to the documentation of this file.
1 
8 #ifndef EXPONENTIALS_H_
9 #define EXPONENTIALS_H_
10 
19 inline void expF1_poly(FP_PRECISION x, FP_PRECISION* expF1) {
20 
21  const FP_PRECISION p0 = 1.;
22  const FP_PRECISION p1 = -4.998618823537523 * 1E-1;
23  const FP_PRECISION p2 = 1.660264339632089 * 1E-1;
24  const FP_PRECISION p3 = -4.0607961247327 * 1E-2;
25  const FP_PRECISION p4 = 7.459558151235148 * 1E-3;
26  const FP_PRECISION p5 = -9.773063624328772 * 1E-4;
27  const FP_PRECISION p6 = 8.004982165323072 * 1E-5;
28  const FP_PRECISION p7 = -3.016010858852142 * 1E-6;
29 
30  *expF1 = p7*x + p6;
31  *expF1 = *expF1*x + p5;
32  *expF1 = *expF1*x + p4;
33  *expF1 = *expF1*x + p3;
34  *expF1 = *expF1*x + p2;
35  *expF1 = *expF1*x + p1;
36  *expF1 = *expF1*x + p0;
37 }
38 
39 
48 inline void expF2_poly(FP_PRECISION x, FP_PRECISION* expF2) {
49 
50  const FP_PRECISION p1 = 1.666656766985825 * 1E-1;
51  const FP_PRECISION p2 = -8.331262914972137 * 1E-2;
52  const FP_PRECISION p3 = 2.492325839109710 * 1E-2;
53  const FP_PRECISION p4 = -5.440953156443790 * 1E-3;
54  const FP_PRECISION p5 = 9.034802244154581 * 1E-4;
55  const FP_PRECISION p6 = -1.091608727341864 * 1E-4;
56  const FP_PRECISION p7 = 8.411465095972204 * 1E-6;
57  const FP_PRECISION p8 = -3.029020287833148 * 1E-7;
58 
59  *expF2 = p8*x + p7;
60  *expF2 = *expF2*x + p6;
61  *expF2 = *expF2*x + p5;
62  *expF2 = *expF2*x + p4;
63  *expF2 = *expF2*x + p3;
64  *expF2 = *expF2*x + p2;
65  *expF2 = *expF2*x + p1;
66  *expF2 = *expF2*x;
67 }
68 
69 
78 inline void expH_poly(FP_PRECISION x, FP_PRECISION* expH) {
79 
80  const FP_PRECISION p0 = 0.5;
81  const FP_PRECISION p1 = -3.33307480097059 * 1E-1;
82  const FP_PRECISION p2 = 1.248597190564637 * 1E-1;
83  const FP_PRECISION p3 = -3.305711771115656 * 1E-2;
84  const FP_PRECISION p4 = 6.667701492411682 * 1E-3;
85  const FP_PRECISION p5 = -1.028726890908856 * 1E-3;
86  const FP_PRECISION p6 = 1.145870114989106 * 1E-4;
87  const FP_PRECISION p7 = -8.079532805720403 * 1E-6;
88  const FP_PRECISION p8 = 2.657707693467421 * 1E-7;
89 
90  *expH = p8*x + p7;
91  *expH = *expH*x + p6;
92  *expH = *expH*x + p5;
93  *expH = *expH*x + p4;
94  *expH = *expH*x + p3;
95  *expH = *expH*x + p2;
96  *expH = *expH*x + p1;
97  *expH = *expH*x + p0;
98 }
99 
100 
110 inline void expG_fractional(FP_PRECISION x, FP_PRECISION* expG) {
111 
112  /* Coefficients for numerator */
113  const FP_PRECISION p0 = 0.5;
114  const FP_PRECISION p1 = 1.76558112351595 * 1E-1;
115  const FP_PRECISION p2 = 4.041584305811143 * 1E-2;
116  const FP_PRECISION p3 = 6.178333902037397 * 1E-3;
117  const FP_PRECISION p4 = 6.429894635552992 * 1E-4;
118  const FP_PRECISION p5 = 6.064409107557148 * 1E-5;
119 
120  /* Coefficients for denominator */
121  const FP_PRECISION d0 = 1.0;
122  const FP_PRECISION d1 = 6.864462055546078 * 1E-1;
123  const FP_PRECISION d2 = 2.263358514260129 * 1E-1;
124  const FP_PRECISION d3 = 4.721469893686252 * 1E-2;
125  const FP_PRECISION d4 = 6.883236664917246 * 1E-3;
126  const FP_PRECISION d5 = 7.036272419147752 * 1E-4;
127  const FP_PRECISION d6 = 6.064409107557148 * 1E-5;
128 
129  FP_PRECISION num, den;
130  den = d6*x + d5;
131  den = den*x + d4;
132  den = den*x + d3;
133  den = den*x + d2;
134  den = den*x + d1;
135  den = den*x + d0;
136  den = 1.f/den;
137 
138  num = p5*x + p4;
139  num = num*x + p3;
140  num = num*x + p2;
141  num = num*x + p1;
142  num = num*x + p0;
143 
144  *expG = num*den;
145 }
146 
147 
156 inline void expF1_fractional(FP_PRECISION x, FP_PRECISION* expF1) {
157 
158  /* Coefficients for numerator */
159  const FP_PRECISION p0 = 1.0;
160  const FP_PRECISION p1 = 2.4172687328033081 * 1E-1;
161  const FP_PRECISION p2 = 6.2804790965268531 * 1E-2;
162  const FP_PRECISION p3 = 1.0567595009016521 * 1E-2;
163  const FP_PRECISION p4 = 1.0059468082903561 * 1E-3;
164  const FP_PRECISION p5 = 1.9309063097411041 * 1E-4;
165 
166  /* Coefficients for denominator */
167  const FP_PRECISION d0 = 1.0;
168  const FP_PRECISION d1 = 7.4169266112320541 * 1E-1;
169  const FP_PRECISION d2 = 2.6722515319494311 * 1E-1;
170  const FP_PRECISION d3 = 6.1643725066901411 * 1E-2;
171  const FP_PRECISION d4 = 1.0590759992367811 * 1E-2;
172  const FP_PRECISION d5 = 1.0057980007137651 * 1E-3;
173  const FP_PRECISION d6 = 1.9309063097411041 * 1E-4;
174 
175  FP_PRECISION num, den;
176 
177  den = d6*x + d5;
178  den = den*x + d4;
179  den = den*x + d3;
180  den = den*x + d2;
181  den = den*x + d1;
182  den = den*x + d0;
183  den = 1.f / den;
184 
185  num = p5*x + p4;
186  num = num*x + p3;
187  num = num*x + p2;
188  num = num*x + p1;
189  num = num*x + p0;
190 
191  *expF1 = num*den;
192 }
193 
194 
203 inline void expF2_fractional(FP_PRECISION x, FP_PRECISION* expF2) {
204 
205  /* Coefficients for numerator */
206  const FP_PRECISION p1 = 1.666661470036759 * 1E-1;
207  const FP_PRECISION p2 = 3.59041632356318 * 1E-2;
208  const FP_PRECISION p3 = 7.675127136944033 * 1E-3;
209  const FP_PRECISION p4 = 6.408491755085618 * 1E-4;
210  const FP_PRECISION p5 = 1.367575707015872 * 1E-4;
211 
212  /* Coefficients for denominator */
213  const FP_PRECISION d0 = 1.0;
214  const FP_PRECISION d1 = 7.153333128932897 * 1E-1;
215  const FP_PRECISION d2 = 2.541555663123697 * 1E-1;
216  const FP_PRECISION d3 = 5.613392571426973 * 1E-2;
217  const FP_PRECISION d4 = 9.476002327852898 * 1E-3;
218  const FP_PRECISION d5 = 9.145637477815584 * 1E-4;
219  const FP_PRECISION d6 = 1.367575707015872 * 1E-4;
220 
221  FP_PRECISION num, den;
222  num = p5*x + p4;
223  num = num*x + p3;
224  num = num*x + p2;
225  num = num*x + p1;
226  num = num*x;
227 
228  den = d6*x + d5;
229  den = den*x + d4;
230  den = den*x + d3;
231  den = den*x + d2;
232  den = den*x + d1;
233  den = den*x + d0;
234  *expF2 = num/den;
235 }
236 
237 
246 inline void expH_fractional(FP_PRECISION x, FP_PRECISION* expH) {
247 
248  /* Coefficients for numerator */
249  const FP_PRECISION p0 = 0.5;
250  const FP_PRECISION p1 = 5.599412483229184 * 1E-2;
251  const FP_PRECISION p2 = 1.294939509305754 * 1E-2;
252  const FP_PRECISION p3 = 2.341166644220405 * 1E-3;
253  const FP_PRECISION p4 = 3.686858969421769 * 1E-5;
254  const FP_PRECISION p5 = 4.220477028150503 * 1E-5;
255 
256  /* Coefficients for denominator */
257  const FP_PRECISION d0 = 1.0;
258  const FP_PRECISION d1 = 7.787274561075199 * 1E-1;
259  const FP_PRECISION d2 = 2.945145030273455 * 1E-1;
260  const FP_PRECISION d3 = 7.440380752801196 * 1E-2;
261  const FP_PRECISION d4 = 1.220791761275212 * 1E-2;
262  const FP_PRECISION d5 = 2.354181374425252 * 1E-3;
263  const FP_PRECISION d6 = 3.679462493221416 * 1E-5;
264  const FP_PRECISION d7 = 4.220477028150503 * 1E-5;
265 
266  FP_PRECISION num, den;
267  num = p5*x + p4;
268  num = num*x + p3;
269  num = num*x + p2;
270  num = num*x + p1;
271  num = num*x + p0;
272 
273  den = d7*x + d6;
274  den = den*x + d5;
275  den = den*x + d4;
276  den = den*x + d3;
277  den = den*x + d2;
278  den = den*x + d1;
279  den = den*x + d0;
280  *expH = num/den;
281 }
282 
283 
293 inline void expG2_fractional(FP_PRECISION x, FP_PRECISION* expG2) {
294 
295  /* Coefficients for numerator */
296  const FP_PRECISION a1 = -8.335775885589858 * 1E-2;
297  const FP_PRECISION a2 = -3.603942303847604 * 1E-3;
298  const FP_PRECISION a3 = 3.7673183263550827 * 1E-3;
299  const FP_PRECISION a4 = 1.124183494990467 * 1E-5;
300  const FP_PRECISION a5 = 1.6837426505799449 * 1E-4;
301 
302  /* Coefficients for denominator */
303  const FP_PRECISION b1 = 7.454048371823628 * 1E-1;
304  const FP_PRECISION b2 = 2.3794300531408347 * 1E-1;
305  const FP_PRECISION b3= 5.367250964303789 * 1E-2;
306  const FP_PRECISION b4 = 6.125197988351906 * 1E-3;
307  const FP_PRECISION b5 = 1.0102514456857377 * 1E-3;
308 
309  FP_PRECISION num, den;
310  num = a5 * x + a4;
311  num = num * x + a3;
312  num = num * x + a2;
313  num = num * x + a1;
314  num *= x;
315 
316  den = b5 * x + b4;
317  den = den * x + b3;
318  den = den * x + b2;
319  den = den * x + b1;
320  den = den * x + 1.f;
321 
322  *expG2 = num / den;
323 }
324 
325 
335 inline void cram7(FP_PRECISION x, FP_PRECISION* expv) {
336 
337  /* Generated in Mathematica, accurate to 6.18 digits (single precision),
338  tau [-1.5e6,0] */
339  const FP_PRECISION c1n = -1.00000014302666667201396424463;
340  const FP_PRECISION c2n = 2.34841040052684510704433796447 * 1E-1;
341  const FP_PRECISION c3n = -6.24785939603762121316592924635 * 1E-2;
342  const FP_PRECISION c4n = 1.00434102711342948752684759736 * 1E-2;
343  const FP_PRECISION c5n = -1.35724435934263932676353754751 * 1E-3;
344  const FP_PRECISION c6n = 9.51474224366003625378414851577 * 1E-5;
345  const FP_PRECISION c7n = -1.60076055315534285575266516209 * 1E-5;
346 
347  const FP_PRECISION c0d = 1.;
348  const FP_PRECISION c1d = -7.34847118148952339633322706422 * 1E-1;
349  const FP_PRECISION c2d = 2.63193362386411901729092564316 * 1E-1;
350  const FP_PRECISION c3d = -6.09467155163113059870970359654 * 1E-2;
351  const FP_PRECISION c4d = 1.00863490579686697359577926719 * 1E-2;
352  const FP_PRECISION c5d = -1.35667018708833025497446407598 * 1E-3;
353  const FP_PRECISION c6d = 9.51502816434275317085698085885 * 1E-5;
354  const FP_PRECISION c7d = -1.60076032420105715765981718742 * 1E-5;
355 
356  FP_PRECISION num, den;
357 
358  den = c7d;
359  den = den * x + c6d;
360  den = den * x + c5d;
361  den = den * x + c4d;
362  den = den * x + c3d;
363  den = den * x + c2d;
364  den = den * x + c1d;
365  den = den * x + c0d;
366 
367  num = c7n;
368  num = num * x + c6n;
369  num = num * x + c5n;
370  num = num * x + c4n;
371  num = num * x + c3n;
372  num = num * x + c2n;
373  num = num * x + c1n;
374  num = num * x;
375 
376  *expv = num / den;
377 }
378 
379 
380 #endif // EXPONENTIALS_H_
void expF1_fractional(FP_PRECISION x, FP_PRECISION *expF1)
Computes the F1 exponential term.
Definition: exponentials.h:156
void cram7(FP_PRECISION x, FP_PRECISION *expv)
Evaluates the function 1 - exp(x) for x negative. Vectorises well. Accurate to 6.18 digits (single pr...
Definition: exponentials.h:335
void expF2_fractional(FP_PRECISION x, FP_PRECISION *expF2)
Computes the F2 exponential term.
Definition: exponentials.h:203
void expH_poly(FP_PRECISION x, FP_PRECISION *expH)
Computes the H exponential term.
Definition: exponentials.h:78
void expF1_poly(FP_PRECISION x, FP_PRECISION *expF1)
Computes the F1 exponential term.
Definition: exponentials.h:19
void expF2_poly(FP_PRECISION x, FP_PRECISION *expF2)
Computes the F2 exponential term.
Definition: exponentials.h:48
void expG2_fractional(FP_PRECISION x, FP_PRECISION *expG2)
Computes the G2 exponential term.
Definition: exponentials.h:293
void expH_fractional(FP_PRECISION x, FP_PRECISION *expH)
Computes the H exponential term.
Definition: exponentials.h:246
void expG_fractional(FP_PRECISION x, FP_PRECISION *expG)
Computes an exponential G term, used to reconstruct F1, F2 and H.
Definition: exponentials.h:110