34177814562

Committed 08 Sep 2026 01:47AM UTC coverage: 81.364% (+0.004%) from 81.36%

Build # 34177814562

Build Type

Pull #4113

github

Committed by

web-flow

Commit Message

Merge c53f9baff into d7d3284a1

Pull Request Pull Request #4113: Extract the combined k-effective estimator and fix its two-estimate branch

Coverage Stats

18701 of 27185 branches covered (68.79%)

Branch coverage included in aggregate %.

64 of 66 new or added lines in 2 files covered. (96.97%)

1 existing line in 1 file now uncovered.

60597 of 70276 relevant lines covered (86.23%)

49874814.6 hits per line

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

94.77

/src/math_functions.cpp

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

openmc-dev / openmc / 34177814562

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