MathFunc.cpp File Reference

#include "geoslib_define.h"
#include "geoslib_d.h"
#include "Matrix/MatrixRectangular.hpp"
#include "Basic/MathFunc.hpp"
#include "Basic/Law.hpp"
#include "Basic/WarningMacro.hpp"
#include "Core/fftn.hpp"
#include "Basic/Memory.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
static VectorDouble	_corputVector (int n, int b)
static double	st_mvnphi (const double *z)
static void	st_mvnlms (double *a, double *b, const int *infin, double *lower, double *upper)
static void	st_dkswap (double *x, double *y)
static void	st_rcswp (const int *p, const int *q, double *a, double *b, int *infin, const int *n, double *c)
static void	st_covsrt (int *n, double *lower, double *upper, double *correl, int *infin, double *y, int *infis, double *a, double *b, double *cov, int *infi)
static double	st_phinvs (const double *p)
static double	st_bvu (const double *sh, const double *sk, const double *r)
static double	st_bvnmvn (double *lower, double *upper, int *infin, double *correl)
static double	st_mvndfn_0 (int n__, int *n, double *w, double *correl, double *lower, double *upper, int *infin, int *infis, double *d, double *e)
static double	st_mvndnt (int *n, double *correl, double *lower, double *upper, int *infin, int *infis, double *d, double *e)
static void	st_dkrcht (const int *s, double *quasi)
static void	st_dksmrc (int *ndim, const int *klim, double *sumkro, const int *prime, double *vk, double(*functn)(int *, double *), double *x)
static void	st_dkbvrc (int *ndim, int *minvls, const int *maxvls, double(*functn)(int *, double *), const double *abseps, const double *releps, double *abserr, double *finest, int *inform)
static double	st_mvndfn (int *n, double *w)
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, const double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
void	mvndst2n (const double *lower, const double *upper, const double *means, double *correl, int maxpts, double abseps, double releps, double *error, double *value, int *inform)
int	mvndst_infin (double low, double sup)
double	besselj (double x, int n)
int	besselj_table (double x, double alpha, int nb, double *b)
int	besselk (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)
static void	st_init_rotation (double *ct, double *st, double *a)
void	ut_vandercorput (int n, int flag_sym, int flag_rot, int *ntri_arg, double **coor_arg)
static void	st_addTriangle (const double v1[3], const double v2[3], const double v3[3], Reg_Coor *R_coor)
static void	st_normalize (double v[3])
void	st_subdivide (double v1[3], double v2[3], double v3[3], int depth, Reg_Coor *R_coor)
static int	st_already_present (Reg_Coor *R_coor, int i0, int ntri, const double *coord, double eps=EPSILON3)
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)
int	solve_P2 (double a, double b, double c, VectorDouble &x)
int	solve_P3 (double a, double b, double c, double d, VectorDouble &x)

Variables
static double	c_b11 = 1.

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

_corputVector()

static VectorDouble _corputVector ( int  n, int  b  )

static

Function to compute the Van der Corput sequence

Parameters

n	The number of values to be computed
b	The base in which the numbers are represented

Returns

A vector of first n values of the the Van Der Corput sequence

Note

The base should be a prime number

besselj()

double besselj	(	double	x,
		int	n
	)

besselj_table()

int besselj_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.

besselk()

