MercuryDPM  Beta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 >
square (T val)
 squares a number More...
 
template<typename 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().

71 {
72 
73  Mdouble prefactor = pow(2, k / 2.0) * gamma(k / 2.0);
74  Mdouble mainfactor = pow(x, k / 2.0 - 1) * exp(x / -2.0);
75 
76  return mainfactor / prefactor;
77 
78 }
double Mdouble
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:47
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().

87 {
88 
89 //The current value was picked by tried were it stopped effect the 4 d.p.
90  const int num_steps_per_unit = 100;
91  Mdouble sum = 0;
92  Mdouble x = 0;
93  long int num_steps = static_cast<int>(num_steps_per_unit * x_max);
94 //Use trapezional rule, but ignoring the ends
95  for (int i = 0; i < num_steps; i++)
96  {
97  x = x_max / num_steps * (i + 0.5);
98  sum = sum + chi_squared(x, k);
99  }
100  return 1.0 - sum * x_max / num_steps;
101 
102 }
double Mdouble
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.
Definition: ExtendedMath.cc:70
template<typename T >
T mathsFunc::cubic ( val)

calculates the cube of a number

Definition at line 99 of file ExtendedMath.h.

Referenced by ChuteBottom::setupInitialConditions().

100  {
101  return val * val * val;
102  }
template<typename T >
constexpr T mathsFunc::factorial ( const T  t)

factorial function

Definition at line 124 of file ExtendedMath.h.

Referenced by HGridOptimiser::histNumberParticlesPerCell().

125  {
126  return (t == 0) ? 1 : t * factorial(t - 1);
127  }
constexpr T factorial(const T t)
factorial function
Definition: ExtendedMath.h:124
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().

48 {
49  const Mdouble ep = 1e-5;
50 
51  if (gamma_in > 1.0 + ep)
52  {
53  return ((gamma_in - 1) * gamma(gamma_in - 1));
54  }
55  else
56  {
57 
58  if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
59  return 1.0;
60  else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
61  return constants::sqrt_pi;
62  else
63  return std::numeric_limits<Mdouble>::quiet_NaN();
64  }
65 } //end func gamma
const Mdouble sqrt_pi
Definition: ExtendedMath.h:43
double Mdouble
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
Definition: ExtendedMath.cc:47
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]functionA function pointer to the function of which you want to calculate the location of its minimum.
[in]minThe minimum location
[in]curThe current location
[in]maxThe maximum location
[in]endConditionThe algorithm terminates when abs(max - min) < endCondition
[in]curValThe value of the function at the current location (on default this value is calculated internally)

Definition at line 104 of file ExtendedMath.cc.

105 {
106  if (std::abs(max - min) < endCondition)
107  {
108  return 0.5 * (min + max);
109  }
110  std::cout << "Min=" << min << " Max=" << max << " diff=" << max - min << std::endl;
111  double resphi = 2 - 0.5 * (1 + std::sqrt(5));
112  double x;
113  if (max - cur > cur - min)
114  {
115  x = cur + resphi * (max - cur);
116  }
117  else
118  {
119  x = cur - resphi * (cur - min);
120  }
121  if(std::isnan(curVal))
122  curVal = function(cur);
123  double xVal = function(x);
124  if (xVal < curVal)
125  {
126  if (max - cur > cur - min)
127  {
128  return goldenSectionSearch(function, cur, x, max, endCondition, xVal);
129  }
130  else
131  {
132  return goldenSectionSearch(function, min, x, cur, endCondition, xVal);
133  }
134  }
135  else
136  {
137  if (max - cur > cur - min)
138  {
139  return goldenSectionSearch(function, min, cur, x, endCondition, curVal);
140  }
141  else
142  {
143  return goldenSectionSearch(function, x, cur, max, endCondition, curVal);
144  }
145  }
146 }
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...
bool mathsFunc::isEqual ( Mdouble  v1,
Mdouble  v2,
double  absError 
)

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

Parameters
[in]v1The first Mdouble
[in]v2The second Mdouble
[in]absErrorThe 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().

149 {
150  return std::abs(v1 - v2) <= absError;
151 }
bool mathsFunc::isEqual ( Vec3D  v1,
Vec3D  v2,
double  absError 
)

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

Parameters
[in]v1The first Vec3D
[in]v2The second Vec3D
[in]absErrorThe 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.

155 {
156  return isEqual(v1.X,v2.X,absError)&&isEqual(v1.Y,v2.Y,absError)&&isEqual(v1.Z,v2.Z,absError);
157 }
Mdouble X
the vector components
Definition: Vector.h:52
bool isEqual(Mdouble v1, Mdouble v2, double absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Mdouble Y
Definition: Vector.h:52
Mdouble Z
Definition: Vector.h:52
template<typename T >
int mathsFunc::sign ( 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().

84  {
85  return (T(0) < val) - (val < T(0));
86  }
template<typename T >
T mathsFunc::square ( 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().

92  {
93  return val * val;
94  }