Namespace for some extra maths function that are often needed. More...

Functions
Mdouble	gamma (Mdouble gamma_in)
	This is the gamma function returns the true value for the half integer value. More...

Mdouble	chi_squared (const Mdouble x, const unsigned int k)
	This is a chi_squared function return the value x and degrees of freedom k. More...

Mdouble	chi_squared_prob (const Mdouble x, const unsigned int k)
	This is the function which actually gives the probability back using a chi squared test. More...

Mdouble	goldenSectionSearch (Mdouble(*function)(const Mdouble), Mdouble min, Mdouble cur, Mdouble max, Mdouble endCondition, Mdouble curVal=std::numeric_limits< Mdouble >::quiet_NaN())
	This function performs a golden section search to find the location of the minimum of a function. More...

template<typename T >
int	sign (T val)
	This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0. More...

template<typename T >
T	square (T val)
	squares a number More...

template<typename T >
T	cubic (T val)
	calculates the cube of a number More...

bool	isEqual (Mdouble v1, Mdouble v2, Mdouble absError)
	Compares the difference of two Mdouble with an absolute error, useful in UnitTests. More...

bool	isEqual (Vec3D v1, Vec3D v2, Mdouble absError)
	Compares the difference of two Vec3D with an absolute error, useful in UnitTests. More...

template<typename T >
constexpr T	factorial (const T t)
	factorial function More...

Mdouble	sin (Mdouble x)

Mdouble	cos (Mdouble x)

Mdouble	exp (Mdouble Exponent)

Mdouble	log (Mdouble Power)

template<typename T >
T	tan (T x)

Mdouble	chebyshev (Mdouble x, const Mdouble coef[], int N)
	Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc. More...

Mdouble	I0_exp (Mdouble x)

Mdouble	I0 (Mdouble x)

Detailed Description

Namespace for some extra maths function that are often needed.

Function Documentation

Mdouble mathsFunc::chebyshev	(	Mdouble	x,
		const Mdouble	coef[],
		int	N
	)

Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc.

Definition at line 242 of file ExtendedMath.cc.

References assert.

Referenced by I0_exp().

 {
     const Mdouble *p = coef;
     Mdouble b0 = *p++;
     Mdouble b1 = 0, b2;
     int i = N - 1;
     
     assert(i > 0);
     do
     {
         b2 = b1;
         b1 = b0;
         b0 = x * b1 - b2 + *p++;
     }
     while ( --i);
     
     return (0.5 * (b0 - b2));
 }

Mdouble mathsFunc::chi_squared	(	const Mdouble	x,
		const unsigned int	k
	)

This is a chi_squared function return the value x and degrees of freedom k.

Definition at line 153 of file ExtendedMath.cc.

References exp(), and gamma().

Referenced by chi_squared_prob().

 {
     
     Mdouble prefactor = pow(2, k / 2.0) * gamma(k / 2.0);
     Mdouble mainfactor = pow(x, k / 2.0 - 1) * exp(x / -2.0);
     
     return mainfactor / prefactor;
     
 }

Mdouble mathsFunc::chi_squared_prob	(	const Mdouble	x_max,
		const unsigned int	k
	)

This is the function which actually gives the probability back using a chi squared test.

This calulates the probabity based on a chi squared test First we calculated the cummelative chi_squared function.

This is the function which actually gives the probability back It is calculated by calling the normal chi_squated function and using the trapezoidal rule. The final results is 1-the cummulative chi_squared function

Definition at line 169 of file ExtendedMath.cc.

References chi_squared().

Referenced by RNG::test().

 {
     
 //The current value was picked by tried were it stopped effect the 4 d.p.
     const int num_steps_per_unit = 100;
     Mdouble sum = 0;
     Mdouble x = 0;
     long int num_steps = static_cast<int>(num_steps_per_unit * x_max);
 //Use trapezional rule, but ignoring the ends
     for (int i = 0; i < num_steps; i++)
     {
         x = x_max / num_steps * (i + 0.5);
         sum = sum + chi_squared(x, k);
     }
     return 1.0 - sum * x_max / num_steps;
     
 }

Mdouble mathsFunc::cos ( Mdouble x )

Definition at line 60 of file ExtendedMath.cc.

References constants::pi.

