#ifndef VECTOR_HPP_
#define VECTOR_HPP_

/** Single vector.
 * basically, just a x[3] but with helpful functions
 */
class vector {   
  public:
    double x[NDIM];

  vector();
  vector(double x1, double x2, double x3);
  ~vector();

  double Distance(const vector *y) const;
  double PeriodicDistance(const vector *y, const double *cell_size) const;
  double ScalarProduct(const vector *y) const;
  double Projection(const vector *y) const;
  double Norm() const ;
  double Angle(vector *y) const;

  void AddVector(const vector *y);
  void SubtractVector(const vector *y);
  void CopyVector(const vector *y);
  void RotateVector(const vector *y, const double alpha);
  void ProjectOntoPlane(const vector *y);
  void Zero(); 
  void One(double one);
  void Init(double x1, double x2, double x3);
  void Normalize();
  void Translate(const vector *x);
  void Mirror(const vector *x);
  void Scale(double **factor);
  void Scale(double *factor);
  void Scale(double factor);
  void MatrixMultiplication(double *M);
  void InverseMatrixMultiplication(double *M);
  void KeepPeriodic(ofstream *out, double *matrix);
  void LinearCombinationOfVectors(const vector *x1, const vector *x2, const vector *x3, double *factors);
  
  double CutsPlaneAt(vector *A, vector *B, vector *C);
  bool GetOneNormalVector(const vector *x1);
  bool MakeNormalVector(const vector *y1);
  bool MakeNormalVector(const vector *y1, const vector *y2);
  bool MakeNormalVector(const vector *x1, const vector *x2, const vector *x3);
  bool SolveSystem(vector *x1, vector *x2, vector *y, double alpha, double beta, double c);
  bool LSQdistance(vector **vectors, int dim);

  void AskPosition(double *cell_size, bool check);
  bool Output(ofstream *out) const;
};

ofstream& operator<<(ofstream& ost, vector& m);
vector& operator+=(vector& a, const vector& b);
vector& operator*=(vector& a, const double m);
vector& operator*(const vector& a, const double m);
vector& operator+(const vector& a, const vector& b);

#endif /*VECTOR_HPP_*/