15567165838

Committed 10 Jun 2025 06:18PM UTC coverage: 85.1% (-0.06%) from 85.158%

Build # 15567165838

Build Type

Pull #3413

github

Committed by

web-flow

Commit Message

Merge ec68c734c into f796fa04e

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

Run Details

69 of 121 new or added lines in 3 files covered. (57.02%)

34 existing lines in 2 files now uncovered.

52396 of 61570 relevant lines covered (85.1%)

36584924.15 hits per line

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

97.23

/src/math_functions.cpp

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

openmc-dev / openmc / 15567165838

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