#ifndef _MATRIX3_
#define _MATRIX3_

#include <math.h>
#include <float.h>

static const double m_pi=3.1415926535897932,twopi=2*m_pi,halfpi=m_pi/2;

struct Vector3 {
  float value[3];
  Vector3() {}
  Vector3(float x, float y, float z) {value[0]=x; value[1]=y; value[2]=z;}
  Vector3(float, float);
  void set_unit(float, float);
  bool operator==(const Vector3& a) const;
  bool operator!=(const Vector3& a) const {return !((*this)==a);}
  float operator*(const Vector3& a) const {return operator*(a.value);}
  float operator*(const float x[3]) const;
  void cross(const Vector3&, const Vector3&);
  float stp(const Vector3&, const Vector3&) const;
  void operator*=(float a) {value[0]*=a; value[1]*=a; value[2]*=a;}
  void operator/=(float a) {(*this)*=1/a;}
  void operator-=(const Vector3&);
  void operator+=(const Vector3&);
  void orthog(const Vector3&);
  float normalize() {float a=(float)sqrt((*this)*(*this)); (*this)/=a; return a;}
  float theta() const
    {return (float)atan2(value[2],sqrt(value[0]*value[0]+value[1]*value[1]));}
  float cotheta() const
    {return (float)atan2(sqrt(value[0]*value[0]+value[1]*value[1]),value[2]);}
  float minus_z() const {return -value[2];}
  float phi() const {return (float)atan2(value[1],value[0]);}
  float phi_positive() const {float p=phi(); return p>=0 ? p : p+(float)twopi;}
};

#endif // ndef _MATRIX3_