#include "vector3.h"

#include "auxmath.h"
#include "matrix3.h"

const evector3 ux3(1.0,0.0,0.0);
const evector3 uy3(0.0,1.0,0.0);
const evector3 uz3(0.0,0.0,1.0);
const evector3 u03(0.0,0.0,0.0);

const evector3f ux3f(1.0,0.0,0.0);
const evector3f uy3f(0.0,1.0,0.0);
const evector3f uz3f(0.0,0.0,1.0);
const evector3f u03f(0.0,0.0,0.0);


std::ostream &operator<<(std::ostream &stream,const evector3 &vec)
{
  stream << "("<<vec.x<<","<<vec.y<<","<<vec.z<<")";
  return(stream);
}


evector3::evector3(): x(0.0), y(0.0), z(0.0)
{
}

evector3::evector3(const evector3& a): x(a.x),y(a.y),z(a.z)
{
}

//evector3::evector3(float fx, float fy, float fz): x(fx),y(fy),z(fz)
//{
//}

evector3::evector3(double fx, double fy, double fz): x(fx),y(fy),z(fz)
{
}

double evector3::len() const
{
  return(sqrt(x*x+y*y+z*z));
}

void evector3::normalize()
{
  double l = 1.0/len();
  x = x*l;
  y = y*l;
  z = z*l;   
}  

evector3 evector3::unit() const
{
  double l = 1.0/len();
  return(evector3(x*l,y*l,z*l));
}

evector3 evector3::proj(const evector3& a) const
{
  return((*this)*((*this)*a));
}

double evector3::operator *(const evector3& a) const
{
  return(a.x*x + a.y*y + a.z*z);
}

evector3 evector3::operator ^(const evector3& a) const
{
  return(evector3(a.y*z-a.z*y , a.z*x-a.x*z  , a.x*y-a.y*x));
//  return(evector3(y*a.z-z*a.y , z*a.x-x*a.z  , x*a.y-y*a.x));  // correct form of cross product
}

evector3 evector3::operator +(const evector3& a) const
{
  return(evector3(x+a.x,y+a.y,z+a.z));
}

evector3 evector3::operator -(const evector3& a) const
{
  return(evector3(x-a.x,y-a.y,z-a.z));
}

evector3 evector3::operator *(double a) const
{
  return(evector3(x*a,y*a,z*a));
}

evector3 evector3::operator /(double a) const
{
  double t = 1.0/a;
  return(evector3(x*t,y*t,z*t));
}

evector3& evector3::operator =(const evector3& a)
{
  x=a.x; y=a.y; z=a.z; return(*this);
}

evector3& evector3::operator +=(const evector3& a)
{
  x+=a.x;
  y+=a.y;
  z+=a.z;
  return(*this);
}

evector3& evector3::operator -=(const evector3& a)
{
  x-=a.x;
  y-=a.y;
  z-=a.z;
  return(*this);
}


evector3 evector3::eulerRot(double psi, double phi, double teta) const
{
  return(mEuler(psi,phi,teta)*(*this));
}

evector3 evector3::rot(double ang,const evector3& vec) const
{
  return(mRot(ang,vec)*(*this));
}

//const float DEGRAD = MPI/180.0;

evector3 evector3::rotx(double ang) const
{
  return(mRotX(ang)*(*this)); 
}  

evector3 evector3::roty(double ang) const
{
  return(mRotY(ang)*(*this)); 
}  

evector3 evector3::rotz(double ang) const
{
  return(mRotZ(ang)*(*this)); 
}  


evector3 polar3(const evector3& rot, double len)
{
  return(mRotY(rot.x)*mRotX(rot.y)*mRotZ(rot.z)*evector3(0.0,0.0,-len));
}  

evector3 polar3(double r, double teta, double phi)
{
  return(mRotZ(phi)*mRotX(teta)*evector3(0.0,0.0,-r)); 
}

evector3fp::evector3fp(): px(0x00) {}
evector3fp::evector3fp(float *vf): px(vf) {}


std::ostream &operator<<(std::ostream &stream,const evector3fp &vec)
{
  stream << "("<<vec.px[0]<<","<<vec.px[1]<<","<<vec.px[2]<<")";
  return(stream);
}

