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.
Mdouble	chi_squared (Mdouble x, int k)
	This is a chi_squared function return the value x and degrees of freedom k.
Mdouble	chi_squared_prob (Mdouble x, int k)
	This is the function which actually gives the probability back using a chi squared test.

---

Detailed Description

Namespace for some extra maths function that are often needed.

---

Function Documentation

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

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

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	(	Mdouble	x_max,
		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

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::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.

References constants::sqrt_pi.

Referenced by chi_squared(), and ChuteWithHopper::create_inflow_particle().

{
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