int besselk	(	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	(	const double *	lower,
		const double *	upper,
		const 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,
		const 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

solve_P2()

int solve_P2	(	double	a,
		double	b,
		double	c,
		VectorDouble &	x
	)

Find the roots of a polynomial of order 2: ax^2 + bx + c = 0

Returns

Number of real solutions

Parameters

[in]	a,b,c	Coefficients of the polynomial
[out]	x	Array of real solutions (Dimension: 2)

Remarks

When the solution is double, the returned number os 1.

solve_P3()

int solve_P3	(	double	a,
		double	b,
		double	c,
		double	d,
		VectorDouble &	x
	)

Find the roots of a polynomial of order 3: a*x^3 + b*x^2 + c*x + d = 0

Returns

Number of real solutions

Parameters

[in]	a,b,c,d	Coefficients of the polynomial
[out]	x	Array of real solutions (Dimension: 3)

Remarks

When the solution is double, the returned number os 1.

st_addTriangle()

static void st_addTriangle ( const double  v1[3], const double  v2[3], const double  v3[3], Reg_Coor *  R_coor  )

static

st_already_present()

static int st_already_present ( Reg_Coor *  R_coor, int  i0, int  ntri, const double *  coord, double  eps = EPSILON3  )

static

st_bvnmvn()

static double st_bvnmvn ( double *  lower, double *  upper, int *  infin, double *  correl  )

static

A function for computing bivariate normal probabilities

Returns

bivariate normal probability

Parameters

[in]	lower	array of lower integration limits
[in]	upper	array of upper integration limits
[in]	infin	array of integration limits flag (0,1,2)
[in]	correl	correlation coefficient

st_bvu()

static double st_bvu ( const double *  sh, const double *  sk, const double *  r  )

static

A function for computing bivariate normal probabilities. Yihong Ge Department of Computer Science and Electrical Engineering Washington State University Pullman, WA 99164-2752 and Alan Genz Department of Mathematics Washington State University Pullman, WA 99164-3113 Email : alang.nosp@m.enz@.nosp@m.wsu.e.nosp@m.du

Returns

Calculate the probability that X is larger than SH and Y is

larger than SK.

Parameters

[in]	sh	integration limit
[in]	sk	integration limit
[in]	r	REAL, correlation coefficient

st_covsrt()

static void st_covsrt ( int *  n, double *  lower, double *  upper, double *  correl, int *  infin, double *  y, int *  infis, double *  a, double *  b, double *  cov, int *  infi  )

static

Subroutine to sort integration limits and determine Cholesky factor

st_dkbvrc()

static void st_dkbvrc ( int *  ndim, int *  minvls, const int *  maxvls, double(*)(int *, double *)  functn, const double *  abseps, const double *  releps, double *  abserr, double *  finest, int *  inform  )

static

Automatic Multidimensional Integration Subroutine AUTHOR: Alan Genz Department of Mathematics Washington State University Pulman, WA 99164-3113 Email: AlanG.nosp@m.enz@.nosp@m.wsu.e.nosp@m.du Last Change: 5/15/98 KRBVRC computes an approximation to the integral 1 1 1 I I ... I F(X) dx(NDIM)...dx(2)dx(1) 0 0 0 DKBVRC uses randomized Korobov rules for the first 20 variables. The primary references are "Randomization of Number Theoretic Methods for Multiple Integration"

• R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13, pp. 904-14,
• and "Optimal Parameters for Multidimensional Integration", P. Keast, SIAM J Numer Anal, 10, pp.831-838. If there are more than 20 variables, the remaining variables are integrated using Richtmeyer rules. A reference is "Methods of Numerical Integration", P.J. Davis and P. Rabinowitz, Academic Press, 1984, pp. 482-483.

Parameters

[in]	ndim	Number of variables, must exceed 1, but not exceed 40 Integer minimum number of function evaluations allowed.
[in,out]	minvls	must not exceed MAXVLS. If MINVLS < 0 then the routine assumes a previous call has been made with the same integrand and continues that calculation. Actual number of function evaluations used.
[in]	maxvls	Integer maximum number of function evaluations allowed.
[in]	functn	EXTERNALLY declared user defined function to be integrated. It must have parameters (NDIM,Z), where Z is a real array of dimension NDIM.
[in]	abseps	Required absolute accuracy. Estimated absolute accuracy of FINEST.
[in]	releps	Required relative accuracy.
[out]	abserr	Absolute returned accuracy
[out]	finest	Estimated value of integral.
[out]	inform	= 0 for normal status, when ABSERR <= MAX(ABSEPS, RELEPS*ABS(FINEST)) and INTVLS <= MAXCLS. = 1 If MAXVLS was too small to obtain the required accuracy. In this case a value FINEST is returned with estimated absolute accuracy ABSERR.

st_dkrcht()

static void st_dkrcht ( const int *  s, double *  quasi  )

static

This subroutine generates a new quasi-random Richtmeyer vector. A reference is "Methods of Numerical Integration", P.J. Davis and P. Rabinowitz, Academic Press, 1984, pp. 482-483.

Parameters

[in]	s	the number of dimensions
[out]	quasi	a new quasi-random S-vector

st_dksmrc()

static void st_dksmrc ( int *  ndim, const int *  klim, double *  sumkro, const int *  prime, double *  vk, double(*)(int *, double *)  functn, double *  x  )

static

st_dksmrc

st_dkswap()

static void st_dkswap ( double *  x, double *  y  )

static

st_dkswap

st_init_rotation()

static void st_init_rotation ( double *  ct, double *  st, double *  a  )

static

Generate a random rotation axis

Parameters

[out]	ct	Cosine of the random rotation angle
[out]	st	Sine of the random rotation angle
[out]	a	Random direction vector

st_mvndfn()

static double st_mvndfn ( int *  n, double *  w  )

static

st_mvndfn

st_mvndfn_0()

static double st_mvndfn_0 ( int  n__, int *  n, double *  w, double *  correl, double *  lower, double *  upper, int *  infin, int *  infis, double *  d, double *  e  )

static

Integrand subroutine

st_mvndnt()

static double st_mvndnt ( int *  n, double *  correl, double *  lower, double *  upper, int *  infin, int *  infis, double *  d, double *  e  )

static

st_mvndnt

st_mvnlms()

static void st_mvnlms ( double *  a, double *  b, const int *  infin, double *  lower, double *  upper  )

static

Multivariate Normal Probability (local function)

st_mvnphi()

static double st_mvnphi ( const double *  z )

static

Normal distribution probabilities accurate to 1.e-15. Z = no. of standard deviations from the mean. Based upon algorithm 5666 for the error function, from: Hart, J.F. et al, 'Computer Approximations', Wiley 1968

st_normalize()

static void st_normalize ( double  v[3] )

static

st_phinvs()

static double st_phinvs ( const double *  p )

static

Produces the normal deviate Z corresponding to a given lower tail area of P. The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

st_rcswp()

static void st_rcswp ( const int *  p, const int *  q, double *  a, double *  b, int *  infin, const int *  n, double *  c  )

static

Swaps rows and columns P and Q in situ, with P <= Q

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.

Variable Documentation

c_b11

double c_b11 = 1.

static