// -*- C++ -*-

#ifndef __SCHOOL_THREEVECTOR_H__
#define __SCHOOL_THREEVECTOR_H__ 1


//   Standard includes
#include <cmath>
#include <iostream>


namespace school {
  


  class threevector
  {
    
  protected:
    //        data member
    double _M_x, _M_y, _M_z;

  public:
    //  constructors
    threevector(double x = 0.0, double y=0.0, double z=0.0) : _M_x(x), _M_y(y), _M_z(z) {}
    
    //  copy
    threevector(const threevector&) = default;
    threevector& operator=(const threevector&) = default;

    //  destructor
    ~threevector() = default;
    
    //  elements access
    const double& X() const { return _M_x;}
    const double& Y() const { return _M_y;}
    const double& Z() const { return _M_z;}
    
    double& X() { return _M_x;}
    double& Y() { return _M_y;}
    double& Z() { return _M_z;}
    
    //   computed assignments
    threevector& operator+=(const threevector& a) {
      _M_x += a._M_x; _M_y += a._M_y; _M_z += a._M_z;
      return *this;
    }
    
    threevector& operator-=(const threevector& a) {
      _M_x -= a._M_x; _M_y -= a._M_y; _M_z -= a._M_z;
      return *this;
    }
    
    threevector& operator*=(double a) {
      _M_x *= a; _M_y *= a; _M_z *= a;
      return *this;
    }
    
    threevector& operator/=(double a) {
      _M_x /= a; _M_y /= a; _M_z /= a;
      return *this ;
    }
    
    //        magnitude and the transvers component squared
    double mag2 () const { return _M_x*_M_x + _M_y*_M_y + _M_z*_M_z;}
    double perp2() const { return _M_x*_M_x + _M_y*_M_y;}

    //        magnitude and the transverse component
    double mag () const { return std::sqrt(this -> mag2());}
    double perp() const { return std::sqrt(this -> perp2());}
    
    //        azimuth and polar angles
    double phi() const {
      return _M_x == 0.0 && _M_y == 0.0 ? 0.0 : std::atan2(_M_y,_M_x);
    }
    
    double theta() const {
    double p = this -> perp();
      return p == 0.0 && _M_z == 0.0 ? 0.0 : std::atan2(p, _M_z);
    }
  };
  
  
  //  unary operations
  inline threevector operator+(const threevector& x) {
    return x;
  }
  
  inline threevector operator-(const threevector& x) {
    return threevector(-x.X(), -x.Y(), -x.Z());
  }
  
  //     other operators
  inline threevector operator+(const threevector& a, const threevector& b) {
    return threevector(a) += b;
  }
  
  inline threevector operator-(const threevector& a, const threevector& b) {
    return threevector(a) -= b;
  }
  
  inline threevector operator*(const threevector& a, double b) {
    return threevector(a) *= b;
  }
  
  inline threevector operator*(double b, const threevector& a) {
    return threevector(a) *= b;
  }
  
  inline threevector operator/(const threevector& a, double b) {
    return threevector(a) /= b;
  }
  
  inline double operator*(const threevector& a, const threevector& b) {
    return a.X()*b.X() + a.Y()*b.Y() + a.Z()*b.Z();
  }
  
  inline bool operator==(const threevector& a, const threevector& b) {
    return a.X() == b.X() && a.Y() == b.Y() && a.Z() == b.Z();
  }
  
  inline bool operator!=(const threevector& a, const threevector& b) {
    return a.X()!= b.X() || a.Y() != b.Y() || a.Z() != b.Z();
  }
  
  inline double dot(const threevector& a, const threevector& b) {
    return a.X()*b.X() + a.Y()*b.Y() + a.Z()*b.Z();
  }
  
  inline threevector cross(const threevector& a, const threevector& b)
  {
    return threevector(a.Y()*b.Z() - a.Z()*b.Y(),
                       a.Z()*b.X() - a.X()*b.Z(),
                       a.X()*b.Y() - a.Y()*b.X());
  }
  
  //        Specializations
  inline double cosAngle(const threevector& a, const threevector& b)
  {
    double ptot2 = a.mag2()*b.mag2();
    return ptot2 <= 0.0 ? 1.0 : dot(a, b)/std::sqrt(ptot2);
  }
  
  
  inline double angle(const threevector& a, const threevector& b)
  {
    double ptot2 = a.mag2()*b.mag2();
    return ptot2 <= 0.0 ? 0.0 : std::acos(dot(a, b)/std::sqrt(ptot2));
  }
  
  
  //    I/O operations
  inline std::ostream& operator<<(std::ostream& os, const threevector& q) {
	return os<<"("<<q.X()<<","<<q.Y()<<","<<q.Z()<<")";
  }
  
}

#endif