Referenced by ChuteWithHopper::addHopper(), LiquidBridgeWilletInteraction::computeAdhesionForce(), LiquidMigrationWilletInteraction::computeAdhesionForce(), CircularPeriodicBoundary::createPeriodicParticles(), HopperInsertionBoundary::generateParticle(), Screw::getDistanceAndNormal(), Coil::getDistanceAndNormal(), InfiniteWallWithHole::getDistanceAndNormal(), Vec3D::getFromCylindricalCoordinates(), Screw::getTriangulation(), Quarternion::integrate(), helpers::objectivenessTest(), Quarternion::Quarternion(), CircularPeriodicBoundary::rotateParticle(), Chute::setChuteAngleAndMagnitudeOfGravity(), tan(), and AxisymmetricIntersectionOfWalls::writeVTK().

                                 {
     Mdouble N = floor(x/(2.0*constants::pi)+0.5);
     x -= N*2.0*constants::pi;
     Mdouble Sum = 1.0;
     Mdouble Power = 1;
     Mdouble Sign = 1;
     const Mdouble x2 = x * x;
     Mdouble Fact = 1.0;
     for (unsigned int i = 2; i < 25; i += 2) {
         Power *= x2;
         Fact *= i * (i - 1);
         Sign *= -1.0;
         Sum += Sign * Power / Fact;
     }
     return Sum;
 }

template<typename T >

T mathsFunc::cubic ( T val )

calculates the cube of a number

Definition at line 99 of file ExtendedMath.h.

Referenced by ChuteBottom::setupInitialConditions().

     {
         return val * val * val;
     }

Mdouble mathsFunc::exp ( Mdouble Exponent )

Definition at line 78 of file ExtendedMath.cc.

References X.

Referenced by chi_squared(), helpers::computeRestitutionCoefficientFromCollisionTimeAndDispAndEffectiveMass(), helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(), helpers::computeRestitutionCoefficientFromKAndDispAndEffectiveMass(), helpers::getMaximumVelocity(), LinearViscoelasticNormalSpecies::getRestitutionCoefficient(), HGridOptimiser::histNumberParticlesPerCell(), I0(), log(), ParhamiMcMeekingSinterSpecies::set(), and SinterNormalSpecies::setParhamiMcKeeping().

                                        {
     Mdouble X, P, Frac, I, L;
     X = Exponent;
     Frac = X;
     P = (1.0 + X);
     I = 1.0;
 
     do
     {
         I++;
         Frac *= (X / I);
         L = P;
         P += Frac;
     }while(L != P);
 
     return P;
 }

template<typename T >

constexpr T mathsFunc::factorial ( const T t )

factorial function

Definition at line 124 of file ExtendedMath.h.

Referenced by HGridOptimiser::histNumberParticlesPerCell().

     {
         return (t == 0) ? 1 : t * factorial(t - 1);
     }

Mdouble mathsFunc::gamma ( Mdouble gamma_in )

This is the gamma function returns the true value for the half integer value.

This is the gamma function, gives 'exact' answers for the half integer values This is done using the recussion relation and the known values for 1 and 0.5 Note, return NaN for non-half integer values.

Definition at line 130 of file ExtendedMath.cc.

References constants::sqrt_pi.

Referenced by chi_squared(), and HopperInsertionBoundary::generateParticle().

 {
     const Mdouble ep = 1e-5;
     
     if (gamma_in > 1.0 + ep)
     {
         return ((gamma_in - 1) * gamma(gamma_in - 1));
     }
     else
     {
         
         if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
             return 1.0;
         else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
             return constants::sqrt_pi;
         else
             return std::numeric_limits<Mdouble>::quiet_NaN();
     }
 } //end func gamma

Mdouble mathsFunc::goldenSectionSearch	(	Mdouble(*)(const Mdouble)	function,
		Mdouble	min,
		Mdouble	cur,
		Mdouble	max,
		Mdouble	endCondition,
		Mdouble	curVal = `std::numeric_limits<Mdouble>::quiet_NaN()`
	)

This function performs a golden section search to find the location of the minimum of a function.

Parameters

[in]	function	A function pointer to the function of which you want to calculate the location of its minimum.
[in]	min	The minimum location
[in]	cur	The current location
[in]	max	The maximum location
[in]	endCondition	The algorithm terminates when abs(max - min) < endCondition
[in]	curVal	The value of the function at the current location (on default this value is calculated internally)

