21960815265

Committed 12 Feb 2026 07:17PM UTC coverage: 81.875% (+0.02%) from 81.851%

Build # 21960815265

Build Type

Pull #3751

github

Committed by

web-flow

Commit Message

Merge 942f5e16f into 8198e6021

Pull Request Pull Request #3751: Resolve conflict with weight windows and global russian roulette

Coverage Stats

17421 of 24368 branches covered (71.49%)

Branch coverage included in aggregate %.

85 of 90 new or added lines in 5 files covered. (94.44%)

121 existing lines in 4 files now uncovered.

56462 of 65871 relevant lines covered (85.72%)

44576043.34 hits per line

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

94.85

/src/math_functions.cpp

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

openmc-dev / openmc / 21960815265

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