#ifndef MATHEMATICALTOOLS_SEEN
#define MATHEMATICALTOOLS_SEEN
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>



using namespace std;

 


class mathTools
{


 public :
  // return the index K such that : 
  // vec[K] <= xx < vec[K+1]
  // vec is assumed to be increasingly ordered
  // if xx is outside , return -1 
  static int locateInVector(const vector<double>& vec, double xx)
  {
    int n= vec.size();
    int jl = 0;
    int ju = n+1;
    int jm;
    if ( vec[0] == xx ) return 0;
    if ( vec.back() <= xx || vec[0] > xx) return -1;
    while ( ju - jl > 1 ) 
      {
        jm = (ju + jl)/2;
	if ( vec[jm-1] <= xx ) jl = jm;
	else ju = jm;
      }
    return jl - 1; 
  }




  // algo. Boris (Birdsall p. 356) in 2D
  //  'tgArcMoitie' is equal to tg(theta/2), if theta is 
  // the desired rotation angle.
  // if tg(theta/2) is given exactly, the rotation is 'exact'
  // 
  static void borisBunemanRotation(double tgArcMoitie, double& vz, double& vx)
  {
    double t2 = tgArcMoitie * tgArcMoitie;
    double sinAngle = 2.0 * tgArcMoitie / ( 1 + t2 );
    double vauxz = vz + vx*tgArcMoitie;
    vx = vx - vauxz*sinAngle;
    vz = vauxz + vx*tgArcMoitie;
  }

};




class TRIDVECTOR 
{
  double vec_[3];

  public :

  TRIDVECTOR() 
    {
      int k;
      for (k=0; k<3; k++) vec_[k] = 0.0;
    }

TRIDVECTOR( const TRIDVECTOR& triv) 
{
  setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]); 
}



TRIDVECTOR(double x,double y,double z)
{
    vec_[0] = x;
    vec_[1] = y;
    vec_[2] = z;
}

 inline double& operator() (int index) { return vec_[index]; }
 inline const double& operator() (int index) const { return vec_[index]; }

inline TRIDVECTOR& operator= (const TRIDVECTOR& triv) 
{
  setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]);
  return *this;
}


inline TRIDVECTOR& operator = (double value) 
{
  setComponents( value,value, value); 
  return *this;
}

 inline void operator *= (const double& factor)
 {
   int k;
   for (k=0; k < 3; k++) vec_[k] *= factor;
   //   return *this;
 }



 inline TRIDVECTOR operator + (const TRIDVECTOR& v2)
 {
   //   int k;
   return TRIDVECTOR(vec_[0] + v2.vec_[0], vec_[1] + v2.vec_[1], vec_[2] + v2.vec_[2]);
     //   for (k=0; k < 3; k++) vec_[k] += v2.vec_[k];
     //   return *this;
 }

 inline void operator += (const TRIDVECTOR& v2)
 {
   int k;
   for (k=0; k < 3; k++) vec_[k] += v2.vec_[k];
   //   return *this;
 }


 inline TRIDVECTOR operator * (const double& facteur) const
 {
   return TRIDVECTOR( vec_[0] * facteur,  vec_[1] * facteur,  vec_[2] * facteur);
 }


 inline TRIDVECTOR operator * (const TRIDVECTOR& v2)
 {
   double auxx, auxy, auxz;
//   for ( k=0; k<3; k++) aux[k] = vec_[k];
   auxx  = vec_[1] * v2.vec_[2]  - vec_[2] * v2.vec_[1];
   auxy  = vec_[2] * v2.vec_[0]  - vec_[0] * v2.vec_[2];
   auxz  = vec_[0] * v2.vec_[1]  - vec_[1] * v2.vec_[0];
//   return *this;
return TRIDVECTOR(auxx, auxy, auxz);
 }



inline double getComponent(int index ) const {return vec_[index];}

 inline void setComponent(int index, double value )  {vec_[index] = value;}

 inline void incrementComponent(int index, double value )  {vec_[index] += value;}


 inline void getComponents(double& x,double& y,double& z) const
   {
     x = vec_[0];
     y = vec_[1];
     z = vec_[2];
   }

 inline void setComponents(double x, double y, double z) 
   {
     vec_[0] = x;
     vec_[1] = y;
     vec_[2] = z;
   }



  // algo. Boris (Birdsall p. 356)
  // the input vector 'tgArcMoitie' is in the direction of the axis of 
  // rotation. Its module is equal to tg(theta/2), if theta is 
  // the desired rotation angle.
  // if tg(theta/2) is given exactly, the rotation is 'exact'
  // 
  void borisBunemanRotation(const TRIDVECTOR& tgArcMoitie)
  {
    double t2 = tgArcMoitie.norm2();
    double fac = 2.0 / ( 1.0 + t2);
    TRIDVECTOR sinArc = tgArcMoitie * fac;
    TRIDVECTOR vecPrim(*this);
    vecPrim += (*this) * tgArcMoitie;
    (*this) += vecPrim * sinArc;
  }




// rotation in ZX plane
inline void rotateZXComponentsBuneman(double tgAngleMoitie )
{
  mathTools::borisBunemanRotation(tgAngleMoitie, vec_[2], vec_[0]);
}

// transform the vector Er, Etheta, Ez to a vector Ex,Ey,Ez, assuming the end
// of the vector at (x,y) in the plane (X,Y)
  inline void fromCylindricalToCartesian(double costet, double sintet)
  {
     double auxr = vec_[0];
     double auxtet = vec_[1];
     vec_[0] = auxr * costet - auxtet * sintet;
     vec_[1] = auxr * sintet + auxtet * costet;
  }


inline void clear()
{
     vec_[0] = 0.0;
     vec_[1] = 0.0;
     vec_[2] = 0.0;
}

inline double norm2() const
{
  return vec_[0]*vec_[0] + vec_[1]*vec_[1] + vec_[2]*vec_[2];
}
inline double norm() const
{
  return sqrt(abs(norm2()));
}

inline void renormalize()
{
  int k;
  double normeInv = 1.0/norm();
  for (k=0; k< 3 ; k++) vec_[k] *= normeInv;
}

inline void opposite()
{
  int k;
  for (k=0; k< 3 ; k++) vec_[k] = -vec_[k];
}

inline void print() const
{
  cout << " x comp. = " << vec_[0] << " y comp. = " << vec_[1] << " z comp. = " << vec_[2] << endl;
  cout << " norme = " << norm() << endl;
}

  string output_flow() const
  {
    ostringstream sortie;
    sortie << " "  << vec_[0] << " " <<  vec_[1] << "  " << vec_[2];
    return sortie.str();
  }

bool input_flow( ifstream& ifs)
 {
    bool test;
    if ( ifs >> vec_[0] >> vec_[1] >> vec_[2]) test= true;
    else test =  false;
    return test;
 }


};







#endif