Definition at line 187 of file ExtendedMath.cc.

 {
     if (std::abs(max - min) < endCondition)
     {
         return 0.5 * (min + max);
     }
     std::cout << "Min=" << min << " Max=" << max << " diff=" << max - min << std::endl;
     Mdouble resphi = 2 - 0.5 * (1 + std::sqrt(5));
     Mdouble x;
     if (max - cur > cur - min)
     {
         x = cur + resphi * (max - cur);
     }
     else
     {
         x = cur - resphi * (cur - min);
     }
     if(std::isnan(curVal))
         curVal = function(cur);
     Mdouble xVal = function(x);
     if (xVal < curVal)
     {
         if (max - cur > cur - min)
         {
             return goldenSectionSearch(function, cur, x, max, endCondition, xVal);
         }
         else
         {
             return goldenSectionSearch(function, min, x, cur, endCondition, xVal);
         }
     }
     else
     {
         if (max - cur > cur - min)
         {
             return goldenSectionSearch(function, min, cur, x, endCondition, curVal);
         }
         else
         {
             return goldenSectionSearch(function, x, cur, max, endCondition, curVal);
         }
     }
 }

Mdouble mathsFunc::I0 ( Mdouble x )

Definition at line 341 of file ExtendedMath.cc.

References exp(), and I0_exp().

 {
     if (x < 0)
         x = -x;
     return exp(x) * I0_exp(x);
 }

Mdouble mathsFunc::I0_exp ( Mdouble x )

Definition at line 261 of file ExtendedMath.cc.

References A, and chebyshev().

