// TRFMath.h

#ifndef TRFMath_H
#define TRFMath_H

#include <cmath>
#include <cstdlib>
#include <limits>

namespace trf {

//**********************************************************************

// Mathematical constants -- values are in cpp file.

// pi
extern double PI;
extern double TWOPI;
extern double PI2;

//**********************************************************************

// Physical constants.

// speed of light in cm/sec
extern double CLIGHT;

// Factor connecting curvature to momentum:
// C = 1/Rc = BFAC * B * q/pT
// where B = magnetic field in Tesla, q = charge in natural units 
// (electron has -1) and pT = transverse momentum in GeV/c.
extern double BFAC;

//**********************************************************************

// Modulus functions.

// Similar to fmod but return value in positive range only.
double fmod1( double value, double range );

// Similar to fmod1 but returns range (-R/2,R/2) instead of (0,R).
inline double fmod2( double value, double range ) {
  return fmod1( value+0.5*range, range ) - 0.5*fabs(range);
}

//**********************************************************************

// Equality

// Function is_equal returns true if values are exactly the same
// or differ by a very small amount.  We set that small amount to
// be twice the machine precision.  Values differing by one tick
// will return true, those differing by three or more will return
// false.  The behavior for those differing by two depends on the value
// and the platform.

inline bool is_equal(double x1, double x2) {
  if ( x1 == 0.0 ) return x2 == 0.0;
  if ( x2 == 0.0 ) return false;
  double eps = std::numeric_limits<double>::epsilon();
  // 1 tick = 0.5, 2tick = 1, ...
  // Choose an intermediate value
  double maxdif = 0.7*eps;
  double num = x1 - x2;
  double den = fabs(x1) + fabs(x2);
  double rat = fabs(num/den);
  return rat < maxdif;
}

//**********************************************************************

// Zero check

// Function is_zero returns true if the value is within epsilon
// of zero.  By definition, 1.0+epsilon is the closest value to 1.0.
// Thus, this function is sensible for variables whose natural scale
// or range is 1.0.

inline bool is_zero(double x) {
  double eps = std::numeric_limits<double>::epsilon();
  return fabs(x) < eps;
}

//**********************************************************************

// asinrat = asin(x)/x

double asinrat(double x);

//**********************************************************************

}  // end namespace trf

#endif