/////////////////////////////////////////////////////////////
//                                                         //
// Copyright (c) 2007-2014 by The University of Queensland //
// Centre for Geoscience Computing                         //
// http://earth.uq.edu.au/centre-geoscience-computing      //
//                                                         //
// Primary Business: Brisbane, Queensland, Australia       //
// Licensed under the Open Software License version 3.0    //
// http://www.opensource.org/licenses/osl-3.0.php          //
//                                                         //
/////////////////////////////////////////////////////////////


#ifndef __VECTOR3_HH
#define __VECTOR3_HH

//#include "Foundation/Matrix3.h"

//the error...
/*
VECTOR3_INLINE VecErr::VecErr(const string& m):MError(m)
{
  message.insert(0,"Vector3 "); 
}
*/

// constructors
VECTOR3_INLINE Vector3::Vector3()
{
  data[0]=0;
  data[1]=0;
  data[2]=0;
}

VECTOR3_INLINE Vector3::Vector3(double s)
{
  data[0]=s;
  data[1]=s;
  data[2]=s;
}

VECTOR3_INLINE Vector3::Vector3(double a,double b,double c)
{
  data[0]=a;
  data[1]=b;
  data[2]=c;
}

VECTOR3_INLINE Vector3::Vector3(const Vector3& rhs)
{
  data[0]=rhs.data[0];
  data[1]=rhs.data[1];
  data[2]=rhs.data[2];
}

// operators

VECTOR3_INLINE Vector3& Vector3::operator=(const Vector3& rhs)
{
  data[0]=rhs.data[0];
  data[1]=rhs.data[1];
  data[2]=rhs.data[2];
  return *this;
}

VECTOR3_INLINE Vector3& Vector3::operator=(double s)
{
  data[0]=s;
  data[1]=s;
  data[2]=s;
  return *this;
}

VECTOR3_INLINE Vector3& Vector3::operator-=(const Vector3& rhs)
{
  data[0]-=rhs.data[0];
  data[1]-=rhs.data[1];
  data[2]-=rhs.data[2];
  return *this;
}

VECTOR3_INLINE Vector3& Vector3::operator+=(const Vector3& rhs)
{
  data[0]+=rhs.data[0];
  data[1]+=rhs.data[1];
  data[2]+=rhs.data[2];
  return *this;
}

VECTOR3_INLINE Vector3 Vector3::operator+(const Vector3& rhs) const
{
  return Vector3(data[0]+rhs.data[0], data[1]+rhs.data[1], data[2]+rhs.data[2]);
}

VECTOR3_INLINE Vector3 Vector3::operator-(const Vector3& rhs) const
{
  return Vector3(data[0]-rhs.data[0], data[1]-rhs.data[1], data[2]-rhs.data[2]);
}

VECTOR3_INLINE Vector3 Vector3::operator-() const
{
  return Vector3( -data[0],-data[1],-data[2] );
}

/*
VECTOR3_INLINE Vector3 Vector3::operator*(const Matrix3 &m) const
{
  const double x = m(0,0)*data[0] + m(1,0)*data[1] + m(2,0)*data[2];
  const double y = m(0,1)*data[0] + m(1,1)*data[1] + m(2,1)*data[2];
  const double z = m(0,2)*data[0] + m(1,2)*data[1] + m(2,2)*data[2];

  return Vector3(x,y,z);
}
*/

VECTOR3_INLINE double Vector3::operator*(const Vector3& rhs) const
{
  return data[0]*rhs.data[0]+data[1]*rhs.data[1]+data[2]*rhs.data[2];
}

VECTOR3_INLINE Vector3 Vector3::operator*(double s) const 
{ 
   return Vector3(data[0]*s,data[1]*s,data[2]*s) ; 
} 

VECTOR3_INLINE Vector3& Vector3::operator*=(double rhs)
{
  data[0]*=rhs;
  data[1]*=rhs;
  data[2]*=rhs;
  return *this;
}

VECTOR3_INLINE Vector3& Vector3::operator/=(double c)
{
  data[0] /= c;
  data[1] /= c;
  data[2] /= c;

  return *this;
}

VECTOR3_INLINE Vector3 Vector3::operator/(double s) const 
{ 
   return Vector3(data[0]/s,data[1]/s,data[2]/s) ; 
}

VECTOR3_INLINE Vector3 Vector3::operator+(double s) const 
{
   return Vector3(data[0]+s, data[1]+s, data[2]+s) ; 
}

VECTOR3_INLINE Vector3 Vector3::operator-(double s) const 
{
   return Vector3(data[0]-s, data[1]-s, data[2]-s); 
}

VECTOR3_INLINE Vector3 Vector3::rotate(const Vector3 &axis, const Vector3 &axisPt) const
{
  const double phi = axis.norm();
  if (phi > 0.0)
  {
    const Vector3 r = *this - axisPt;
    const Vector3 n = axis/phi;
    const double cosPhi = cos(phi);
    const Vector3 rotatedR =
      r*cosPhi + n*((dot(n, r))*(1-cosPhi)) + cross(r, n)*sin(phi);
    return rotatedR + axisPt;
  }
  return *this;
}

VECTOR3_INLINE Vector3 &Vector3::operator+=(double s)
{
   data[0] += s;
   data[1] += s;
   data[2] += s;
   return *this; 
}

VECTOR3_INLINE Vector3 &Vector3::operator-=(double s)
{
   data[0] -= s;
   data[1] -= s;
   data[2] -= s;
   return *this; 
}

