/********************************************************************************
 *                                                                              *
 * ---------------------------  Vector Class  ----------------------------------*
 *                                                                              *
 *                                                   XianMing Chen              *
 *                                                   modified on Jan 7, 2002   *
 ********************************************************************************/

#ifndef _VECTOR_H_
#define _VECTOR_H_

#include <iostream.h>
#include <math.h>
#include <assert.h>

#define EPSILON 1.0E-6
#define PI 3.141592

class Vector {
 public:
  Vector(int d=3) : dims(d), data(0) { 
      if(dims) { data = new double[dims]; for(int i=0; i<dims; i++) data[i]=0; }
  }
  Vector(double x, double y, double z) : dims(3) { 
      data = new double[3]; data[0] = x; data[1] = y; data[2] = z; 
  }
  Vector(double *pv, int dimens=3) 
    {data=new double[dims=dimens]; for(int i=0;i<dims;i++) data[i]=pv[i]; }
  Vector(const Vector & other) : data(NULL), dims(0) { copy(other); }
  virtual  ~Vector() { delete [] data; }

  void copy(const Vector& org); 
  Vector& operator= (const Vector& rhs) { copy(rhs); return *this; }

  int  dimensions() const             { return dims; }
  int  dimensions(int dms);           // reset dims, copy old values.
  double& operator[] (int idx)             { assert(idx>=0&&idx<dims); return data[idx]; }
  double operator[] (int idx) const        { assert(idx>=0&&idx<dims); return data[idx]; }

  double norm() const {
      double ret = 0; for(int i=0; i<dims; i++) ret += data[i]*data[i]; return sqrt(ret); }
  double sqNorm() const {
      double ret = 0; for(int i=0; i<dims; i++) ret += data[i]*data[i]; return ret; }
  Vector& normalize() { double len = norm();
      for(int i=0; i<dims; i++) data[i] /= len; return *this; }  


  //;;;;;;;;;;;;;;;;;;;;;;;;;;; basic vector operations ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  Vector operator+ (const Vector &r) const { assert(dims==r.dims);
       Vector ret(dims); for(int i=0;i<dims;i++) ret[i] = data[i]+r.data[i]; return Vector(ret); }
  Vector& operator+= (const Vector &r) { assert(dims==r.dims);
       for(int i=0;i<dims;i++) data[i] += r.data[i]; return *this; }

  Vector operator- (const Vector &r) const { assert(dims==r.dims);
       Vector ret(dims); for(int i=0;i<dims;i++) ret[i] = data[i]-r.data[i]; return Vector(ret); }
  Vector& operator-= (const Vector &r) { assert(dims==r.dims);
       for(int i=0;i<dims;i++) data[i] -= r.data[i]; return *this; }

  //unary minus
  Vector operator- () const { Vector ret(dims); 
       for(int i=0;i<dims;i++) ret[i] = - data[i]; return Vector(ret); }

  //scaling
  Vector  operator* (const double & s) const { Vector ret(dims);
       for(int i=0;i<dims;i++) ret[i] = data[i] * s; return Vector(ret); }
  friend Vector operator* (const double& s, const Vector &v);

  Vector& operator*= (const double & s) {
       for(int i=0;i<dims;i++) data[i] *= s; return *this; }
  Vector  operator/ (const double & s) const { Vector ret(dims);
       for(int i=0;i<dims;i++) ret[i] = data[i] / s; return Vector(ret); }
  Vector& operator/= (const double & s) {
       for(int i=0;i<dims;i++) data[i] /= s; return *this; }


  //cross product of 3D vectors
  Vector operator* (const Vector &r) const { assert(dims==3 && dims==r.dims);
       Vector ret(dims);
       for(int i=0; i<dims; i++) 
	   ret[i] = data[(i+1)%3]*r.data[(i+2)%3] - r.data[(i+1)%3]*data[(i+2)%3];
       return Vector(ret); }
  Vector& operator*= (const Vector &r) { assert(dims==3 && dims==r.dims);
       for(int i=0; i<dims; i++) 
	   data[i] = data[(i+1)%3]*r.data[(i+2)%3] - r.data[(i+1)%3]*data[(i+2)%3];
       return *this; }

  //component multiplication of vectors
  Vector operator& (const Vector &r) const { assert(dims==r.dims);
       Vector ret(dims);
       for(int i=0; i<dims; i++) 
	   ret[i] = data[i]*r.data[i];
       return Vector(ret); }
  Vector& operator&= (const Vector &r) { assert(dims==r.dims);
       for(int i=0; i<dims; i++) 
	   data[i] = data[i]*r.data[i];
       return *this; }

  //dot product of 3D vectors
  double operator%  (const Vector &r) const { assert(dims==r.dims); double ret = 0;
       for(int i=0;i<dims;i++) ret += data[i] * r.data[i]; return ret; }

  //angle relative to another vector
  double angle(const Vector &r) const {
      return ( acos(*this%r/ (norm()*r.norm())) * 180/PI); }

  Vector projecTo(const Vector &r) const {
      return ( r/r.norm() * (*this % r) ); }

  friend ostream& operator<<   (ostream& os, const Vector& v);

 protected:
  double* data;
  int dims;
};  

#define RGBcolor Vector
#endif