19072640351

Committed 04 Nov 2025 02:49PM UTC coverage: 81.932% (+0.06%) from 81.872%

Build # 19072640351

Build Type

Pull #3413

github

Committed by

web-flow

Commit Message

Merge 043dbf5f9 into bd76fc056

Pull Request Pull Request #3413: More interpolation types in Tabular.

Coverage Stats

16712 of 23260 branches covered (71.85%)

Branch coverage included in aggregate %.

72 of 129 new or added lines in 3 files covered. (55.81%)

125 existing lines in 2 files now uncovered.

54241 of 63340 relevant lines covered (85.63%)

44677982.28 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

94.78

/src/math_functions.cpp

1	#include "openmc/math_functions.h"
2
3	#include "openmc/external/Faddeeva.hh"
4	#include "openmc/external/LambertW.hpp"
5
6	#include "openmc/constants.h"
7	#include "openmc/random_lcg.h"
8
9	namespace openmc {
10
11	//==============================================================================
12	// Mathematical methods
13	//==============================================================================
14
15	double normal_percentile(double p)	168 ✔
16	{
17	constexpr double p_low = 0.02425;	168 ✔
18	constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
19	-2.759285104469687e2, 1.383577518672690e2, -3.066479806614716e1,
20	2.506628277459239e0};
21	constexpr double b[5] = {-5.447609879822406e1, 1.615858368580409e2,
22	-1.556989798598866e2, 6.680131188771972e1, -1.328068155288572e1};
23	constexpr double c[6] = {-7.784894002430293e-3, -3.223964580411365e-1,
24	-2.400758277161838, -2.549732539343734, 4.374664141464968,
25	2.938163982698783};
26	constexpr double d[4] = {7.784695709041462e-3, 3.224671290700398e-1,
27	2.445134137142996, 3.754408661907416};
28
29	// The rational approximation used here is from an unpublished work at
30	// http://home.online.no/~pjacklam/notes/invnorm/
31
32	double z;
33	double q;
34
35	if (p < p_low) {	168 ✔
36	// Rational approximation for lower region.
37
38	q = std::sqrt(-2.0 * std::log(p));	11 ✔
39	z = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /	11 ✔
40	((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);	11 ✔
41
42	} else if (p <= 1.0 - p_low) {	157 ✔
43	// Rational approximation for central region
44	q = p - 0.5;	146 ✔
45	double r = q * q;	146 ✔
46	z = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) *	146 ✔
47	q /
48	(((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1.0);	146 ✔
49
50	} else {
51	// Rational approximation for upper region
52
53	q = std::sqrt(-2.0 * std::log(1.0 - p));	11 ✔
54	z = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) /	11 ✔
55	((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1.0);	11 ✔
56	}
57
58	// Refinement based on Newton's method
59
60	z = z - (0.5 * std::erfc(-z / std::sqrt(2.0)) - p) * std::sqrt(2.0 * PI) *	168 ✔
61	std::exp(0.5 * z * z);	168 ✔
62
63	return z;	168 ✔
64	}
65
66	double t_percentile(double p, int df)	310 ✔
67	{
68	double t;
69
70	if (df == 1) {	310 ✔
71	// For one degree of freedom, the t-distribution becomes a Cauchy
72	// distribution whose cdf we can invert directly
73
74	t = std::tan(PI * (p - 0.5));	71 ✔
75	} else if (df == 2) {	239 ✔
76	// For two degrees of freedom, the cdf is given by 1/2 + x/(2*sqrt(x^2 +
77	// 2)). This can be directly inverted to yield the solution below
78
79	t = 2.0 * std::sqrt(2.0) * (p - 0.5) /	71 ✔
80	std::sqrt(1. - 4. * std::pow(p - 0.5, 2.));	71 ✔
81	} else {
82	// This approximation is from E. Olusegun George and Meenakshi Sivaram, "A
83	// modification of the Fisher-Cornish approximation for the student t
84	// percentiles," Communication in Statistics - Simulation and Computation,
85	// 16 (4), pp. 1123-1132 (1987).
86	double n = df;	168 ✔
87	double k = 1. / (n - 2.);	168 ✔
88	double z = normal_percentile(p);	168 ✔
89	double z2 = z * z;	168 ✔
90	t = std::sqrt(n * k) *	168 ✔
91	(z + (z2 - 3.) * z * k / 4. +	168 ✔
92	((5. * z2 - 56.) * z2 + 75.) * z * k * k / 96. +	168 ✔
93	(((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) * z * k * k * k / 384.);	168 ✔
94	}
95
96	return t;	310 ✔
97	}
98
99	void calc_pn_c(int n, double x, double pnx[])	222,916,793 ✔
100	{
101	pnx[0] = 1.;	222,916,793 ✔
102	if (n >= 1) {	222,916,793 ✔
103	pnx[1] = x;	154,709,494 ✔
104	}
105
106	// Use recursion relation to build the higher orders
107	for (int l = 1; l < n; l++) {	228,490,130 ✔
108	pnx[l + 1] = ((2 * l + 1) * x * pnx[l] - l * pnx[l - 1]) / (l + 1);	5,573,337 ✔
109	}
110	}	222,916,793 ✔
111
112	double evaluate_legendre(int n, const double data[], double x)	96,380,554 ✔
113	{
114	double* pnx = new double[n + 1];	96,380,554 !
115	double val = 0.0;	96,380,554 ✔
116	calc_pn_c(n, x, pnx);	96,380,554 ✔
117	for (int l = 0; l <= n; l++) {	241,994,875 ✔
118	val += (l + 0.5) * data[l] * pnx[l];	145,614,321 ✔
119	}
120	delete[] pnx;	96,380,554 !
121	return val;	96,380,554 ✔
122	}
123
124	void calc_rn_c(int n, const double uvw[3], double rn[])	11 ✔
125	{
126	Direction u {uvw};	11 ✔
127	calc_rn(n, u, rn);	11 ✔
128	}	11 ✔
129
130	void calc_rn(int n, Direction u, double rn[])	4,568,344 ✔
131	{
132	// rn[] is assumed to have already been allocated to the correct size
133
134	// Store the cosine of the polar angle and the azimuthal angle
135	double w = u.z;	4,568,344 ✔
136	double phi;
137	if (u.x == 0.) {	4,568,344 !
138	phi = 0.;	×
139	} else {
140	phi = std::atan2(u.y, u.x);	4,568,344 ✔
141	}
142
143	// Store the shorthand of 1-w * w
144	double w2m1 = 1. - w * w;	4,568,344 ✔
145
146	// Now evaluate the spherical harmonics function
147	rn[0] = 1.;	4,568,344 ✔
148	int i = 0;	4,568,344 ✔
149	for (int l = 1; l <= n; l++) {	22,844,778 ✔
150	// Set the index to the start of this order
151	i += 2 * (l - 1) + 1;	18,276,434 ✔
152
153	// Now evaluate each
154	switch (l) {	18,276,434 !
155	case 1:	4,568,344 ✔
156	// l = 1, m = -1
157	rn[i] = -(std::sqrt(w2m1) * std::sin(phi));	4,568,344 ✔
158	// l = 1, m = 0
159	rn[i + 1] = w;	4,568,344 ✔
160	// l = 1, m = 1
161	rn[i + 2] = -(std::sqrt(w2m1) * std::cos(phi));	4,568,344 ✔
162	break;	4,568,344 ✔
163	case 2:	4,568,344 ✔
164	// l = 2, m = -2
165	rn[i] = 0.288675134594813 * (-3. * w * w + 3.) * std::sin(2. * phi);	4,568,344 ✔
166	// l = 2, m = -1
167	rn[i + 1] = -(1.73205080756888 * w * std::sqrt(w2m1) * std::sin(phi));	4,568,344 ✔
168	// l = 2, m = 0
169	rn[i + 2] = 1.5 * w * w - 0.5;	4,568,344 ✔
170	// l = 2, m = 1
171	rn[i + 3] = -(1.73205080756888 * w * std::sqrt(w2m1) * std::cos(phi));	4,568,344 ✔
172	// l = 2, m = 2
173	rn[i + 4] = 0.288675134594813 * (-3. * w * w + 3.) * std::cos(2. * phi);	4,568,344 ✔
174	break;	4,568,344 ✔
175	case 3:	4,568,344 ✔
176	// l = 3, m = -3
177	rn[i] = -(0.790569415042095 * std::pow(w2m1, 1.5) * std::sin(3. * phi));	4,568,344 ✔
178	// l = 3, m = -2
179	rn[i + 1] = 1.93649167310371 * w * (w2m1)std::sin(2. phi);	4,568,344 ✔
180	// l = 3, m = -1
181	rn[i + 2] = -(0.408248290463863 * std::sqrt(w2m1) *	4,568,344 ✔
182	((7.5) * w * w - 3. / 2.) * std::sin(phi));	4,568,344 ✔
183	// l = 3, m = 0
184	rn[i + 3] = 2.5 * std::pow(w, 3) - 1.5 * w;	4,568,344 ✔
185	// l = 3, m = 1
186	rn[i + 4] = -(0.408248290463863 * std::sqrt(w2m1) *	4,568,344 ✔
187	((7.5) * w * w - 3. / 2.) * std::cos(phi));	4,568,344 ✔
188	// l = 3, m = 2
189	rn[i + 5] = 1.93649167310371 * w * (w2m1)std::cos(2. phi);	4,568,344 ✔
190	// l = 3, m = 3
191	rn[i + 6] =	4,568,344 ✔
192	-(0.790569415042095 * std::pow(w2m1, 1.5) * std::cos(3. * phi));	4,568,344 ✔
193	break;	4,568,344 ✔
194	case 4:	4,568,344 ✔
195	// l = 4, m = -4
196	rn[i] = 0.739509972887452 * (w2m1 * w2m1) * std::sin(4.0 * phi);	4,568,344 ✔
197	// l = 4, m = -3
198	rn[i + 1] =	4,568,344 ✔
199	-(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::sin(3. * phi));	4,568,344 ✔
200	// l = 4, m = -2
201	rn[i + 2] =	4,568,344 ✔
202	0.074535599249993 * (w2m1) * (52.5 * w * w - 7.5) * std::sin(2. * phi);	4,568,344 ✔
203	// l = 4, m = -1
204	rn[i + 3] = -(0.316227766016838 * std::sqrt(w2m1) *	9,136,688 ✔
205	(17.5 * std::pow(w, 3) - 7.5 * w) * std::sin(phi));	4,568,344 ✔
206	// l = 4, m = 0
207	rn[i + 4] = 4.375 * std::pow(w, 4) - 3.75 * w * w + 0.375;	4,568,344 ✔
208	// l = 4, m = 1
209	rn[i + 5] = -(0.316227766016838 * std::sqrt(w2m1) *	9,136,688 ✔
210	(17.5 * std::pow(w, 3) - 7.5 * w) * std::cos(phi));	4,568,344 ✔
211	// l = 4, m = 2
212	rn[i + 6] =	4,568,344 ✔
213	0.074535599249993 * (w2m1) * (52.5 * w * w - 7.5) * std::cos(2. * phi);	4,568,344 ✔
214	// l = 4, m = 3
215	rn[i + 7] =	4,568,344 ✔
216	-(2.09165006633519 * w * std::pow(w2m1, 1.5) * std::cos(3. * phi));	4,568,344 ✔
217	// l = 4, m = 4
218	rn[i + 8] = 0.739509972887452 * w2m1 * w2m1 * std::cos(4.0 * phi);	4,568,344 ✔
219	break;	4,568,344 ✔
220	case 5:	3,003 ✔
221	// l = 5, m = -5
222	rn[i] = -(0.701560760020114 * std::pow(w2m1, 2.5) * std::sin(5.0 * phi));	3,003 ✔
223	// l = 5, m = -4
224	rn[i + 1] = 2.21852991866236 * w * w2m1 * w2m1 * std::sin(4.0 * phi);	3,003 ✔
225	// l = 5, m = -3
226	rn[i + 2] = -(0.00996023841111995 * std::pow(w2m1, 1.5) *	3,003 ✔
227	((945.0 / 2.) * w * w - 52.5) * std::sin(3. * phi));	3,003 ✔
228	// l = 5, m = -2
229	rn[i + 3] = 0.0487950036474267 * (w2m1) *	6,006 ✔
230	((315.0 / 2.) * std::pow(w, 3) - 52.5 * w) *	3,003 ✔
231	std::sin(2. * phi);	3,003 ✔
232	// l = 5, m = -1
233	rn[i + 4] =	3,003 ✔
234	-(0.258198889747161 * std::sqrt(w2m1) *	6,006 ✔
235	(39.375 * std::pow(w, 4) - 105.0 / 4.0 * w * w + 15.0 / 8.0) *	3,003 ✔
236	std::sin(phi));	3,003 ✔
237	// l = 5, m = 0
238	rn[i + 5] = 7.875 * std::pow(w, 5) - 8.75 * std::pow(w, 3) + 1.875 * w;	3,003 ✔
239	// l = 5, m = 1
240	rn[i + 6] =	3,003 ✔
241	-(0.258198889747161 * std::sqrt(w2m1) *	6,006 ✔
242	(39.375 * std::pow(w, 4) - 105.0 / 4.0 * w * w + 15.0 / 8.0) *	3,003 ✔
243	std::cos(phi));	3,003 ✔
244	// l = 5, m = 2
245	rn[i + 7] = 0.0487950036474267 * (w2m1) *	6,006 ✔
246	((315.0 / 2.) * std::pow(w, 3) - 52.5 * w) *	3,003 ✔
247	std::cos(2. * phi);	3,003 ✔
248	// l = 5, m = 3
249	rn[i + 8] = -(0.00996023841111995 * std::pow(w2m1, 1.5) *	3,003 ✔
250	((945.0 / 2.) * w * w - 52.5) * std::cos(3. * phi));	3,003 ✔
251	// l = 5, m = 4
252	rn[i + 9] = 2.21852991866236 * w * w2m1 * w2m1 * std::cos(4.0 * phi);	3,003 ✔
253	// l = 5, m = 5
254	rn[i + 10] =	3,003 ✔
255	-(0.701560760020114 * std::pow(w2m1, 2.5) * std::cos(5.0 * phi));	3,003 ✔
256	break;	3,003 ✔
257	case 6:	11 ✔
258	// l = 6, m = -6
259	rn[i] = 0.671693289381396 * std::pow(w2m1, 3) * std::sin(6.0 * phi);	11 ✔
260	// l = 6, m = -5
261	rn[i + 1] =	11 ✔
262	-(2.32681380862329 * w * std::pow(w2m1, 2.5) * std::sin(5.0 * phi));	11 ✔
263	// l = 6, m = -4
264	rn[i + 2] = 0.00104990131391452 * w2m1 * w2m1 *	11 ✔
265	((10395.0 / 2.) * w * w - 945.0 / 2.) * std::sin(4.0 * phi);	11 ✔
266	// l = 6, m = -3
267	rn[i + 3] = -(0.00575054632785295 * std::pow(w2m1, 1.5) *	22 ✔
268	((3465.0 / 2.) * std::pow(w, 3) - 945.0 / 2. * w) *	11 ✔
269	std::sin(3. * phi));	11 ✔
270	// l = 6, m = -2
271	rn[i + 4] =	11 ✔
272	0.0345032779671177 * (w2m1) *	22 ✔
273	((3465.0 / 8.0) * std::pow(w, 4) - 945.0 / 4.0 * w * w + 105.0 / 8.0) *	11 ✔
274	std::sin(2. * phi);	11 ✔
275	// l = 6, m = -1
276	rn[i + 5] = -(0.218217890235992 * std::sqrt(w2m1) *	22 ✔
277	((693.0 / 8.0) * std::pow(w, 5) -	11 ✔
278	315.0 / 4.0 * std::pow(w, 3) + (105.0 / 8.0) * w) *	11 ✔
279	std::sin(phi));	11 ✔
280	// l = 6, m = 0
281	rn[i + 6] = 14.4375 * std::pow(w, 6) - 19.6875 * std::pow(w, 4) +	11 ✔
282	6.5625 * w * w - 0.3125;	11 ✔
283	// l = 6, m = 1
284	rn[i + 7] = -(0.218217890235992 * std::sqrt(w2m1) *	22 ✔
285	((693.0 / 8.0) * std::pow(w, 5) -	11 ✔
286	315.0 / 4.0 * std::pow(w, 3) + (105.0 / 8.0) * w) *	11 ✔
287	std::cos(phi));	11 ✔
288	// l = 6, m = 2
289	rn[i + 8] =	11 ✔
290	0.0345032779671177 * w2m1 *	22 ✔
291	((3465.0 / 8.0) * std::pow(w, 4) - 945.0 / 4.0 * w * w + 105.0 / 8.0) *	11 ✔
292	std::cos(2. * phi);	11 ✔
293	// l = 6, m = 3
294	rn[i + 9] = -(0.00575054632785295 * std::pow(w2m1, 1.5) *	22 ✔
295	((3465.0 / 2.) * std::pow(w, 3) - 945.0 / 2. * w) *	11 ✔
296	std::cos(3. * phi));	11 ✔
297	// l = 6, m = 4
298	rn[i + 10] = 0.00104990131391452 * w2m1 * w2m1 *	11 ✔
299	((10395.0 / 2.) * w * w - 945.0 / 2.) * std::cos(4.0 * phi);	11 ✔
300	// l = 6, m = 5
301	rn[i + 11] =	11 ✔
302	-(2.32681380862329 * w * std::pow(w2m1, 2.5) * std::cos(5.0 * phi));	11 ✔
303	// l = 6, m = 6
304	rn[i + 12] = 0.671693289381396 * std::pow(w2m1, 3) * std::cos(6.0 * phi);	11 ✔
305	break;	11 ✔
306	case 7:	11 ✔
307	// l = 7, m = -7
308	rn[i] = -(0.647259849287749 * std::pow(w2m1, 3.5) * std::sin(7.0 * phi));	11 ✔
309	// l = 7, m = -6
310	rn[i + 1] =	11 ✔
311	2.42182459624969 * w * std::pow(w2m1, 3) * std::sin(6.0 * phi);	11 ✔
312	// l = 7, m = -5
313	rn[i + 2] =	11 ✔
314	-(9.13821798555235e-5 * std::pow(w2m1, 2.5) *	11 ✔
315	((135135.0 / 2.) * w * w - 10395.0 / 2.) * std::sin(5.0 * phi));	11 ✔
316	// l = 7, m = -4
317	rn[i + 3] = 0.000548293079133141 * w2m1 * w2m1 *	22 ✔
318	((45045.0 / 2.) * std::pow(w, 3) - 10395.0 / 2. * w) *	11 ✔
319	std::sin(4.0 * phi);	11 ✔
320	// l = 7, m = -3
321	rn[i + 4] = -(0.00363696483726654 * std::pow(w2m1, 1.5) *	22 ✔
322	((45045.0 / 8.0) * std::pow(w, 4) - 10395.0 / 4.0 * w * w +	11 ✔
323	945.0 / 8.0) *	11 ✔
324	std::sin(3. * phi));	11 ✔
325	// l = 7, m = -2
326	rn[i + 5] = 0.025717224993682 * (w2m1) *	22 ✔
327	((9009.0 / 8.0) * std::pow(w, 5) -	11 ✔
328	3465.0 / 4.0 * std::pow(w, 3) + (945.0 / 8.0) * w) *	11 ✔
329	std::sin(2. * phi);	11 ✔
330	// l = 7, m = -1
331	rn[i + 6] =	11 ✔
332	-(0.188982236504614 * std::sqrt(w2m1) *	22 ✔
333	((3003.0 / 16.0) * std::pow(w, 6) - 3465.0 / 16.0 * std::pow(w, 4) +	11 ✔
334	(945.0 / 16.0) * w * w - 35.0 / 16.0) *	11 ✔
335	std::sin(phi));	11 ✔
336	// l = 7, m = 0
337	rn[i + 7] = 26.8125 * std::pow(w, 7) - 43.3125 * std::pow(w, 5) +	11 ✔
338	19.6875 * std::pow(w, 3) - 2.1875 * w;	11 ✔
339	// l = 7, m = 1
340	rn[i + 8] =	11 ✔
341	-(0.188982236504614 * std::sqrt(w2m1) *	22 ✔
342	((3003.0 / 16.0) * std::pow(w, 6) - 3465.0 / 16.0 * std::pow(w, 4) +	11 ✔
343	(945.0 / 16.0) * w * w - 35.0 / 16.0) *	11 ✔
344	std::cos(phi));	11 ✔
345	// l = 7, m = 2
346	rn[i + 9] = 0.025717224993682 * (w2m1) *	22 ✔
347	((9009.0 / 8.0) * std::pow(w, 5) -	11 ✔
348	3465.0 / 4.0 * std::pow(w, 3) + (945.0 / 8.0) * w) *	11 ✔
349	std::cos(2. * phi);	11 ✔
350	// l = 7, m = 3
351	rn[i + 10] = -(0.00363696483726654 * std::pow(w2m1, 1.5) *	22 ✔
352	((45045.0 / 8.0) * std::pow(w, 4) - 10395.0 / 4.0 * w * w +	11 ✔
353	945.0 / 8.0) *	11 ✔
354	std::cos(3. * phi));	11 ✔
355	// l = 7, m = 4
356	rn[i + 11] = 0.000548293079133141 * w2m1 * w2m1 *	22 ✔
357	((45045.0 / 2.) * std::pow(w, 3) - 10395.0 / 2. * w) *	11 ✔
358	std::cos(4.0 * phi);	11 ✔
359	// l = 7, m = 5
360	rn[i + 12] =	11 ✔
361	-(9.13821798555235e-5 * std::pow(w2m1, 2.5) *	11 ✔
362	((135135.0 / 2.) * w * w - 10395.0 / 2.) * std::cos(5.0 * phi));	11 ✔
363	// l = 7, m = 6
364	rn[i + 13] =	11 ✔
365	2.42182459624969 * w * std::pow(w2m1, 3) * std::cos(6.0 * phi);	11 ✔
366	// l = 7, m = 7
367	rn[i + 14] =	11 ✔
368	-(0.647259849287749 * std::pow(w2m1, 3.5) * std::cos(7.0 * phi));	11 ✔
369	break;	11 ✔
370	case 8:	11 ✔
371	// l = 8, m = -8
372	rn[i] = 0.626706654240044 * std::pow(w2m1, 4) * std::sin(8.0 * phi);	11 ✔
373	// l = 8, m = -7
374	rn[i + 1] =	11 ✔
375	-(2.50682661696018 * w * std::pow(w2m1, 3.5) * std::sin(7.0 * phi));	11 ✔
376	// l = 8, m = -6
377	rn[i + 2] = 6.77369783729086e-6 * std::pow(w2m1, 3) *	11 ✔
378	((2027025.0 / 2.) * w * w - 135135.0 / 2.) *	11 ✔
379	std::sin(6.0 * phi);	11 ✔
380	// l = 8, m = -5
381	rn[i + 3] = -(4.38985792528482e-5 * std::pow(w2m1, 2.5) *	22 ✔
382	((675675.0 / 2.) * std::pow(w, 3) - 135135.0 / 2. * w) *	11 ✔
383	std::sin(5.0 * phi));	11 ✔
384	// l = 8, m = -4
385	rn[i + 4] = 0.000316557156832328 * w2m1 * w2m1 *	22 ✔
386	((675675.0 / 8.0) * std::pow(w, 4) - 135135.0 / 4.0 * w * w +	11 ✔
387	10395.0 / 8.0) *	11 ✔
388	std::sin(4.0 * phi);	11 ✔
389	// l = 8, m = -3
390	rn[i + 5] = -(0.00245204119306875 * std::pow(w2m1, 1.5) *	22 ✔
391	((135135.0 / 8.0) * std::pow(w, 5) -	11 ✔
392	45045.0 / 4.0 * std::pow(w, 3) + (10395.0 / 8.0) * w) *	11 ✔
393	std::sin(3. * phi));	11 ✔
394	// l = 8, m = -2
395	rn[i + 6] =	11 ✔
396	0.0199204768222399 * (w2m1) *	22 ✔
397	((45045.0 / 16.0) * std::pow(w, 6) - 45045.0 / 16.0 * std::pow(w, 4) +	11 ✔
398	(10395.0 / 16.0) * w * w - 315.0 / 16.0) *	11 ✔
399	std::sin(2. * phi);	11 ✔
400	// l = 8, m = -1
401	rn[i + 7] =	11 ✔
402	-(0.166666666666667 * std::sqrt(w2m1) *	22 ✔
403	((6435.0 / 16.0) * std::pow(w, 7) - 9009.0 / 16.0 * std::pow(w, 5) +	11 ✔
404	(3465.0 / 16.0) * std::pow(w, 3) - 315.0 / 16.0 * w) *	11 ✔
405	std::sin(phi));	11 ✔
406	// l = 8, m = 0
407	rn[i + 8] = 50.2734375 * std::pow(w, 8) - 93.84375 * std::pow(w, 6) +	11 ✔
408	54.140625 * std::pow(w, 4) - 9.84375 * w * w + 0.2734375;	11 ✔
409	// l = 8, m = 1
410	rn[i + 9] =	11 ✔
411	-(0.166666666666667 * std::sqrt(w2m1) *	22 ✔
412	((6435.0 / 16.0) * std::pow(w, 7) - 9009.0 / 16.0 * std::pow(w, 5) +	11 ✔
413	(3465.0 / 16.0) * std::pow(w, 3) - 315.0 / 16.0 * w) *	11 ✔
414	std::cos(phi));	11 ✔
415	// l = 8, m = 2
416	rn[i + 10] =	11 ✔
417	0.0199204768222399 * (w2m1) *	22 ✔
418	((45045.0 / 16.0) * std::pow(w, 6) - 45045.0 / 16.0 * std::pow(w, 4) +	11 ✔
419	(10395.0 / 16.0) * w * w - 315.0 / 16.0) *	11 ✔
420	std::cos(2. * phi);	11 ✔
421	// l = 8, m = 3
422	rn[i + 11] = -(0.00245204119306875 * std::pow(w2m1, 1.5) *	22 ✔
423	((135135.0 / 8.0) * std::pow(w, 5) -	11 ✔
424	45045.0 / 4.0 * std::pow(w, 3) + (10395.0 / 8.0) * w) *	11 ✔
425	std::cos(3. * phi));	11 ✔
426	// l = 8, m = 4
427	rn[i + 12] = 0.000316557156832328 * w2m1 * w2m1 *	22 ✔
428	((675675.0 / 8.0) * std::pow(w, 4) - 135135.0 / 4.0 * w * w +	11 ✔
429	10395.0 / 8.0) *	11 ✔
430	std::cos(4.0 * phi);	11 ✔
431	// l = 8, m = 5
432	rn[i + 13] = -(4.38985792528482e-5 * std::pow(w2m1, 2.5) *	22 ✔
433	((675675.0 / 2.) * std::pow(w, 3) - 135135.0 / 2. * w) *	11 ✔
434	std::cos(5.0 * phi));	11 ✔
435	// l = 8, m = 6
436	rn[i + 14] = 6.77369783729086e-6 * std::pow(w2m1, 3) *	11 ✔
437	((2027025.0 / 2.) * w * w - 135135.0 / 2.) *	11 ✔
438	std::cos(6.0 * phi);	11 ✔
439	// l = 8, m = 7
440	rn[i + 15] =	11 ✔
441	-(2.50682661696018 * w * std::pow(w2m1, 3.5) * std::cos(7.0 * phi));	11 ✔
442	// l = 8, m = 8
443	rn[i + 16] = 0.626706654240044 * std::pow(w2m1, 4) * std::cos(8.0 * phi);	11 ✔
444	break;	11 ✔
445	case 9:	11 ✔
446	// l = 9, m = -9
447	rn[i] = -(0.609049392175524 * std::pow(w2m1, 4.5) * std::sin(9.0 * phi));	11 ✔
448	// l = 9, m = -8
449	rn[i + 1] =	11 ✔
450	2.58397773170915 * w * std::pow(w2m1, 4) * std::sin(8.0 * phi);	11 ✔
451	// l = 9, m = -7
452	rn[i + 2] =	11 ✔
453	-(4.37240315267812e-7 * std::pow(w2m1, 3.5) *	11 ✔
454	((34459425.0 / 2.) * w * w - 2027025.0 / 2.) * std::sin(7.0 * phi));	11 ✔
455	// l = 9, m = -6
456	rn[i + 3] = 3.02928976464514e-6 * std::pow(w2m1, 3) *	11 ✔
457	((11486475.0 / 2.) * std::pow(w, 3) - 2027025.0 / 2. * w) *	11 ✔
458	std::sin(6.0 * phi);	11 ✔
459	// l = 9, m = -5
460	rn[i + 4] = -(2.34647776186144e-5 * std::pow(w2m1, 2.5) *	22 ✔
461	((11486475.0 / 8.0) * std::pow(w, 4) -	11 ✔
462	2027025.0 / 4.0 * w * w + 135135.0 / 8.0) *	11 ✔
463	std::sin(5.0 * phi));	11 ✔
464	// l = 9, m = -4
465	rn[i + 5] = 0.000196320414650061 * w2m1 * w2m1 *	22 ✔
466	((2297295.0 / 8.0) * std::pow(w, 5) -	11 ✔
467	675675.0 / 4.0 * std::pow(w, 3) + (135135.0 / 8.0) * w) *	11 ✔
468	std::sin(4.0 * phi);	11 ✔
469	// l = 9, m = -3
470	rn[i + 6] = -(	11 ✔
471	0.00173385495536766 * std::pow(w2m1, 1.5) *	22 ✔
472	((765765.0 / 16.0) * std::pow(w, 6) - 675675.0 / 16.0 * std::pow(w, 4) +	11 ✔
473	(135135.0 / 16.0) * w * w - 3465.0 / 16.0) *	11 ✔
474	std::sin(3. * phi));	11 ✔
475	// l = 9, m = -2
476	rn[i + 7] =	11 ✔
477	0.0158910431540932 * (w2m1) *	22 ✔
478	((109395.0 / 16.0) * std::pow(w, 7) - 135135.0 / 16.0 * std::pow(w, 5) +	11 ✔
479	(45045.0 / 16.0) * std::pow(w, 3) - 3465.0 / 16.0 * w) *	11 ✔
480	std::sin(2. * phi);	11 ✔
481	// l = 9, m = -1
482	rn[i + 8] = -(	11 ✔
483	0.149071198499986 * std::sqrt(w2m1) *	22 ✔
484	((109395.0 / 128.0) * std::pow(w, 8) - 45045.0 / 32.0 * std::pow(w, 6) +	11 ✔
485	(45045.0 / 64.0) * std::pow(w, 4) - 3465.0 / 32.0 * w * w +	11 ✔
486	315.0 / 128.0) *	11 ✔
487	std::sin(phi));	11 ✔
488	// l = 9, m = 0
489	rn[i + 9] = 94.9609375 * std::pow(w, 9) - 201.09375 * std::pow(w, 7) +	11 ✔
490	140.765625 * std::pow(w, 5) - 36.09375 * std::pow(w, 3) +	11 ✔
491	2.4609375 * w;	11 ✔
492	// l = 9, m = 1
493	rn[i + 10] = -(	11 ✔
494	0.149071198499986 * std::sqrt(w2m1) *	22 ✔
495	((109395.0 / 128.0) * std::pow(w, 8) - 45045.0 / 32.0 * std::pow(w, 6) +	11 ✔
496	(45045.0 / 64.0) * std::pow(w, 4) - 3465.0 / 32.0 * w * w +	11 ✔
497	315.0 / 128.0) *	11 ✔
498	std::cos(phi));	11 ✔
499	// l = 9, m = 2
500	rn[i + 11] =	11 ✔
501	0.0158910431540932 * (w2m1) *	22 ✔
502	((109395.0 / 16.0) * std::pow(w, 7) - 135135.0 / 16.0 * std::pow(w, 5) +	11 ✔
503	(45045.0 / 16.0) * std::pow(w, 3) - 3465.0 / 16.0 * w) *	11 ✔
504	std::cos(2. * phi);	11 ✔
505	// l = 9, m = 3
506	rn[i + 12] = -(	11 ✔
507	0.00173385495536766 * std::pow(w2m1, 1.5) *	22 ✔
508	((765765.0 / 16.0) * std::pow(w, 6) - 675675.0 / 16.0 * std::pow(w, 4) +	11 ✔
509	(135135.0 / 16.0) * w * w - 3465.0 / 16.0) *	11 ✔
510	std::cos(3. * phi));	11 ✔
511	// l = 9, m = 4
512	rn[i + 13] = 0.000196320414650061 * w2m1 * w2m1 *	22 ✔
513	((2297295.0 / 8.0) * std::pow(w, 5) -	11 ✔
514	675675.0 / 4.0 * std::pow(w, 3) + (135135.0 / 8.0) * w) *	11 ✔
515	std::cos(4.0 * phi);	11 ✔
516	// l = 9, m = 5
517	rn[i + 14] = -(2.34647776186144e-5 * std::pow(w2m1, 2.5) *	22 ✔
518	((11486475.0 / 8.0) * std::pow(w, 4) -	11 ✔
519	2027025.0 / 4.0 * w * w + 135135.0 / 8.0) *	11 ✔
520	std::cos(5.0 * phi));	11 ✔
521	// l = 9, m = 6
522	rn[i + 15] = 3.02928976464514e-6 * std::pow(w2m1, 3) *	11 ✔
523	((11486475.0 / 2.) * std::pow(w, 3) - 2027025.0 / 2. * w) *	11 ✔
524	std::cos(6.0 * phi);	11 ✔
525	// l = 9, m = 7
526	rn[i + 16] =	11 ✔
527	-(4.37240315267812e-7 * std::pow(w2m1, 3.5) *	11 ✔
528	((34459425.0 / 2.) * w * w - 2027025.0 / 2.) * std::cos(7.0 * phi));	11 ✔
529	// l = 9, m = 8
530	rn[i + 17] =	11 ✔
531	2.58397773170915 * w * std::pow(w2m1, 4) * std::cos(8.0 * phi);	11 ✔
532	// l = 9, m = 9
533	rn[i + 18] =	11 ✔
534	-(0.609049392175524 * std::pow(w2m1, 4.5) * std::cos(9.0 * phi));	11 ✔
535	break;	11 ✔
536	case 10:	11 ✔
537	// l = 10, m = -10
538	rn[i] = 0.593627917136573 * std::pow(w2m1, 5) * std::sin(10.0 * phi);	11 ✔
539	// l = 10, m = -9
540	rn[i + 1] =	11 ✔
541	-(2.65478475211798 * w * std::pow(w2m1, 4.5) * std::sin(9.0 * phi));	11 ✔
542	// l = 10, m = -8
543	rn[i + 2] = 2.49953651452314e-8 * std::pow(w2m1, 4) *	11 ✔
544	((654729075.0 / 2.) * w * w - 34459425.0 / 2.) *	11 ✔
545	std::sin(8.0 * phi);	11 ✔
546	// l = 10, m = -7
547	rn[i + 3] =	11 ✔
548	-(1.83677671621093e-7 * std::pow(w2m1, 3.5) *	22 ✔
549	((218243025.0 / 2.) * std::pow(w, 3) - 34459425.0 / 2. * w) *	11 ✔
550	std::sin(7.0 * phi));	11 ✔
551	// l = 10, m = -6
552	rn[i + 4] = 1.51464488232257e-6 * std::pow(w2m1, 3) *	11 ✔
553	((218243025.0 / 8.0) * std::pow(w, 4) -	11 ✔
554	34459425.0 / 4.0 * w * w + 2027025.0 / 8.0) *	11 ✔
555	std::sin(6.0 * phi);	11 ✔
556	// l = 10, m = -5
557	rn[i + 5] =	11 ✔
558	-(1.35473956745817e-5 * std::pow(w2m1, 2.5) *	22 ✔
559	((43648605.0 / 8.0) * std::pow(w, 5) -	11 ✔
560	11486475.0 / 4.0 * std::pow(w, 3) + (2027025.0 / 8.0) * w) *	11 ✔
561	std::sin(5.0 * phi));	11 ✔
562	// l = 10, m = -4
563	rn[i + 6] = 0.000128521880085575 * w2m1 * w2m1 *	22 ✔
564	((14549535.0 / 16.0) * std::pow(w, 6) -	11 ✔
565	11486475.0 / 16.0 * std::pow(w, 4) +	11 ✔
566	(2027025.0 / 16.0) * w * w - 45045.0 / 16.0) *	11 ✔
567	std::sin(4.0 * phi);	11 ✔
568	// l = 10, m = -3
569	rn[i + 7] = -(0.00127230170115096 * std::pow(w2m1, 1.5) *	22 ✔
570	((2078505.0 / 16.0) * std::pow(w, 7) -	11 ✔
571	2297295.0 / 16.0 * std::pow(w, 5) +	11 ✔
572	(675675.0 / 16.0) * std::pow(w, 3) - 45045.0 / 16.0 * w) *	11 ✔
573	std::sin(3. * phi));	11 ✔
574	// l = 10, m = -2
575	rn[i + 8] = 0.012974982402692 * (w2m1) *	22 ✔
576	((2078505.0 / 128.0) * std::pow(w, 8) -	11 ✔
577	765765.0 / 32.0 * std::pow(w, 6) +	11 ✔
578	(675675.0 / 64.0) * std::pow(w, 4) -	11 ✔
579	45045.0 / 32.0 * w * w + 3465.0 / 128.0) *	11 ✔
580	std::sin(2. * phi);	11 ✔
581	// l = 10, m = -1
582	rn[i + 9] = -(0.134839972492648 * std::sqrt(w2m1) *	22 ✔
583	((230945.0 / 128.0) * std::pow(w, 9) -	11 ✔
584	109395.0 / 32.0 * std::pow(w, 7) +	11 ✔
585	(135135.0 / 64.0) * std::pow(w, 5) -	11 ✔
586	15015.0 / 32.0 * std::pow(w, 3) + (3465.0 / 128.0) * w) *	11 ✔
587	std::sin(phi));	11 ✔
588	// l = 10, m = 0
589	rn[i + 10] = 180.42578125 * std::pow(w, 10) -	11 ✔
590	427.32421875 * std::pow(w, 8) +	11 ✔
591	351.9140625 * std::pow(w, 6) - 117.3046875 * std::pow(w, 4) +	11 ✔
592	13.53515625 * w * w - 0.24609375;	11 ✔
593	// l = 10, m = 1
594	rn[i + 11] = -(0.134839972492648 * std::sqrt(w2m1) *	22 ✔
595	((230945.0 / 128.0) * std::pow(w, 9) -	11 ✔
596	109395.0 / 32.0 * std::pow(w, 7) +	11 ✔
597	(135135.0 / 64.0) * std::pow(w, 5) -	11 ✔
598	15015.0 / 32.0 * std::pow(w, 3) + (3465.0 / 128.0) * w) *	11 ✔
599	std::cos(phi));	11 ✔
600	// l = 10, m = 2
601	rn[i + 12] = 0.012974982402692 * (w2m1) *	22 ✔
602	((2078505.0 / 128.0) * std::pow(w, 8) -	11 ✔
603	765765.0 / 32.0 * std::pow(w, 6) +	11 ✔
604	(675675.0 / 64.0) * std::pow(w, 4) -	11 ✔
605	45045.0 / 32.0 * w * w + 3465.0 / 128.0) *	11 ✔
606	std::cos(2. * phi);	11 ✔
607	// l = 10, m = 3
608	rn[i + 13] =	11 ✔
609	-(0.00127230170115096 * std::pow(w2m1, 1.5) *	22 ✔
610	((2078505.0 / 16.0) * std::pow(w, 7) -	11 ✔
611	2297295.0 / 16.0 * std::pow(w, 5) +	11 ✔
612	(675675.0 / 16.0) * std::pow(w, 3) - 45045.0 / 16.0 * w) *	11 ✔
613	std::cos(3. * phi));	11 ✔
614	// l = 10, m = 4
615	rn[i + 14] = 0.000128521880085575 * w2m1 * w2m1 *	22 ✔
616	((14549535.0 / 16.0) * std::pow(w, 6) -	11 ✔
617	11486475.0 / 16.0 * std::pow(w, 4) +	11 ✔
618	(2027025.0 / 16.0) * w * w - 45045.0 / 16.0) *	11 ✔
619	std::cos(4.0 * phi);	11 ✔
620	// l = 10, m = 5
621	rn[i + 15] =	11 ✔
622	-(1.35473956745817e-5 * std::pow(w2m1, 2.5) *	22 ✔
623	((43648605.0 / 8.0) * std::pow(w, 5) -	11 ✔
624	11486475.0 / 4.0 * std::pow(w, 3) + (2027025.0 / 8.0) * w) *	11 ✔
625	std::cos(5.0 * phi));	11 ✔
626	// l = 10, m = 6
627	rn[i + 16] = 1.51464488232257e-6 * std::pow(w2m1, 3) *	11 ✔
628	((218243025.0 / 8.0) * std::pow(w, 4) -	11 ✔
629	34459425.0 / 4.0 * w * w + 2027025.0 / 8.0) *	11 ✔
630	std::cos(6.0 * phi);	11 ✔
631	// l = 10, m = 7
632	rn[i + 17] =	11 ✔
633	-(1.83677671621093e-7 * std::pow(w2m1, 3.5) *	22 ✔
634	((218243025.0 / 2.) * std::pow(w, 3) - 34459425.0 / 2. * w) *	11 ✔
635	std::cos(7.0 * phi));	11 ✔
636	// l = 10, m = 8
637	rn[i + 18] = 2.49953651452314e-8 * std::pow(w2m1, 4) *	11 ✔
638	((654729075.0 / 2.) * w * w - 34459425.0 / 2.) *	11 ✔
639	std::cos(8.0 * phi);	11 ✔
640	// l = 10, m = 9
641	rn[i + 19] =	11 ✔
642	-(2.65478475211798 * w * std::pow(w2m1, 4.5) * std::cos(9.0 * phi));	11 ✔
643	// l = 10, m = 10
644	rn[i + 20] = 0.593627917136573 * std::pow(w2m1, 5) * std::cos(10.0 * phi);	11 ✔
645	}
646	}
647	}	4,568,344 ✔
648
649	void calc_zn(int n, double rho, double phi, double zn[])	4,670,416 ✔
650	{
651	// ===========================================================================
652	// Determine vector of sin(nphi) and cos(nphi). This takes advantage of the
653	// following recurrence relations so that only a single sin/cos have to be
654	// evaluated (https://mathworld.wolfram.com/Multiple-AngleFormulas.html)
655	//
656	// sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
657	// cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
658
659	double sin_phi = std::sin(phi);	4,670,416 ✔
660	double cos_phi = std::cos(phi);	4,670,416 ✔
661
662	vector<double> sin_phi_vec(n + 1); // Sin[n * phi]	4,670,416 ✔
663	vector<double> cos_phi_vec(n + 1); // Cos[n * phi]	4,670,416 ✔
664	sin_phi_vec[0] = 1.0;	4,670,416 ✔
665	cos_phi_vec[0] = 1.0;	4,670,416 ✔
666	sin_phi_vec[1] = 2.0 * cos_phi;	4,670,416 ✔
667	cos_phi_vec[1] = cos_phi;	4,670,416 ✔
668
669	for (int i = 2; i <= n; i++) {	23,348,819 ✔
670	sin_phi_vec[i] = 2. * cos_phi * sin_phi_vec[i - 1] - sin_phi_vec[i - 2];	18,678,403 ✔
671	cos_phi_vec[i] = 2. * cos_phi * cos_phi_vec[i - 1] - cos_phi_vec[i - 2];	18,678,403 ✔
672	}
673
674	for (int i = 0; i <= n; i++) {	32,689,651 ✔
675	sin_phi_vec[i] *= sin_phi;	28,019,235 ✔
676	}
677
678	// ===========================================================================
679	// Calculate R_pq(rho)
680	// Matrix forms of the coefficients which are easier to work with
681	vector<vector<double>> zn_mat(n + 1, vector<double>(n + 1));	9,340,832 ✔
682
683	// Fill the main diagonal first (Eq 3.9 in Chong)
684	for (int p = 0; p <= n; p++) {	32,689,651 ✔
685	zn_mat[p][p] = std::pow(rho, p);	28,019,235 ✔
686	}
687
688	// Fill the 2nd diagonal (Eq 3.10 in Chong)
689	for (int q = 0; q <= n - 2; q++) {	23,348,819 ✔
690	zn_mat[q][q + 2] = (q + 2) * zn_mat[q + 2][q + 2] - (q + 1) * zn_mat[q][q];	18,678,403 ✔
691	}
692
693	// Fill in the rest of the values using the original results (Eq. 3.8 in
694	// Chong)
695	for (int p = 4; p <= n; p++) {	14,007,987 ✔
696	double k2 = 2 * p * (p - 1) * (p - 2);	9,337,571 ✔
697	for (int q = p - 4; q >= 0; q -= 2) {	18,675,241 ✔
698	double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;	9,337,670 ✔
699	double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);	9,337,670 ✔
700	double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;	9,337,670 ✔
701	zn_mat[q][p] =	9,337,670 ✔
702	((k2 * rho * rho + k3) * zn_mat[q][p - 2] + k4 * zn_mat[q][p - 4]) / k1;	9,337,670 ✔
703	}
704	}
705
706	// Roll into a single vector for easier computation later
707	// The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
708	// (2, 2), .... in (n,m) indices
709	// Note that the cos and sin vectors are offset by one
710	// sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
711	// cos_phi_vec = [1.0, cos(x), cos(2x)... ]
712	int i = 0;	4,670,416 ✔
713	for (int p = 0; p <= n; p++) {	32,689,651 ✔
714	for (int q = -p; q <= p; q += 2) {	126,080,066 ✔
715	if (q < 0) {	98,060,831 ✔
716	zn[i] = zn_mat[std::abs(q)][p] * sin_phi_vec[std::abs(q) - 1];	42,025,523 ✔
717	} else if (q == 0) {	56,035,308 ✔
718	zn[i] = zn_mat[q][p];	14,009,785 ✔
719	} else {
720	zn[i] = zn_mat[q][p] * cos_phi_vec[q];	42,025,523 ✔
721	}
722	i++;	98,060,831 ✔
723	}
724	}
725	}	4,670,416 ✔
726
727	void calc_zn_rad(int n, double rho, double zn_rad[])	65 ✔
728	{
729	// Calculate R_p0(rho) as Zn_p0(rho)
730	// Set up the array of the coefficients
731
732	double q = 0;	65 ✔
733
734	// R_00 is always 1
735	zn_rad[0] = 1;	65 ✔
736
737	// Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong)
738	for (int p = 2; p <= n; p += 2) {	390 ✔
739	int index = int(p / 2);	325 ✔
740	if (p == 2) {	325 ✔
741	// Setting up R_22 to calculate R_20 (Eq 3.10 in Chong)
742	double R_22 = rho * rho;	65 ✔
743	zn_rad[index] = 2 * R_22 - zn_rad[0];	65 ✔
744	} else {
745	double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;	260 ✔
746	double k2 = 2 * p * (p - 1) * (p - 2);	260 ✔
747	double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);	260 ✔
748	double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;	260 ✔
749	zn_rad[index] =	260 ✔
750	((k2 * rho * rho + k3) * zn_rad[index - 1] + k4 * zn_rad[index - 2]) /	260 ✔
751	k1;
752	}
753	}
754	}	65 ✔
755
756	void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* seed)	33 ✔
757	{
758	Direction u = rotate_angle({uvw}, mu, phi, seed);	33 ✔
759	uvw[0] = u.x;	33 ✔
760	uvw[1] = u.y;	33 ✔
761	uvw[2] = u.z;	33 ✔
762	}	33 ✔
763
764	Direction rotate_angle(	2,147,483,647 ✔
765	Direction u, double mu, const double* phi, uint64_t* seed)
766	{
767	// Sample azimuthal angle in [0,2pi) if none provided
768	double phi_;
769	if (phi != nullptr) {	2,147,483,647 ✔
770	phi_ = (*phi);	22,500,941 ✔
771	} else {
772	phi_ = 2.0 * PI * prn(seed);	2,147,483,647 ✔
773	}
774
775	// Precompute factors to save flops
776	double sinphi = std::sin(phi_);	2,147,483,647 ✔
777	double cosphi = std::cos(phi_);	2,147,483,647 ✔
778	double a = std::sqrt(std::fmax(0., 1. - mu * mu));	2,147,483,647 ✔
779	double b = std::sqrt(std::fmax(0., 1. - u.z * u.z));	2,147,483,647 ✔
780
781	// Need to treat special case where sqrt(1 - w**2) is close to zero by
782	// expanding about the v component rather than the w component
783	if (b > 1e-10) {	2,147,483,647 ✔
784	return {mu * u.x + a * (u.x * u.z * cosphi - u.y * sinphi) / b,	2,147,483,647 ✔
785	mu * u.y + a * (u.y * u.z * cosphi + u.x * sinphi) / b,	2,147,483,647 ✔
786	mu * u.z - a * b * cosphi};	2,147,483,647 ✔
787	} else {
788	b = std::sqrt(1. - u.y * u.y);	390,907 ✔
789	return {mu * u.x + a * (-u.x * u.y * sinphi + u.z * cosphi) / b,	390,907 ✔
790	mu * u.y + a * b * sinphi,	390,907 ✔
791	mu * u.z - a * (u.y * u.z * sinphi + u.x * cosphi) / b};	390,907 ✔
792	}
793	}
794
795	void spline(int n, const double x[], const double y[], double z[])	125,970 ✔
796	{
797	vector<double> c_new(n - 1);	125,970 ✔
798
799	// Set natural boundary conditions
800	c_new[0] = 0.0;	125,970 ✔
801	z[0] = 0.0;	125,970 ✔
802	z[n - 1] = 0.0;	125,970 ✔
803
804	// Solve using tridiagonal matrix algorithm; first do forward sweep
805	for (int i = 1; i < n - 1; i++) {	12,628,090 ✔
806	double a = x[i] - x[i - 1];	12,502,120 ✔
807	double c = x[i + 1] - x[i];	12,502,120 ✔
808	double b = 2.0 * (a + c);	12,502,120 ✔
809	double d = 6.0 * ((y[i + 1] - y[i]) / c - (y[i] - y[i - 1]) / a);	12,502,120 ✔
810
811	c_new[i] = c / (b - a * c_new[i - 1]);	12,502,120 ✔
812	z[i] = (d - a * z[i - 1]) / (b - a * c_new[i - 1]);	12,502,120 ✔
813	}
814
815	// Back substitution
816	for (int i = n - 2; i >= 0; i--) {	12,754,060 ✔
817	z[i] = z[i] - c_new[i] * z[i + 1];	12,628,090 ✔
818	}
819	}	125,970 ✔
820
821	double spline_interpolate(	×
822	int n, const double x[], const double y[], const double z[], double xint)
823	{
824	// Find the lower bounding index in x of xint
825	int i = n - 1;	×
826	while (--i) {	×
827	if (xint >= x[i])	×
828	break;	×
829	}
830
831	double h = x[i + 1] - x[i];	×
832	double r = xint - x[i];	×
833
834	// Compute the coefficients
835	double b = (y[i + 1] - y[i]) / h - (h / 6.0) * (z[i + 1] + 2.0 * z[i]);	×
836	double c = z[i] / 2.0;	×
837	double d = (z[i + 1] - z[i]) / (h * 6.0);	×
838
839	return y[i] + b * r + c * r * r + d * r * r * r;	×
840	}
841
842	double spline_integrate(int n, const double x[], const double y[],	12,628,090 ✔
843	const double z[], double xa, double xb)
844	{
845	// Find the lower bounding index in x of the lower limit of integration.
846	int ia = n - 1;	12,628,090 ✔
847	while (--ia) {	844,174,998 ✔
848	if (xa >= x[ia])	844,049,028 ✔
849	break;	12,502,120 ✔
850	}
851
852	// Find the lower bounding index in x of the upper limit of integration.
853	int ib = n - 1;	12,628,090 ✔
854	while (--ib) {	831,672,878 !
855	if (xb >= x[ib])	831,672,878 ✔
856	break;	12,628,090 ✔
857	}
858
859	// Evaluate the integral
860	double s = 0.0;	12,628,090 ✔
861	for (int i = ia; i <= ib; i++) {	37,758,300 ✔
862	double h = x[i + 1] - x[i];	25,130,210 ✔
863
864	// Compute the coefficients
865	double b = (y[i + 1] - y[i]) / h - (h / 6.0) * (z[i + 1] + 2.0 * z[i]);	25,130,210 ✔
866	double c = z[i] / 2.0;	25,130,210 ✔
867	double d = (z[i + 1] - z[i]) / (h * 6.0);	25,130,210 ✔
868
869	// Subtract the integral from x[ia] to xa
870	if (i == ia) {	25,130,210 ✔
871	double r = xa - x[ia];	12,628,090 ✔
872	s = s - (y[i] * r + b / 2.0 * r * r + c / 3.0 * r * r * r +	12,628,090 ✔
873	d / 4.0 * r * r * r * r);	12,628,090 ✔
874	}
875
876	// Integrate from x[ib] to xb in final interval
877	if (i == ib) {	25,130,210 ✔
878	h = xb - x[ib];	12,628,090 ✔
879	}
880
881	// Accumulate the integral
882	s = s + y[i] * h + b / 2.0 * h * h + c / 3.0 * h * h * h +	25,130,210 ✔
883	d / 4.0 * h * h * h * h;	25,130,210 ✔
884	}
885
886	return s;	12,628,090 ✔
887	}
888
889	std::complex<double> faddeeva(std::complex<double> z)	3,505,777 ✔
890	{
891	// Technically, the value we want is given by the equation:
892	// w(z) = I/pi * Integrate[Exp[-t^2]/(z-t), {t, -Infinity, Infinity}]
893	// as shown in Equation 63 from Hwang, R. N. "A rigorous pole
894	// representation of multilevel cross sections and its practical
895	// applications." Nucl. Sci. Eng. 96.3 (1987): 192-209.
896	//
897	// The MIT Faddeeva function evaluates w(z) = exp(-z^2)erfc(-iz). These
898	// two forms of the Faddeeva function are related by a transformation.
899	//
900	// If we call the integral form w_int, and the function form w_fun:
901	// For imag(z) > 0, w_int(z) = w_fun(z)
902	// For imag(z) < 0, w_int(z) = -conjg(w_fun(conjg(z)))
903
904	// Note that Faddeeva::w will interpret zero as machine epsilon
905	return z.imag() > 0.0 ? Faddeeva::w(z)	3,505,777 !
906	: -std::conj(Faddeeva::w(std::conj(z)));	7,011,554 !
907	}
908
909	std::complex<double> w_derivative(std::complex<double> z, int order)	4,232,547 ✔
910	{
911	using namespace std::complex_literals;
912	switch (order) {	4,232,547 ✔
913	case 0:	1,410,849 ✔
914	return faddeeva(z);	1,410,849 ✔
915	case 1:	1,410,849 ✔
916	return -2.0 * z * faddeeva(z) + 2.0i / SQRT_PI;	2,821,698 ✔
917	default:	1,410,849 ✔
918	return -2.0 * z * w_derivative(z, order - 1) -	1,410,849 ✔
919	2.0 * (order - 1) * w_derivative(z, order - 2);	4,232,547 ✔
920	}
921	}
922
NEW 923	double exprel(double x)	×
924	{
NEW 925	if (std::abs(x) < 1e-16)	×
NEW 926	return 1.0;	×
927	else {
NEW 928	return std::expm1(x) / x;	×
929	}
930	}
931
NEW 932	double log1prel(double x)	×
933	{
NEW 934	if (std::abs(x) < 1e-16)	×
NEW 935	return 1.0;	×
936	else {
NEW 937	return std::log1p(x) / x;	×
938	}
939	}
940
941	double lambert_w0(double x)	143 ✔
942	{
943	return LambertW::lambert_w0(x);	143 ✔
944	}
945
946	double lambert_wm1(double x)	33 ✔
947	{
948	return LambertW::lambert_wm1(x);	33 ✔
949	}
950
951	} // namespace openmc

openmc-dev / openmc / 19072640351

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous