// TRFMath_t.cpp

// Test TRFMath.

#include "TRFMath.h"
#include <string>
#include <cassert>
#include <iostream>
#include <vector>
#include <cmath>

using std::string;
using std::cout;
using std::cerr;
using std::endl;
using std::vector;
using std::asin;
using std::abs;

using trf::PI;
using trf::TWOPI;
using trf::PI2;
using trf::CLIGHT;
using trf::BFAC;
using trf::fmod1;
using trf::fmod2;
using trf::asinrat;
using trf::is_equal;
using trf::is_zero;

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

int main( ) {

  string ok_prefix = "TRFMath (I): ";
  string error_prefix = "TRFMath test (E): ";

  cout << ok_prefix
       << "-------- Testing component TRFMath. --------" << endl;

  // Make sure assert is enabled.
  bool assert_flag = false;
  assert ( ( assert_flag = true, assert_flag ) );
  if ( ! assert_flag ) {
    cerr << "Assert is disabled" << endl;
    return 1;
  }

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

  cout << ok_prefix << "Test mathematical constants." << endl;
  cout << "   pi: " << PI << endl;
  cout << " 2*pi: " << TWOPI << endl;
  cout << " pi/2: " << PI2 << endl;
  assert( TWOPI = 2.0*PI );
  assert( PI2 = 0.5*PI );

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

  cout << ok_prefix << "Test physical constants." << endl;
  cout << "c (cm/sec): " << CLIGHT << endl;
  cout << "pT/(q*B*Rc): " << BFAC << endl;
  assert( CLIGHT > 0.0 );
  assert( BFAC > 0.0 );
  assert( fabs( (CLIGHT - 3.e10) / CLIGHT ) < 0.01 ); 
  assert( fabs( (BFAC - 0.003) / BFAC ) < 0.01 );

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

  cout << ok_prefix << "Test fmod1 and fmod2." << endl;
  double val[] =  { -2.7, -1.7, -0.7, 0.3,  1.3, 2.3, 200.3 };
  double val2[] = { -0.7,  0.3, -0.7, 0.3, -0.7, 0.3,   0.3 };
  double range1 = 1.0;
  double range2 = 2.0;
  double maxdiff = 1.e-10;
  for ( int i=0; i<7; ++i ) {
    double x1 = fmod1( val[i], range1 );
    double x2 = fmod2( val[i], range2 );
    cout << x1 << ' ' << x2 << endl;
    assert( fabs(x1-0.3) < maxdiff );
    assert( fabs(x2-val2[i]) < maxdiff );
    double x1m = fmod1( val[i], -range1 );
    double x2m = fmod2( val[i], -range2 );
    cout << x1m << ' ' << x2m << endl;
    assert( fabs(x1m-x1) < maxdiff );
    assert( fabs(x2m-x2) < maxdiff );
    cout << "----------------" << endl;
  }

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

  cout << ok_prefix << "Test is_equal for doubles." << endl;
  double eps = std::numeric_limits<double>::epsilon();
  cout << "eps = " << eps << endl;
  vector<double> vals;
  vals.push_back(1.0);
  vals.push_back(20.0);
  vals.push_back(3.e20);
  vals.push_back(4.e-20);
  vals.push_back(eps);
  vector<double>::const_iterator ival;
  for ( ival=vals.begin(); ival!=vals.end(); ++ival ) {
    double x0 = *ival;
    assert( is_equal(x0,x0) );
    cout << "Testing value " << x0 << endl;
    double dif = eps*x0;
    double x1 = x0 + dif;
    assert( x1 != x0 );
    double x2 = x1 + dif;
    assert( x2 != x1 );
    double x3 = x2 + dif;
    assert( x3 != x2 );
    cout << ok_prefix << "1*eps" << endl;
    cout << x1-x0 << " " << is_equal(x0,x1) << endl;
    assert( is_equal(x0,x1) );
    cout << ok_prefix << "2*eps" << endl;
    cout << x2-x0 << " " << is_equal(x0,x2) << endl;
    //assert( ! is_equal(x0,x2) );
    cout << ok_prefix << "3*eps" << endl;
    cout << x3-x0 << " " << is_equal(x0,x3) << endl;
    assert( ! is_equal(x0,x3) );
  }

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

  cout << ok_prefix << "Check is_equal for zero" << endl;
  double x0 = 0.0;
  assert( is_equal(x0,x0) );
  assert( ! is_equal(x0,eps) );
  assert( ! is_equal(eps,x0) );

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

  cout << ok_prefix << "Test is_zero." << endl;
  double zero = 0.0;
  assert( is_zero(zero) );
  assert( is_zero(0.5*eps) );
  assert( ! is_zero(eps) );
  assert( ! is_zero(1.5*eps) );

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

  cout << ok_prefix << "Test asinrat." << endl;
  { 
    assert( asinrat(0.0) == 1.0 );
    vector<double> vals;
    vals.push_back(1.e-15);
    vals.push_back(0.00004);
    vals.push_back(0.003);
    vals.push_back(1.0);
    vals.push_back(-.137);
    double close = 1.e-12;
    for ( vector<double>::const_iterator ival=vals.begin();
          ival!=vals.end(); ++ival ) {
       double x =*ival;
       double asx1 = asin(x);
       double asx2 = asinrat(x)*x;
       double dif = abs(asx2-asx1);
       cout << x << " " << asx1 << " " << asx2 << " " << dif << endl;
       assert( dif < close );
    }
  }

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

  cout << ok_prefix
       << "------------- All tests passed. -------------" << endl;
  return 0;

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

}