29833350562

Committed 21 Jul 2026 01:11PM UTC coverage: 81.257% (-0.04%) from 81.301%

Build # 29833350562

Build Type

Pull #4012

github

Committed by

web-flow

Commit Message

Merge a7791fc9c into 61c8a59cf

Pull Request Pull Request #4012: Geometry debug function

Coverage Stats

18354 of 26615 branches covered (68.96%)

Branch coverage included in aggregate %.

5 of 103 new or added lines in 1 file covered. (4.85%)

227 existing lines in 11 files now uncovered.

59875 of 69659 relevant lines covered (85.95%)

47832103.0 hits per line

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

95.76

/src/math_functions.cpp

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

openmc-dev / openmc / 29833350562

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