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 MatrixRectangular *	vanDerCorput (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]	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 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]	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.

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_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

double loggamma

(

double

parameter

)

Calculation of the logarithm of the gamma function

Returns

Logarithm of the gamma function

Parameters

[in]

parameter

raw 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]	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

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

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

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

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]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	blin	Array 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]	func	Function to be approximated
[in]	cheb_elem	Cheb_Elem structure
[in]	x	Sampling value
[in]	blin	Array of coefficients for polynomial expansion

double ut_factorial

(

int

k

)

Calculates the factorial coefficient

Returns

Returned value

Parameters

[in]

k

Value

double ut_flegendre	(	int	flag_norm,
		int	n,
		int	k0,
		double	theta
	)

Returns the Spherical Legendre normalized function

Parameters

[in]	flag_norm	1 for normalized and 0 otherwise
[in]	n	Degree
[in]	k0	Order (ABS(k0) <= n)
[in]	theta	Theta 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]	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

double ut_legendre	(	int	flag_norm,
		int	n,
		double	v
	)

Returns the Associated Legendre Function

Parameters

[in]	flag_norm	1 for normalized and 0 otherwise
[in]	n	Degree
[in]	v	Value

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

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)

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.