evector3f operator*(float v,const evector3fp &vec)
{
  return(vec*v);
}


evector3f::evector3f(): evector3fp(new float[3])
{
}

evector3f::evector3f(const evector3fp& a)
{
  px=new float[3];
  px[0]=a.px[0];
  px[1]=a.px[1];
  px[2]=a.px[2];
}

//evector3::evector3(float fx, float fy, float fz): x(fx),y(fy),z(fz)
//{
//}

evector3f::evector3f(float fx, float fy, float fz)
{
  px=new float[3];
  px[0]=fx; px[1]=fy; px[2]=fz;
}

evector3f::~evector3f()
{
  delete[] px;
}

float evector3fp::len() const
{
  return(sqrt(px[0]*px[0]+px[1]*px[1]+px[2]*px[2]));
}

void evector3fp::normalize()
{
  float l = 1.0/len();
  px[0]*=l;
  px[1]*=l;
  px[2]*=l;   
}  

evector3f evector3fp::unit() const
{
  float l = 1.0/len();
  return(evector3f(px[0]*l,px[1]*l,px[2]*l));
}

evector3f evector3fp::proj(const evector3fp& a) const
{
  return((*this)*((*this)*a));
}

float evector3fp::operator *(const evector3fp& a) const
{
  return(a.px[0]*px[0] + a.px[1]*px[1] + a.px[2]*px[2]);
}

evector3f evector3fp::operator ^(const evector3fp& a) const
{
  return(evector3f(a.px[1]*px[2]-a.px[2]*px[1] , a.px[2]*px[0]-a.px[0]*px[2]  , a.px[0]*px[1]-a.px[1]*px[0]));
//  return(evector3(y*a.z-z*a.y , z*a.x-x*a.z  , x*a.y-y*a.x));  // correct form of cross product
}

evector3f evector3fp::operator +(const evector3fp& a) const
{
  return(evector3f(px[0]+a.px[0],px[1]+a.px[1],px[2]+a.px[2]));
}

evector3f evector3fp::operator -(const evector3fp& a) const
{
  return(evector3f(px[0]-a.px[0],px[1]-a.px[1],px[2]-a.px[2]));
}

evector3f evector3fp::operator *(float a) const
{
  return(evector3f(px[0]*a,px[1]*a,px[2]*a));
}

evector3f evector3fp::operator /(float a) const
{
  float t = 1.0/a;
  return(evector3f(px[0]*t,px[1]*t,px[2]*t));
}

evector3fp& evector3fp::operator =(const evector3fp& a)
{
  px[0]=a.px[0]; px[1]=a.px[1]; px[2]=a.px[2]; return(*this);
}

evector3fp& evector3fp::operator +=(const evector3fp& a)
{
  px[0]+=a.px[0];
  px[1]+=a.px[1];
  px[2]+=a.px[2];
  return(*this);
}

evector3fp& evector3fp::operator -=(const evector3fp& a)
{
  px[0]-=a.px[0];
  px[1]-=a.px[1];
  px[2]-=a.px[2];
  return(*this);
}


evector3f evector3fp::eulerRot(float psi, float phi, float teta) const
{
  return(mEuler(psi,phi,teta)*(*this));
}

evector3f evector3fp::rot(float ang,const evector3fp& vec) const
{
  return(mRot(ang,vec)*(*this));
}

//const float DEGRAD = MPI/180.0;

evector3f evector3fp::rotx(float ang) const
{
  return(mRotX(ang)*(*this)); 
}  

evector3f evector3fp::roty(float ang) const
{
  return(mRotY(ang)*(*this)); 
}  

evector3f evector3fp::rotz(float ang) const
{
  return(mRotZ(ang)*(*this)); 
}  


evector3f polar3f(const evector3fp& rot, float len)
{
  return(mRotY(rot.px[0])*mRotX(rot.px[1])*mRotZ(rot.px[2])*evector3f(0.0,0.0,-len));
}  

evector3f polar3f(float r, float teta, float phi)
{
  return(mRotZ(phi)*mRotX(teta)*evector3f(0.0,0.0,-r)); 
}