1.3.1
CCC
 
MathFunc.cpp File Reference
#include "geoslib_old_f.h"
#include "geoslib_define.h"
#include "Matrix/MatrixRectangular.hpp"
#include "Basic/MathFunc.hpp"
#include "Basic/Law.hpp"
#include "Basic/WarningMacro.hpp"
#include "Basic/VectorHelper.hpp"
#include <math.h>
#include <boost/math/special_functions/legendre.hpp>
#include <boost/math/special_functions/spherical_harmonic.hpp>
#include <Geometry/GeometryHelper.hpp>

Classes

struct  Reg_Coor
 

Macros

#define COORD(i, ip)   (coord[3 * (ip) + (i)])
 
#define RCOORD(i, ip)   (R_coor->coor[3 * (ip) + (i)])
 
#define w   ((double *)&equiv_99)
 
#define x   ((double *)&equiv_100)
 
#define X   0.525731112119133696
 
#define Z   0.850650808352039932
 

Functions

void mvndst (int n, double *lower, double *upper, int *infin, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
 
void mvndst4 (double *lower, double *upper, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
 
void mvndst2n (double *lower, double *upper, double *means, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
 
int mvndst_infin (double low, double sup)
 
double bessel_j (double x, int n)
 
int bessel_j_table (double x, double alpha, int nb, double *b)
 
int bessel_k (double x, double alpha, int nb, double *bk)
 
double loggamma (double parameter)
 
double ut_legendre (int n, double v, bool flagNorm)
 
VectorDouble ut_legendreVec (int n, const VectorDouble &vecin, bool flagNorm)
 
MatrixRectangular ut_legendreMatNorm (int n, const VectorDouble &v)
 
MatrixRectangular ut_legendreAssociatedMat (int l, const VectorDouble &v, bool flagNorm)
 
double ut_flegendre (int n, int k0, double theta, bool flagNorm)
 
double ut_sphericalHarmonic (int n, int k, double theta, double phi)
 
VectorDouble ut_sphericalHarmonicVec (int n, int k, VectorDouble theta, VectorDouble phi)
 
double golden_search (double(*func_evaluate)(double test, void *user_data), void *user_data, double tolstop, double a0, double c0, double *test_loc, double *niter)
 
int ut_chebychev_count (double(*func)(double, double, const VectorDouble &), Cheb_Elem *cheb_elem, double x, const VectorDouble &blin)
 
int ut_chebychev_coeffs (double(*func)(double, double, const VectorDouble &), Cheb_Elem *cheb_elem, const VectorDouble &blin)
 
void ut_vandercorput (int n, int flag_sym, int flag_rot, int *ntri_arg, double **coor_arg)
 
void st_subdivide (double v1[3], double v2[3], double v3[3], int depth, Reg_Coor *R_coor)
 
int ut_icosphere (int n, int flag_rot, int *ntri_arg, double **coor_arg)
 
void ut_log_factorial (int nbpoly, double *factor)
 
double ut_factorial (int k)
 
DISABLE_WARNING_POP MatrixRectangularvanDerCorput (int n, int nd)
 
MatrixRectangular fillLegendreMatrix (const VectorDouble &r, int legendreOrder)
 

Macro Definition Documentation

◆ COORD

#define COORD (   i,
  ip 
)    (coord[3 * (ip) + (i)])

◆ RCOORD

#define RCOORD (   i,
  ip 
)    (R_coor->coor[3 * (ip) + (i)])

◆ w

#define w   ((double *)&equiv_99)

◆ x

#define x   ((double *)&equiv_100)

◆ X

#define X   0.525731112119133696

◆ Z

#define Z   0.850650808352039932

Function Documentation

◆ bessel_j()

double bessel_j ( double  x,
int  n 
)

◆ bessel_j_table()

int bessel_j_table ( double  x,
double  alpha,
int  nb,
double *  b 
)

This routine calculates Bessel functions J SUB(NB+ALPHA) (X) for non-negative argument X, and non-negative order NB+ALPHA.

Returns
Error return code : NCALC
NCALC < -1: An argument is out of range.
1 < NCALC < NB: Not all requested function values could be calculated accurately.
BY(I) contains correct function values
Parameters
[in]xWorking precision non-negative real argument for which J's are to be calculated.
[in]alphaWorking precision fractional part of order for which J's are to be calculated. 0 <= ALPHA < 1.0.
[in]nbInteger number of functions to be calculated, NB > 0 The first function calculated is of order ALPHA, and the last is of order (NB - 1 + ALPHA).
[out]bWorking precision output vector of length NB. If the routine terminates normally (NCALC=NB), the vector b[] contains the functions Y(ALPHA,X), ... ,Y(NB-1+ALPHA,X),
Remarks
This program is based on a program written by David J. Sookne (2)
that computes values of the Bessel functions J or I of real
argument and integer order. Modifications include the
restriction of the computation to the J Bessel function of
non-negative real argument, the extension of the computation to
arbitrary positive order, and the elimination of most underflow.
References: "A Note on Backward Recurrence Algorithms," Olver, F.
W. J., and Sookne, D. J., Math. Comp. 26, 1972, pp 941-947.
"Bessel Functions of Real Argument and Integer Order," Sookne, D.
J., NBS Jour. of Res. B. 77B, 1973, pp 125-132.

◆ bessel_k()

int bessel_k ( double  x,
double  alpha,
int  nb,
double *  bk 
)

This routine calculates modified Bessel functions of the second kind, K SUB(N+ALPHA) (X), for non-negative argument X and non-negative order N+ALPHA

Returns
Error return code :
NCALC < -1: An argument is out of range.
0 < NCALC < NB: Not all requested function values could be calculated accurately.
BK(I) contains correct function values for I <= NCALC, and contains the ratios
K(ALPHA+I-1,X)/K(ALPHA+I-2,X) for the rest of the array.
Parameters
[in]xWorking precision non-negative real argument for which K's are to calculated. If K's are to be calculated, X must not be greater than XMAX.
[in]alphaWorking precision fractional part of order for which K's are to be calculated. 0 <= ALPHA <1.0.
[in]nbInteger number of functions to be calculated, NB > 0. The first function calculated is of order ALPHA, and the last is of order (NB - 1 + ALPHA).
[out]bkWorking precision output vector of length NB. If the routine terminates normally (NCALC=NB), the vector BK contains the functions : K(ALPHA,X), ... , K(NB-1+ALPHA,X),
Remarks
This program is based on a program written by J. B. Campbell (2)
that computes values of the Bessel functions K of real argument
and real order. References: "On Temme's Algorithm for the \remark Modified Bessel Functions of the Third Kind," Campbell, J. B.,
TOMS 6(4), Dec. 1980, pp. 581-586. "A FORTRAN IV Subroutine for
the Modified Bessel Functions of the Third Kind of Real Order and
Real Argument," Campbell, J. B., Report NRC/ERB-925, National
Research Council, Canada.

◆ fillLegendreMatrix()

MatrixRectangular fillLegendreMatrix ( const VectorDouble r,
int  legendreOrder 
)

◆ golden_search()

double golden_search ( double(*)(double test, void *user_data)  func_evaluate,
void *  user_data,
double  tolstop,
double  a0,
double  c0,
double *  test_loc,
double *  niter 
)

Golden Search algorithm

Returns
The estimated value
Parameters
[in]func_evaluateEvaluating function
[in]user_dataUser Data
[in]tolstopTolerance parameter
[in]a0Initial value for lower bound of interval
[in]c0Initial value for upper bound of interval
[out]test_locFinal value of the evaluating function
[out]niterNumber of iterations

◆ loggamma()

double loggamma ( double  parameter)

Calculation of the logarithm of the gamma function

Returns
Logarithm of the gamma function
Parameters
[in]parameterraw value

◆ mvndst()

void mvndst ( int  n,
double *  lower,
double *  upper,
int *  infin,
double *  correl,
int  maxpts,
double  abseps,
double  releps,
double *  error,
double *  value,
int *  inform 
)

Multivariate Normal Probability

Parameters
[in]nNumber of variables
[in]lowerArray of lower integration limits
[in]upperArray of upper integration limits
[in]infinArray of integration limit flags
  • <0 for ]-Infinity; +Infinity[
  • =0 for ]-Infinity; upper]
  • =1 for [lower; +Infinity[
  • =2 for [lower; upper]
[in]correlArray of correlation coefficients
[in]maxptsMaximum number of function values allowed
[in]absepsAbsolute error tolerance
[in]relepsRelative error tolerance
[out]errorEstimated absolute error with 90% confidence level
[out]valueEstimated value for the integral
[out]informReturned code
Remarks
The array correl must be entered as the non-diagonal upper part
of the correlation matrix, entered by line
This subroutine uses an algorithm given in the paper:
"Numerical Computation of Multivariate Normal Probabilities", in
J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
Alan Genz
Department of Mathematics
Washington State University
Pullman, WA 99164-3113
Email : AlanG.nosp@m.enz@.nosp@m.wsu.e.nosp@m.du

◆ mvndst2n()

void mvndst2n ( double *  lower,
double *  upper,
double *  means,
double *  correl,
int  maxpts,
double  abseps,
double  releps,
double *  error,
double *  value,
int *  inform 
)

Calculate the multigaussian integral (non-normalized)

Parameters
[in]lowerArray of lower bounds (Dimension: nvar)
[in]upperArray of upper bounds (Dimension: nvar)
[in]meansArray of means (Dimension: 2)
[in]correlCorrelation matrix (Dimension: 2*2)
[in]maxptsMaximum number of function values allowed
[in]absepsAbsolute error tolerance
[in]relepsRelative error tolerance
[out]errorEstimated absolute error with 90% confidence level
[out]valueEstimated value for the integral
[out]informReturned code

◆ mvndst4()

void mvndst4 ( double *  lower,
double *  upper,
double *  correl,
int  maxpts,
double  abseps,
double  releps,
double *  error,
double *  value,
int *  inform 
)

Calculate the quadri-variable gaussian integral

Parameters
[in]lowerArray of lower bounds
[in]upperArray of upper bounds
[in]correlCorrelation matrix (Dimension: 4*4)
[in]maxptsMaximum number of function values allowed
[in]absepsAbsolute error tolerance
[in]relepsRelative error tolerance
[out]errorEstimated absolute error with 90% confidence level
[out]valueEstimated value for the integral
[out]informReturned code

◆ mvndst_infin()

int mvndst_infin ( double  low,
double  sup 
)

Set the flags for the bound of numerical integration

Returns
Flag for the bound
Parameters
[in]lowLower integration bound
[in]supUpper integration bound

◆ st_subdivide()

void st_subdivide ( double  v1[3],
double  v2[3],
double  v3[3],
int  depth,
Reg_Coor R_coor 
)

◆ ut_chebychev_coeffs()

int ut_chebychev_coeffs ( double(*)(double, double, const VectorDouble &)  func,
Cheb_Elem cheb_elem,
const VectorDouble blin 
)

Calculates the coefficients of the Chebychev polynomial which is an approximation of a given function

Returns
Error return code
Parameters
[in]funcFunction to be approximated
[in]cheb_elemCheb_Elem structure
[in]blinArray of coefficients for polynomial expansion

◆ ut_chebychev_count()

int ut_chebychev_count ( double(*)(double, double, const VectorDouble &)  func,
Cheb_Elem cheb_elem,
double  x,
const VectorDouble blin 
)

Evaluate the number of coefficients necessary to evaluate a function (at a sample location) at a given approximation

Returns
Minimum number of necessary coefficients
Parameters
[in]funcFunction to be approximated
[in]cheb_elemCheb_Elem structure
[in]xSampling value
[in]blinArray of coefficients for polynomial expansion

◆ ut_factorial()

double ut_factorial ( int  k)

Calculates the factorial coefficient

Returns
Returned value
Parameters
[in]kValue

◆ ut_flegendre()

double ut_flegendre ( int  n,
int  k0,
double  theta,
bool  flagNorm 
)

Returns the Spherical Legendre normalized function. According to boost Library documentation, this returns

Y_n^m(theta, phi) = sqrt{{(2n+1)(n-m)!} / {4pi(n+m)!}} P_n^m(cos(theta))e{imphi}

with phi=0 and m = |k0|

If flagNorm is switched ON, the previous result is multiplied by:

sqrt(2 * pi)

Parameters
[in]nDegree
[in]k0Order (ABS(k0) <= n)
[in]thetaTheta angle in radian
[in]flagNormfor normalized and 0 otherwise

◆ ut_icosphere()

int ut_icosphere ( int  n,
int  flag_rot,
int *  ntri_arg,
double **  coor_arg 
)

Generate regular Icosahedron discretization

Returns
Error return code
Parameters
[in]nNumber of discretization steps
[in]flag_rotPerform a random rotation
[out]ntri_argNumber of points
[out]coor_argArray of point coordinates (Dimension: 3*ntri)
Remarks
As random number are used in this function, a specific seed
is fixed here

◆ ut_legendre()

double ut_legendre ( int  n,
double  v,
bool  flagNorm 
)

Returns the Legendre Function described as follows:

legendre_p(n,v) = 1/{2^n n!) d^n/dx^n [x^2-1)^n with |x|<=1

If flagNorm is switched ON, the Boost library cannot be used anymore. We have to rely on the ode provided by X; Freulon.

Parameters
[in]nDegree of the Legendre Polynomial to be computed (n >= 0)
[in]vValue for which the polynomial is evaluated (-1 <= v <= 1)
[in]flagNormTrue for normalized and 0 otherwise
Returns
Value of the Legendre polynomial.

◆ ut_legendreAssociatedMat()

MatrixRectangular ut_legendreAssociatedMat ( int  l,
const VectorDouble v,
bool  flagNorm 
)

Returns the Associated Legendre Function described as follows:

In the case flagNorm=true

Using the relations: P_0^0 (x) = 1 P_{l+1}^{l+1} (x) = - sqrt((2l+3)/(2l+2)) * (1-x^2)^{1/2} * P_{l}^{l}(x) P_{l+1}^{m} (x) = sqrt((2l+3)(2l+1)/((l-m+1)(l+m+1))) * x * P_{l}^{m}(x) - sqrt((2l+3)/(2l-1)*(l+m)/(l+m+1)*(l-m)/(l-m+1)) P_{l-1}^m(x) Changing l to l-1, P_{l}^{l} (x) = - a0 * (1-x^2)^{1/2} * P_{l-1}^{l-1}(x) with a0 = sqrt((2l+1)/(2l)) P_{l}^{m} (x) = a * x * P_{l-1}^{m}(x) - b P_{l-2}^m(x) with a = sqrt((2l+1)(2l-1)/(l-m)/(l+m)) and b = sqrt((2l+1)/(2l-3)*(l+m-1)/(l+m)*(l-m-1)/(l-m))

In the case flagNorm=false

Using the relations: P_0^0 (x) = 1 P_{l+1}^{l+1} (x) = - sqrt((2l+3)/(2l+2)) * (1-x^2)^{1/2} * P_{l}^{l}(x) P_{l+1}^{m} (x) = sqrt((2l+3)(2l+1)/((l-m+1)(l+m+1))) * x * P_{l}^{m}(x) - sqrt((2l+3)/(2l-1)*(l+m)/(l+m+1)*(l-m)/(l-m+1)) P_{l-1}^m(x) Changing l to l-1, P_{l}^{l} (x) = - a0 * (1-x^2)^{1/2} * P_{l-1}^{l-1}(x) with a0 = sqrt((2l+1)/(2l)) P_{l}^{m} (x) = a * x * P_{l-1}^{m}(x) - b P_{l-2}^m(x) with a = sqrt((2l+1)(2l-1)/(l-m)/(l+m)) and b = sqrt((2l+1)/(2l-3)*(l+m-1)/(l+m)*(l-m-1)/(l-m))

Parameters
[in]lDegree of the Legendre Polynomial to be computed (n >= 0)
[in]vValue for which the polynomial is evaluated (-1 <= v <= 1)
[in]flagNormNormalized when True
Returns
A matrix of the (normalized) associated Legendre functions: P_l^m (x) for 0 <= m <= l

◆ ut_legendreMatNorm()

MatrixRectangular ut_legendreMatNorm ( int  n,
const VectorDouble v 
)

◆ ut_legendreVec()

VectorDouble ut_legendreVec ( int  n,
const VectorDouble vecin,
bool  flagNorm 
)

◆ ut_log_factorial()

void ut_log_factorial ( int  nbpoly,
double *  factor 
)

Calculates the nbpoly log-factorial coefficients

Parameters
[in]nbpolyNumber of terms
[out]factorlogarithm of factorials

◆ ut_sphericalHarmonic()

double ut_sphericalHarmonic ( int  n,
int  k,
double  theta,
double  phi 
)

Returns the Spherical harmonic

Parameters
[in]nDegree of HS (n >= 0)
[in]kOrder of the HS (-n <= k <= n)
[in]thetaColatitude angle in radian (0 <= theta <= pi
[in]phiLongitude angle in radian (0 <= phi <= 2* pi)

◆ ut_sphericalHarmonicVec()

VectorDouble ut_sphericalHarmonicVec ( int  n,
int  k,
VectorDouble  theta,
VectorDouble  phi 
)

◆ ut_vandercorput()

void ut_vandercorput ( int  n,
int  flag_sym,
int  flag_rot,
int *  ntri_arg,
double **  coor_arg 
)

Generate a Van Der Corput list of points in R^3

Parameters
[in]nNumber of points
[in]flag_symDuplicate the samples by symmetry
[in]flag_rotPerform a random rotation
[out]ntri_argNumber of points
[out]coor_argArray of point coordinates (Dimension: 3*ntri)

◆ vanDerCorput()

DISABLE_WARNING_POP MatrixRectangular* vanDerCorput ( int  n,
int  nd 
)

Function to compute the vector Van der Corput sequence or the Halton sequence

Parameters
nThe number of values to be computed
ndThe dimension of output sequence
Returns
A matrix of dimensions [n,nd] with the sequence values (between 0 and 1)
Note
The dimension nd should be lower or equal to 50.