1.2.0
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 <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)
 
int bessel_j (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 flag_norm, int n, double v)
 
double ut_flegendre (int flag_norm, int n, int k0, double theta)
 
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

#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

Function Documentation

int bessel_j ( 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 by[] 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.
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 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.
MatrixRectangular fillLegendreMatrix ( const VectorDouble r,
int  legendreOrder 
)
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
double loggamma ( double  parameter)

Calculation of the logarithm of the gamma function

Returns
Logarithm of the gamma function
Parameters
[in]parameterraw value
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
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
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
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
void st_subdivide ( double  v1[3],
double  v2[3],
double  v3[3],
int  depth,
Reg_Coor R_coor 
)
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
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
double ut_factorial ( int  k)

Calculates the factorial coefficient

Returns
Returned value
Parameters
[in]kValue
double ut_flegendre ( int  flag_norm,
int  n,
int  k0,
double  theta 
)

Returns the Spherical Legendre normalized function

Parameters
[in]flag_norm1 for normalized and 0 otherwise
[in]nDegree
[in]k0Order (ABS(k0) <= n)
[in]thetaTheta angle in radian
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
double ut_legendre ( int  flag_norm,
int  n,
double  v 
)

Returns the Associated Legendre Function

Parameters
[in]flag_norm1 for normalized and 0 otherwise
[in]nDegree
[in]vValue
void ut_log_factorial ( int  nbpoly,
double *  factor 
)

Calculates the nbpoly log-factorial coefficients

Parameters
[in]nbpolyNumber of terms
[out]factorlogarithm of factorials
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)
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.