Repository navigation
Expand file tree
/
Copy pathmath_functions.h
More file actions
263 lines (225 loc) · 10.8 KB
/
Copy pathmath_functions.h
File metadata and controls
263 lines (225 loc) · 10.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
//! \file math_functions.h
//! A collection of elementary math functions.
#ifndef OPENMC_MATH_FUNCTIONS_H
#define OPENMC_MATH_FUNCTIONS_H
#include <cmath>
#include <complex>
#include <cstdlib>
#include "openmc/position.h"
#include "openmc/search.h"
namespace openmc {
//==============================================================================
//! Calculate the percentile of the standard normal distribution with a
//! specified probability level.
//!
//! \param p The probability level
//! \return The requested percentile
//==============================================================================
extern "C" double normal_percentile(double p);
//==============================================================================
//! Calculate the percentile of the Student's t distribution with a specified
//! probability level and number of degrees of freedom.
//!
//! \param p The probability level
//! \param df The degrees of freedom
//! \return The requested percentile
//==============================================================================
extern "C" double t_percentile(double p, int df);
//==============================================================================
//! Calculate the n-th order Legendre polynomials at the value of x.
//!
//! \param n The maximum order requested
//! \param x The value to evaluate at; x is expected to be within [-1,1]
//! \param pnx The requested Legendre polynomials of order 0 to n (inclusive)
//! evaluated at x.
//==============================================================================
extern "C" void calc_pn_c(int n, double x, double pnx[]);
//==============================================================================
//! Find the value of f(x) given a set of Legendre coefficients and the value
//! of x.
//!
//! \param n The maximum order of the expansion
//! \param data The polynomial expansion coefficient data; without the (2l+1)/2
//! factor.
//! \param x The value to evaluate at; x is expected to be within [-1,1]
//! \return The requested Legendre polynomials of order 0 to n (inclusive)
//! evaluated at x
//==============================================================================
extern "C" double evaluate_legendre(int n, const double data[], double x);
//==============================================================================
//! Calculate the n-th order real spherical harmonics for a given angle (in
//! terms of (u,v,w)) for all 0<=n and -m<=n<=n.
//!
//! \param n The maximum order requested
//! \param uvw[3] The direction the harmonics are requested at
//! \param rn The requested harmonics of order 0 to n (inclusive)
//! evaluated at uvw.
//==============================================================================
extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
void calc_rn(int n, Direction u, double rn[]);
//==============================================================================
//! Calculate the n-th order modified Zernike polynomial moment for a given
//! angle (rho, theta) location on the unit disk.
//!
//! This procedure uses the modified Kintner's method for calculating Zernike
//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
//! R. (2003). A comparative analysis of algorithms for fast computation of
//! Zernike moments. Pattern Recognition, 36(3), 731-742.
//! The normalization of the polynomials is such that the integral of Z_pq^2
//! over the unit disk is exactly pi.
//!
//! \param n The maximum order requested
//! \param rho The radial parameter to specify location on the unit disk
//! \param phi The angle parameter to specify location on the unit disk
//! \param zn The requested moments of order 0 to n (inclusive)
//! evaluated at rho and phi.
//==============================================================================
extern "C" void calc_zn(int n, double rho, double phi, double zn[]);
//==============================================================================
//! Calculate only the even radial components of n-th order modified Zernike
//! polynomial moment with azimuthal dependency m = 0 for a given angle
//! (rho, theta) location on the unit disk.
//!
//! Since m = 0, n could only be even orders. Z_q0 = R_q0
//!
//! This procedure uses the modified Kintner's method for calculating Zernike
//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
//! R. (2003). A comparative analysis of algorithms for fast computation of
//! Zernike moments. Pattern Recognition, 36(3), 731-742.
//! The normalization of the polynomials is such that the integral of Z_pq^2
//! over the unit disk is exactly pi.
//!
//! \param n The maximum order requested
//! \param rho The radial parameter to specify location on the unit disk
//! \param phi The angle parameter to specify location on the unit disk
//! \param zn_rad The requested moments of order 0 to n (inclusive)
//! evaluated at rho and phi when m = 0.
//==============================================================================
extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
//==============================================================================
//! Rotate the direction cosines through a polar angle whose cosine is mu and
//! through an azimuthal angle sampled uniformly.
//!
//! This is done with direct sampling rather than rejection sampling as is done
//! in MCNP and Serpent.
//!
//! \param uvw[3] The initial, and final, direction vector
//! \param mu The cosine of angle in lab or CM
//! \param phi The azimuthal angle; will randomly chosen angle if a nullptr
//! is passed
//! \param seed A pointer to the pseudorandom seed
//==============================================================================
extern "C" void rotate_angle_c(
double uvw[3], double mu, const double* phi, uint64_t* seed);
Direction rotate_angle(
Direction u, double mu, const double* phi, uint64_t* seed);
//==============================================================================
//! Constructs a natural cubic spline.
//!
//! Given a tabulated function y_i = f(x_i), this computes the second
//! derivative of the interpolating function at each x_i, which can then be
//! used in any subsequent calls to spline_interpolate or spline_integrate for
//! the same set of x and y values.
//!
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! \param y Values of the dependent variable.
//! \param[out] z The second derivative of the interpolating function at each
//! value of x.
//==============================================================================
void spline(int n, const double x[], const double y[], double z[]);
//==============================================================================
//! Determine the cubic spline interpolated y-value for a given x-value.
//!
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! \param y Values of the dependent variable.
//! \param z The second derivative of the interpolating function at each
//! value of x.
//! \param xint Point at which to evaluate the cubic spline polynomial
//! \return Interpolated value
//==============================================================================
double spline_interpolate(
int n, const double x[], const double y[], const double z[], double xint);
//==============================================================================
//! Evaluate the definite integral of the interpolating cubic spline between
//! the given endpoints.
//!
//! \param n Number of points
//! \param x Values of the independent variable, which must be strictly
//! increasing.
//! \param y Values of the dependent variable.
//! \param z The second derivative of the interpolating function at each
//! value of x.
//! \param xa Lower limit of integration
//! \param xb Upper limit of integration
//! \return Integral
//==============================================================================
double spline_integrate(int n, const double x[], const double y[],
const double z[], double xa, double xb);
//! Evaluate the Faddeeva function
//!
//! \param z Complex argument
//! \return Faddeeva function evaluated at z
std::complex<double> faddeeva(std::complex<double> z);
//! Evaluate derivative of the Faddeeva function
//!
//! \param z Complex argument
//! \param order Order of the derivative
//! \return Derivative of Faddeeva function evaluated at z
std::complex<double> w_derivative(std::complex<double> z, int order);
//! Evaluate relative exponential function
//!
//! \param x Real argument
//! \return (exp(x)-1)/x without loss of precision near 0
double exprel(double x);
//! Evaluate relative logarithm function
//!
//! \param x Real argument
//! \return log(1+x)/x without loss of precision near 0
double log1prel(double x);
//! Evaluate the cylindrical Bessel function of the first kind J_n(x)
//!
//! Uses std::cyl_bessel_j where available (e.g., libstdc++). On standard
//! library implementations lacking the C++17 special math functions (e.g.,
//! libc++ on Apple Clang/LLVM), falls back to the ascending power series,
//! which converges to machine precision for the small arguments (|x| <= 2)
//! used in OpenMC. Unlike std::cyl_bessel_j, negative arguments are handled
//! via the parity relation J_n(-x) = (-1)^n J_n(x).
//!
//! \param n Non-negative integer order of the Bessel function
//! \param x Real argument
//! \return J_n(x)
double cyl_bessel_j(int n, double x);
//! Helper function to get index and interpolation function on an incident
//! energy grid
//!
//! \param energies energy grid
//! \param E incident energy
//! \param i grid index
//! \param f interpolation factor
void get_energy_index(
const vector<double>& energies, double E, int& i, double& f);
//==============================================================================
//! Calculate the cumulative distribution function of the standard normal
//! distribution at a given value.
//!
//! \param z The value at which to evaluate the CDF
//! \return Phi(z) = P(X <= z) for X ~ N(0,1)
//==============================================================================
double standard_normal_cdf(double z);
//==============================================================================
//! Return true if two floating-point values are approximately equal within a
//! combined relative and absolute tolerance.
//!
//! \param a first floating point value
//! \param b second floating point value
//! \param rel_tol relative tolerance
//! \param abs_tol absolute tolerance
//! \return true if a and b are approximately equal, false otherwise
//==============================================================================
bool isclose(double a, double b, double rel_tol, double abs_tol);
} // namespace openmc
#endif // OPENMC_MATH_FUNCTIONS_H