19476608580

Committed 18 Nov 2025 06:15PM UTC coverage: 81.528% (-0.4%) from 81.961%

Build # 19476608580

Build Type

Pull #3339

github

Committed by

web-flow

Commit Message

Merge 25c21d84d into 7815d3a68

Pull Request Pull Request #3339: Enable Scikit-Build-Core Support

Run Details

16890 of 23555 branches covered (71.7%)

Branch coverage included in aggregate %.

34 of 34 new or added lines in 3 files covered. (100.0%)

376 existing lines in 7 files now uncovered.

54523 of 64038 relevant lines covered (85.14%)

42168701.72 hits per line

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

47.56

/src/math_functions.cpp

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

openmc-dev / openmc / 19476608580

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