// -*- C++ -*- #ifndef __SCHOOL_THREEVECTOR_H__ #define __SCHOOL_THREEVECTOR_H__ 1 // Standard includes #include #include 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<<"("<