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 MatrixRectangular *	vanDerCorput (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]	x	Working precision non-negative real argument for which J's are to be calculated.
[in]	alpha	Working precision fractional part of order for which J's are to be calculated. 0 <= ALPHA < 1.0.
[in]	nb	Integer 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]	b	Working 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]	x	Working 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]	alpha	Working precision fractional part of order for which K's are to be calculated. 0 <= ALPHA <1.0.
[in]	nb	Integer 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]	bk	Working 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_evaluate	Evaluating function
[in]	user_data	User Data
[in]	tolstop	Tolerance parameter
[in]	a0	Initial value for lower bound of interval
[in]	c0	Initial value for upper bound of interval
[out]	test_loc	Final value of the evaluating function
[out]	niter	Number of iterations

loggamma()

double loggamma

(

double

parameter

)

Calculation of the logarithm of the gamma function

Returns

Logarithm of the gamma function

Parameters

[in]

parameter

raw 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]	n	Number of variables
[in]	lower	Array of lower integration limits
[in]	upper	Array of upper integration limits
[in]	infin	Array of integration limit flags <0 for ]-Infinity; +Infinity[ =0 for ]-Infinity; upper] =1 for [lower; +Infinity[ =2 for [lower; upper]
[in]	correl	Array of correlation coefficients
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned 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]	lower	Array of lower bounds (Dimension: nvar)
[in]	upper	Array of upper bounds (Dimension: nvar)
[in]	means	Array of means (Dimension: 2)
[in]	correl	Correlation matrix (Dimension: 2*2)
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned 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]	lower	Array of lower bounds
[in]	upper	Array of upper bounds
[in]	correl	Correlation matrix (Dimension: 4*4)
[in]	maxpts	Maximum number of function values allowed
[in]	abseps	Absolute error tolerance
[in]	releps	Relative error tolerance
[out]	error	Estimated absolute error with 90% confidence level
[out]	value	Estimated value for the integral
[out]	inform	Returned 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]	low	Lower integration bound
[in]	sup	Upper 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]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	blin	Array 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]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	x	Sampling value
[in]	blin	Array of coefficients for polynomial expansion

ut_factorial()

double ut_factorial

(

int

k

)

Calculates the factorial coefficient

Returns

Returned value

Parameters

[in]

k

Value

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]	n	Degree
[in]	k0	Order (ABS(k0) <= n)
[in]	theta	Theta angle in radian
[in]	flagNorm	for 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]	n	Number of discretization steps
[in]	flag_rot	Perform a random rotation
[out]	ntri_arg	Number of points
[out]	coor_arg	Array 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]	n	Degree of the Legendre Polynomial to be computed (n >= 0)
[in]	v	Value for which the polynomial is evaluated (-1 <= v <= 1)
[in]	flagNorm	True 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]	l	Degree of the Legendre Polynomial to be computed (n >= 0)
[in]	v	Value for which the polynomial is evaluated (-1 <= v <= 1)
[in]	flagNorm	Normalized 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]	nbpoly	Number of terms
[out]	factor	logarithm of factorials

ut_sphericalHarmonic()

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

Returns the Spherical harmonic

Parameters

[in]	n	Degree of HS (n >= 0)
[in]	k	Order of the HS (-n <= k <= n)
[in]	theta	Colatitude angle in radian (0 <= theta <= pi
[in]	phi	Longitude 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]	n	Number of points
[in]	flag_sym	Duplicate the samples by symmetry
[in]	flag_rot	Perform a random rotation
[out]	ntri_arg	Number of points
[out]	coor_arg	Array 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

n	The number of values to be computed
nd	The 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.