mathsFunc Namespace Reference

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...
double	goldenSectionSearch (double(*function)(const double), double min, double cur, double max, double endCondition, double 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, double absError)
	Compares the difference of two Mdouble with an absolute error, useful in UnitTests. More...
bool	isEqual (Vec3D v1, Vec3D v2, double 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...

Detailed Description

Namespace for some extra maths function that are often needed.

Function Documentation

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 70 of file ExtendedMath.cc.

References 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 86 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;
    
}

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;
    }

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 47 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	(	double(*)(const double)	function,
		double	min,
		double	cur,
		double	max,
		double	endCondition,
		double	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 104 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;
    double resphi = 2 - 0.5 * (1 + std::sqrt(5));
    double x;
    if (max - cur > cur - min)
    {
        x = cur + resphi * (max - cur);
    }
    else
    {
        x = cur - resphi * (cur - min);
    }
    if(std::isnan(curVal))
        curVal = function(cur);
    double 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);
        }
    }
}

bool mathsFunc::isEqual	(	Mdouble	v1,
		Mdouble	v2,
		double	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 148 of file ExtendedMath.cc.

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

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

bool mathsFunc::isEqual	(	Vec3D	v1,
		Vec3D	v2,
		double	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 154 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);
}

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));
    }

template<typename T >

T mathsFunc::square

(

T

val

)

squares a number

Definition at line 91 of file ExtendedMath.h.

Referenced by DPMBase::areInContact(), DPMBase::checkParticleForInteraction(), helpers::computeCollisionTimeFromKAndDispAndEffectiveMass(), helpers::computeCollisionTimeFromKAndRestitutionCoefficientAndEffectiveMass(), helpers::computeDispFromKAndCollisionTimeAndEffectiveMass(), helpers::computeDispFromKAndRestitutionCoefficientAndEffectiveMass(), helpers::computeDisptFromCollisionTimeAndRestitutionCoefficientAndTangentialRestitutionCoefficientAndEffectiveMass(), SlidingFrictionInteraction::computeFrictionForce(), FrictionInteraction::computeFrictionForce(), helpers::computeKFromCollisionTimeAndDispAndEffectiveMass(), helpers::computeKFromCollisionTimeAndRestitutionCoefficientAndEffectiveMass(), helpers::computeKFromDispAndRestitutionCoefficientAndEffectiveMass(), ParticleSpecies::computeMass(), helpers::computeRestitutionCoefficientFromKAndCollisionTimeAndEffectiveMass(), helpers::computeRestitutionCoefficientFromKAndDispAndEffectiveMass(), LinearPlasticViscoelasticNormalSpecies::computeTimeStep(), LinearViscoelasticNormalSpecies::getCollisionTime(), StatisticsVector< T >::getCutoff2(), HertzianViscoelasticInteraction::getElasticEnergy(), LinearPlasticViscoelasticInteraction::getElasticEnergy(), LinearViscoelasticInteraction::getElasticEnergy(), IrreversibleAdhesiveInteraction::getElasticEnergy(), helpers::getMaximumVelocity(), StatisticsVector< T >::setCGWidth(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficient(), SlidingFrictionSpecies::setCollisionTimeAndNormalAndTangentialRestitutionCoefficientNoDispt(), LinearPlasticViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setCollisionTimeAndRestitutionCoefficient(), LinearViscoelasticNormalSpecies::setStiffnessAndRestitutionCoefficient(), Matrix3D::square(), and MatrixSymmetric3D::square().

    {
        return val * val;
    }