Referenced by I0().

 {
     // Cooefficients for [0..8]
     const Mdouble A[] =
             {
                     -4.415341646479339379501E-18,
                     3.330794518822238097831E-17,
                     -2.431279846547954693591E-16,
                     1.715391285555133030611E-15,
                     -1.168533287799345168081E-14,
                     7.676185498604935616881E-14,
                     -4.856446783111929460901E-13,
                     2.955052663129639834611E-12,
                     -1.726826291441555707231E-11,
                     9.675809035373236912241E-11,
                     -5.189795601635262906661E-10,
                     2.659823724682386650351E-9,
                     -1.300025009986248042121E-8,
                     6.046995022541918949321E-8,
                     -2.670793853940611733911E-7,
                     1.117387539120103718151E-6,
                     -4.416738358458750563591E-6,
                     1.644844807072889708931E-5,
                     -5.754195010082103703981E-5,
                     1.885028850958416557291E-4,
                     -5.763755745385823658851E-4,
                     1.639475616941335798421E-3,
                     -4.324309995050575944301E-3,
                     1.054646039459499831831E-2,
                     -2.373741480589946881561E-2,
                     4.930528423967070848781E-2,
                     -9.490109704804764442101E-2,
                     1.716209015222087753491E-1,
                     -3.046826723431983986831E-1,
                     6.767952744094760849951E-1
             };
     
     // Cooefficients for [8..infinity]
     const Mdouble B[] =
             {
                     -7.233180487874753954561E-18,
                     -4.830504485944182071261E-18,
                     4.465621420296759999011E-17,
                     3.461222867697461093101E-17,
                     -2.827623980516583484941E-16,
                     -3.425485619677219134621E-16,
                     1.772560133056526383601E-15,
                     3.811680669352622420751E-15,
                     -9.554846698828307648701E-15,
                     -4.150569347287222086631E-14,
                     1.540086217521409826911E-14,
                     3.852778382742142701141E-13,
                     7.180124451383666233671E-13,
                     -1.794178531506806117781E-12,
                     -1.321581184044771311881E-11,
                     -3.149916527963241364541E-11,
                     1.188914710784643834241E-11,
                     4.940602388224969589101E-10,
                     3.396232025708386345151E-9,
                     2.266668990498178064591E-8,
                     2.048918589469063741831E-7,
                     2.891370520834756482971E-6,
                     6.889758346916823984261E-5,
                     3.369116478255694089901E-3,
                     8.044904110141088316081E-1
             };
     
     if (x < 0)
         x = -x;
     
     if (x <= 8.0)
     {
         Mdouble y = (x / 2.0) - 2.0;
         return (chebyshev(y, A, 30));
     }
     
     return (chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
 
 }

bool mathsFunc::isEqual	(	Mdouble	v1,
		Mdouble	v2,
		Mdouble	absError
	)

Compares the difference of two Mdouble with an absolute error, useful in UnitTests.

Parameters

[in]	v1	The first Mdouble
[in]	v2	The second Mdouble
[in]	absError	The allowed maximum absolute error

Returns: True if the two Mdouble are equal

Definition at line 231 of file ExtendedMath.cc.

Referenced by helpers::check(), ChuteInsertionBoundary::generateParticle(), isEqual(), and MercuryBase::setHGridCellOverSizeRatio().

 {
     return std::abs(v1 - v2) <= absError;
 }

bool mathsFunc::isEqual	(	Vec3D	v1,
		Vec3D	v2,
		Mdouble	absError
	)

Compares the difference of two Vec3D with an absolute error, useful in UnitTests.

Parameters

[in]	v1	The first Vec3D
[in]	v2	The second Vec3D
[in]	absError	The allowed maximum absolute error

Returns: true if the two Vec3D are equal

Definition at line 237 of file ExtendedMath.cc.

References isEqual(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

 {
     return isEqual(v1.X,v2.X,absError)&&isEqual(v1.Y,v2.Y,absError)&&isEqual(v1.Z,v2.Z,absError);
 }

Mdouble mathsFunc::log ( Mdouble Power )

Todo:: check if this function works

Definition at line 97 of file ExtendedMath.cc.

References A, exp(), and R.

 {
     Mdouble N, P, L, R, A, E;
     E = 2.71828182845905;
     P = Power;
     N = 0.0;
 
     // This speeds up the convergence by calculating the integral
     while(P >= E)
     {
         P /= E;
         N++;
     }
     N += (P / E);
     P = Power;
     do
     {
         A = N;
         L = (P / (exp(N - 1.0)));
         R = ((N - 1.0) * E);
         N = ((L + R) / E);
     } while (N < A);
     //} while (N != A);
 
     return N;
 }

template<typename T >

int mathsFunc::sign ( T val )

This is a sign function, it returns -1 for negative numbers, 1 for positive numbers and 0 for 0.

Definition at line 83 of file ExtendedMath.h.

Referenced by Screw::getDistanceAndNormal().

     {
         return (T(0) < val) - (val < T(0));
     }

Mdouble mathsFunc::sin ( Mdouble x )

Definition at line 42 of file ExtendedMath.cc.

References constants::pi.

Referenced by ChuteWithHopper::addHopper(), HopperInsertionBoundary::generateParticle(), Quarternion::getAxisAndAngle(), Screw::getDistanceAndNormal(), Coil::getDistanceAndNormal(), InfiniteWallWithHole::getDistanceAndNormal(), Vec3D::getFromCylindricalCoordinates(), ChuteWithHopper::getMaximumVelocityInducedByGravity(), Screw::getTriangulation(), Quarternion::integrate(), helpers::objectivenessTest(), Quarternion::Quarternion(), CircularPeriodicBoundary::rotateParticle(), Chute::setChuteAngleAndMagnitudeOfGravity(), tan(), and AxisymmetricIntersectionOfWalls::writeVTK().

                                 {
     Mdouble N = floor(x/(2.0*constants::pi)+0.5);
     x -= N*2.0*constants::pi;
     Mdouble Sum = 0;
     Mdouble Power = x;
     Mdouble Sign = 1;
     const Mdouble x2 = x * x;
     Mdouble Fact = 1.0;
     for (unsigned int i = 1; i < 25; i += 2) {
         Sum += Sign * Power / Fact;
         Power *= x2;
         Fact *= (i + 1) * (i + 2);
         Sign *= -1.0;
     }
     return Sum;
 }

template<typename T >

T mathsFunc::square ( T val )

squares a number

Definition at line 91 of file ExtendedMath.h.

     {
         return val * val;
     }

template<typename T >

T mathsFunc::tan ( T x )

Todo:: should be properly computed

Definition at line 146 of file ExtendedMath.h.

References cos(), and sin().

Referenced by ChuteWithHopper::addHopper(), and ChuteWithHopper::setHopper().

                                     {
         return sin(x)/cos(x);
     }

Functions

Detailed Description

Function Documentation