39072696583

Committed 17 Oct 2026 09:51PM UTC coverage: 81.441% (-0.005%) from 81.446%

Build # 39072696583

Build Type

Pull #4202

github

Committed by

web-flow

Commit Message

Merge f30434957 into 98bd86f5c

Pull Request Pull Request #4202: Bit-reproducible transport across platforms with an opt-in correctly rounded math backend

Coverage Stats

20118 of 29242 branches covered (68.8%)

Branch coverage included in aggregate %.

571 of 595 new or added lines in 42 files covered. (95.97%)

62699 of 72447 relevant lines covered (86.54%)

42721094.23 hits per line

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

95.63

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

openmc-dev / openmc / 39072696583

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