// vector product
// 9 Flops ( 6 mult, 3 sub ) 
VECTOR3_INLINE Vector3 cross(const Vector3& lhs,const Vector3& rhs)
{
  return Vector3(lhs.data[1]*rhs.data[2]-lhs.data[2]*rhs.data[1],
	      lhs.data[2]*rhs.data[0]-lhs.data[0]*rhs.data[2],
	      lhs.data[0]*rhs.data[1]-lhs.data[1]*rhs.data[0]);
}

//  dot product  

VECTOR3_INLINE double dot(const Vector3& v1, const Vector3& v2)
{
  return v1.data[0] * v2.data[0] + 
         v1.data[1] * v2.data[1] + 
         v1.data[2] * v2.data[2];
}

VECTOR3_INLINE Vector3 operator*(double f,const Vector3& rhs)
{
  return Vector3(f*rhs.data[0], f*rhs.data[1], f*rhs.data[2]);
}


// euclidian norm
// 6 Flops ( 3 mult, 2 add, 1 sqrt )
VECTOR3_INLINE double Vector3::norm() const
{
  return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]);
}

// weighted euclidian norm
VECTOR3_INLINE double Vector3::wnorm(double wx, double wy, double wz) const
{
  double dx=data[0]/wx;
  double dy=data[1]/wy;
  double dz=data[2]/wz;
  
  return sqrt(dx*dx+dy*dy+dz*dz);
}

// square of weighted euclidian norm
VECTOR3_INLINE double Vector3::wnorm2(double wx, double wy, double wz) const
{
  double dx=data[0]/wx;
  double dy=data[1]/wy;
  double dz=data[2]/wz;
  
  return dx*dx+dy*dy+dz*dz;
}

// square of the euclidian norm
// 5 Flops ( 3 mult, 2 add)
VECTOR3_INLINE double Vector3::norm2() const
{
  return data[0]*data[0]+data[1]*data[1]+data[2]*data[2];
}

// returns unit vector in direction of the original vector
// 9 Flops ( 3 mult, 2 add, 3 div, 1 sqrt ) 
VECTOR3_INLINE Vector3 Vector3::unit() const
{
  return (*this)/norm();
}

// per element min/max
VECTOR3_INLINE Vector3 cmax(const Vector3& v1,const Vector3& v2)
{
  Vector3 res;
  res.data[0]=v1.data[0]>v2.data[0] ? v1.data[0] : v2.data[0];
  res.data[1]=v1.data[1]>v2.data[1] ? v1.data[1] : v2.data[1];
  res.data[2]=v1.data[2]>v2.data[2] ? v1.data[2] : v2.data[2];
  return res;
}

VECTOR3_INLINE Vector3 cmin(const Vector3& v1,const Vector3& v2)
{
  Vector3 res;
  res.data[0]=v1.data[0]<v2.data[0] ? v1.data[0] : v2.data[0];
  res.data[1]=v1.data[1]<v2.data[1] ? v1.data[1] : v2.data[1];
  res.data[2]=v1.data[2]<v2.data[2] ? v1.data[2] : v2.data[2];
  return res;
}

// save version, throws exception if norm()==0
/*
VECTOR3_INLINE Vector3 Vector3::unit_s() const
{
  double n=norm();
  if(n==0) throw VecErr("norm() of data[2]ero-vector"); 
  Vector3 res(data[0],data[1],data[2]);
  return res/n;
}
*/

VECTOR3_INLINE double Vector3::max() const
{
  double m = ( data[0]>data[1] ? data[0] : data[1] );

  return ( m>data[2] ? m : data[2] );

}

VECTOR3_INLINE double Vector3::min() const
{
  double m = ( data[0]<data[1] ? data[0] : data[1] );

  return ( m<data[2] ? m : data[2] );
}

VECTOR3_INLINE bool Vector3::operator==(const Vector3& V) const
{
  return((data[0]==V.data[0])&&(data[1]==V.data[1])&&(data[2]==V.data[2]));
}

VECTOR3_INLINE bool Vector3::operator!=(const Vector3& V) const
{
  return((data[0]!=V.data[0])||(data[1]!=V.data[1])||(data[2]!=V.data[2]));
}

// per component min/max

VECTOR3_INLINE Vector3 comp_max(const Vector3& V1,const Vector3& V2)
{
  double x=(V1.x() > V2.x()) ? V1.x() : V2.x();
  double y=(V1.y() > V2.y()) ? V1.y() : V2.y();
  double z=(V1.z() > V2.z()) ? V1.z() : V2.z();
 
  return Vector3(x,y,z);
}

VECTOR3_INLINE Vector3 comp_min(const Vector3& V1,const Vector3& V2)
{
  double x=(V1.x() < V2.x()) ? V1.x() : V2.x();
  double y=(V1.y() < V2.y()) ? V1.y() : V2.y();
  double z=(V1.z() < V2.z()) ? V1.z() : V2.z();
 
  return Vector3(x,y,z);
}

// in/output

VECTOR3_INLINE std::ostream& operator << (std::ostream& ostr,const Vector3& V)
{
  const char delimiter = ' ';
  ostr
    << V.data[0] << delimiter
    << V.data[1] << delimiter
    << V.data[2];

  return ostr;
}

VECTOR3_INLINE std::istream& operator >> (std::istream& istr,Vector3& V)
{
  istr
    >> V.data[0]
    >> V.data[1]
    >> V.data[2];

  return istr;
}

#endif // __VECTOR3_HH