Utility
-
double SBody::absolute_accuracy
Global absolute accuracy.
-
double SBody::relative_accuracy
Global relative accuracy.
-
constexpr double SBody::M_2PI = 6.28318530717958647692528676655900576
\(2\pi\).
-
constexpr double SBody::M_PI2 = 9.86960440108935861883449099987615111
\(\pi^2\).
-
constexpr double SBody::M_SQRT27 = 5.19615242270663188058233902451761710
\(\sqrt{27}\)
-
constexpr double SBody::EPSILON = 1e-10
Epsilon, a small value. \(\varepsilon\).
-
constexpr double SBody::SIN_EPSILON = 1e-10
\(\sin\varepsilon\).
-
constexpr double SBody::COS_EPSILON = 0.999999999999999999995
\(\cos\varepsilon\).
-
constexpr double SBody::EPSILON_CIRCLE_AREA = M_PI * GSL_DBL_EPSILON
Area of a circle with radius of GSL_SQRT_DBL_EPSILON. \(\pi\varepsilon^2\).
-
constexpr int SBody::SAMPLE_NUMBER = 100
Sample number.
-
constexpr double SBody::ANGLE_INTERVAL = M_2PI / SAMPLE_NUMBER
The angle corresponding to the sample number.
-
constexpr double SBody::EPSILON_POLYGON_AREA = 0.06279051952931337 / ANGLE_INTERVAL * EPSILON_CIRCLE_AREA
Area of the regular polygon with SAMPLE_NUMBER edges.
Integrator
-
class Integrator
A wrapper of the gsl_odeiv2_evolve.
Public Functions
-
Integrator(int (*function)(double, const double*, double*, void*), int (*jacobian)(double, const double*, double*, double*, void*), void *params = nullptr, const gsl_odeiv2_step_type *type = gsl_odeiv2_step_rk8pd)
Construct a new Integrator object.
- Parameters:
function – The function calculates
\[\frac{\mathrm dy_i(t)}{\mathrm dt}=f_i\left[t,y_1(t),\dots,y_n(t)\right].\]jacobian – The jacobian of the function
\[J_{ij}=\frac{\partial f_i\left[t,y(t)\right]}{\partial y_j}.\]params – The parameters passed to the function, like the PN parameter, or the spin of the black hole.
type – Type of the algorithms. Explicit embedded Runge-Kutta Prince-Dormand (8, 9) method is set by default.
-
~Integrator()
Destructor.
-
int Apply(double *t, double t1, double *h, double *y)
gsl_odeiv2_evolve_apply.- Parameters:
t – time.
t1 – maximum time not to be exceeded by the time step.
h – step size.
y – values of the integrated system.
- Returns:
status
-
int ApplyStep(double *t, double t1, double *h, double *y)
gsl_odeiv2_evolve_apply.- Parameters:
t – time.
t1 – maximum time not to be exceeded by the time step.
h – step size.
y – values of the integrated system.
- Returns:
status
-
int ApplyFixedStep(double *t, const double h, double *y)
gsl_odeiv2_evolve_apply_fixed_step.- Parameters:
t – time.
h – step size.
y – values of the integrated system.
- Returns:
status
-
int Reset()
Resets the evolution function and the stepping function.
- Returns:
status
-
int CheckCoordinate(double *y)
- Parameters:
y –
- Returns:
status
-
Integrator(int (*function)(double, const double*, double*, void*), int (*jacobian)(double, const double*, double*, double*, void*), void *params = nullptr, const gsl_odeiv2_step_type *type = gsl_odeiv2_step_rk8pd)
Solver
-
class Solver
Subclassed by SBody::DerivativeSolver, SBody::FunctionSolver
FunctionSolver
-
class FunctionSolver : public SBody::Solver
DerivativeSolver
-
class DerivativeSolver : public SBody::Solver
MultiSolver
-
class MultiSolver
Subclassed by SBody::MultiDerivativeSolver, SBody::MultiFunctionSolver
MultiFunctionSolver
-
class MultiFunctionSolver : public SBody::MultiSolver
MultiDerivativeSolver
-
class MultiDerivativeSolver : public SBody::MultiSolver
GslBlock
-
class GslBlock
A wrapper of the
gsl_vector,gsl_matrix, andgsl_permutation.Public Functions
-
~GslBlock()
Destructor.
-
gsl_vector *VectorAlloc(size_t n)
gsl_vector_alloc, creates a vector of lengthn- Parameters:
n – length of the vector
- Returns:
pointer
-
gsl_vector *VectorCalloc(size_t n)
gsl_vector_calloc, creates a vector of lengthnand initializes all the elements of the vector to zero.- Parameters:
n – length of the vector
- Returns:
pointer
-
gsl_vector *VectorAllocFromBlock(gsl_block *block, const size_t offset, const size_t n, const size_t stride = 1)
gsl_vector_alloc_from_block, creates a vector as a slice of an existingblock.- Parameters:
block – target block
offset – offset to the data block
n – length of the vector
stride – step-size of the vector
- Returns:
pointer
-
gsl_vector *VectorAllocRowFromMatrix(gsl_matrix *matrix, const size_t i)
gsl_vector_alloc_row_from_matrix, allocates a newgsl_vectorwhich points to thei-th row of thematrix.- Parameters:
matrix – target matrix
i – row index
- Returns:
pointer
-
gsl_vector *VectorAllocColFromMatrix(gsl_matrix *matrix, const size_t j)
gsl_vector_alloc_col_from_matrix, allocates a newgsl_vectorwhich points to thej-th column of thematrix.- Parameters:
matrix – target matrix
j – column index
- Returns:
pointer
-
gsl_matrix *MatrixAlloc(size_t n1, size_t n2)
gsl_matrix_alloc, creates a matrix of sizen1rows byn2columns.- Parameters:
n1 – rows of the matrix
n2 – columns of the matrix
- Returns:
pointer
-
gsl_matrix *MatrixCalloc(size_t n1, size_t n2)
gsl_matrix_alloc, creates a matrix of sizen1rows byn2columns and initializes all the elements of the matrix to zero- Parameters:
n1 – rows of the matrix
n2 – columns of the matrix
- Returns:
pointer
-
gsl_permutation *PermutationAlloc(size_t n)
gsl_permutation_alloc, allocates memory for a new permutation of sizen- Parameters:
n – size of the permutation
- Returns:
pointer
-
gsl_permutation *PermutationCalloc(size_t n)
gsl_permutation_alloc, allocates memory for a new permutation of sizenand initializes it to the identity- Parameters:
n – size of the permutation
- Returns:
pointer
-
~GslBlock()
Math
-
double SBody::Dot(const double x[], const double y[], size_t dimension = 3)
Dot product of vector
xandy. \(\vec{x}\cdot\vec{y}\).- Parameters:
x – vector
y – vector
dimension – dimension of the vector
- Returns:
result
-
double SBody::Dot(const double x[], size_t dimension = 3)
Dot product of vector
x. \(\vec{x}\cdot\vec{x}\).- Parameters:
x – vector
dimension – dimension of the vector
- Returns:
result
-
double SBody::Norm(const double x[], size_t dimension = 3)
Euclidean norm of
x. \(\|\vec{x}\|_2\).- Parameters:
x – vector
dimension – dimension of the vector
- Returns:
result
-
int SBody::Cross(const double x[], const double y[], double z[])
Cross product of vector
xandy, stored inz. \(\vec{z}=\vec{x}\times\vec{y}\).- Parameters:
x – 3 dimensional vector
y – 3 dimensional vector
z – 3 dimensional vector
- Returns:
status
-
double SBody::DotCross(const double x[], const double y[], const double z[])
Dot product of vector
xand the cross product of vectoryand vectorz. \(\vec{x}\cdot(\vec{y}\times\vec{z})\).- Parameters:
x – 3 dimensional vector
y – 3 dimensional vector
z – 3 dimensional vector
- Returns:
result
-
double SBody::TriangleArea(double a, double b, double c)
-
double SBody::TriangleArea(const double x[], const double y[], const double z[])
-
int SBody::RotateAroundAxis(double x[], Axis axis, double angle)
Rotate vector
xaround theaxisbyangle.- Parameters:
x – 3 dimensional vector
axis – the rotation axis.
angle – in rad
- Returns:
status
-
int SBody::CartesianToSpherical(const double cartesian[], double spherical[], size_t dimension = 8)
- Parameters:
cartesian – 4 or 8 dimensional vector
spherical – 4 or 8 dimensional vector
dimension – 4 or 8.
- Returns:
status
-
int SBody::CartesianToSpherical(double x[], size_t dimension = 8)
- Parameters:
x – 4 or 8 dimensional vector
dimension – 4 or 8
- Returns:
status
-
int SBody::SphericalToCartesian(const double spherical[], double cartesian[], size_t dimension = 8)
- Parameters:
spherical – 4 or 8 dimensional vector
cartesian – 4 or 8 dimensional vector
dimension – 4 or 8.
- Returns:
status
-
int SBody::SphericalToCartesian(double x[], size_t dimension = 8)
- Parameters:
x – 4 or 8 dimensional vector
dimension – 4 or 8.
- Returns:
status
-
int SBody::OppositeSign(double x, double y)
Return
1ifx,yhave opposite signs, else0.- Parameters:
x – number
y – number
- Returns:
result
-
int SBody::MapTheta(const double theta_0, double y[])
The modified spherical coordiante maps the polar angle from \(\theta\in[0,\pi]\) to.
\[\begin{split} \theta'\equiv\begin{cases}\theta & (\theta\leq\pi/2)\\ \theta-\pi & (\theta>\pi/2)\end{cases}. \end{split}\]- Parameters:
theta_0 – theta in the last integration step
y – 8 dimensional vector
- Returns:
status
-
double SBody::ModBy2Pi(double phi)
Return
phiin \([0, 2\pi)\).- Parameters:
phi – \(\phi\)
- Returns:
result
-
double SBody::PhiDifference(double phi)
Similar to
ModBy2Pi, but returnphiin \([-\pi, \pi)\).- Parameters:
phi – \(\phi\)
- Returns:
result
-
double SBody::LinearInterpolation(double x, double x0, double x1, double y0, double y1)
Linear interpolation of points (
x0,y0) and (x1,y1) atx.\[y=\frac{y_0(x_1-x)+y_1(x-x_0)}{x_1-x_0}\]- Parameters:
x – position to evaluate \(y\)
x0 – \(x_0\)
x1 – \(x_1\)
y0 – \(y_0\)
y1 – \(y_1\)
- Returns:
result
-
int SBody::LinearInterpolation(double x, double x0, double x1, const double y0[], const double y1[], double y[], size_t size)
Linear interpolation of vectors (
x0,y0) and (x1,y1) atx, stored iny.\[y=\frac{y_0(x_1-x)+y_1(x-x_0)}{x_1-x_0}\]- Parameters:
x – position to evaluate \(y\)
x0 – \(x_0\)
x1 – \(x_1\)
y0 – \(y_0\)
y1 – \(y_1\)
y – \(y\)
size – size of the vectors
- Returns:
status
-
int SBody::InterpolateSphericalPositionToCartesian(double t, double t0, double t1, const double y0[], const double y1[], double y[])
Linear interpolation of two spherical positions
y0andy1att, stored iny.\[y=\frac{\text{SphericalToCartesian}(y_0)(t_1-t)+\text{SphericalToCartesian}(y_1)(t-t_0)}{t_1-t_0}\]- Parameters:
t – time to evaluate \(y\)
t0 – \(t_0\)
t1 – \(t_1\)
y0 – 8 dimensional vector, spherical, \(y_0\)
y1 – 8 dimensional vector, spherical, \(y_1\)
y – 8 dimensional vector, cartesian, \(y\)
- Returns:
status
-
double SBody::Flux(double luminosity, double magnification, double redshift)
Calculate the bolometric flux. The result is.
\[F_\text{obs}=\frac{L_\text{src}}{(1+z)^2}\frac{\Omega_\text{src}}{\Omega_\text{obs}}.\]The additional \((1+z)^{-1}\) comes from the expansion of the frequency that \((1+z)^{-1}=\frac{\nu_\text{obs}}{\nu_\text{src}}\).- Parameters:
luminosity – intrinsic luminosity of the source, \(L_\text{src}\)
magnification – magnification of the source, including the relativistic redshift and beaming, and gravitational lensing. \(\frac{\Omega_\text{src}}{(1+z)\Omega_\text{obs}}\)
redshift – relativistic redshift of the source, \(1+z\)
- Returns:
result
-
double SBody::FluxDensity(double spectral_density, double magnification)
Calculate the flux density. The result is.
\[F_{\nu,\text{obs}}=\frac{L_{\nu,\text{src}}}{1+z}\frac{\Omega_\text{src}}{\Omega_\text{obs}}.\]- Parameters:
spectral_density – intrinsic spectral density of the source, \(L_{\nu,\text{src}}\)
magnification – magnification of the source, including the relativistic redshift and beaming, and gravitational lensing. \(\frac{\Omega_\text{src}}{(1+z)\Omega_\text{obs}}\)
- Returns:
result
Elliptic Integrals
Calculate the elliptic integrals using the Legendre forms. Further information can be found in GSL Manual, Carlson (1988), Carlson (1989), Carlson (1991), and Carlson (1992).
-
double SBody::EllipticIntegral(int p5, double y, double x, double a5, double b5, double a1, double b1, double a2, double b2, double a3, double b3, double a4 = 1., double b4 = 0.)
- \[\int_y^x(a_5+b_5t)^{p_5/2}\prod_{i=1}^4(a_i+b_it)^{-1/2}dt\].
- Returns:
result
-
double SBody::EllipticIntegral2Complex(int p5, double y, double x, double a5, double b5, double f, double g, double h, double a1, double b1, double a4 = 1., double b4 = 0.)
- \[\int_y^x(a_5+b_5t)^{p_5/2}(f+gt+ht^2)^{-1/2}\prod_{i=1,4}(a_i+b_it)^{-1/2}dt\].
- Returns:
result
-
double SBody::EllipticIntegral4Complex(int p5, double y, double x, double a5, double b5, double f1, double g1, double h1, double f2, double g2, double h2)
- \[\int_y^x(a_5+b_5t)^{p_5/2}\prod_{i=1}^2(f_i+g_it+h_it^2)^{-1/2}dt\].
- Returns:
result
-
double SBody::Carlson_RC(double x, double y, gsl_mode_t mode = GSL_PREC_DOUBLE)
-
double SBody::Carlson_RJ(double x, double y, double z, double p, gsl_mode_t mode = GSL_PREC_DOUBLE)