#if !defined(VECTOR3D)
	#define VECTOR3D
	#include "vector3d.h"
#endif

// ----------------------------------------------
// -------------- vector3d()  -------------------

vector3d::vector3d(void)
{
	x=y=z=0.0;

}

// ----------------------------------------------
// -------------- vector3d(ux,uy,uz)  -------------------

vector3d::vector3d(double ux, double uy, double uz)
{
	x=ux;
	y=uy;
	z=uz;
}

// ----------------------------------------------
// --------------      -u     -------------------

vector3d vector3d::operator-(void)
{
	vector3d tmp;
	
	tmp.x=-x;
	tmp.y=-y;
	tmp.z=-z;
	
	return tmp;
	
}

// ----------------------------------------------
// --------------       u=v     -----------------

vector3d vector3d::operator=(const vector3d&  u)
{
	vector3d tmp;

	if (this==&u)
		return *this;
	else
	{ 
		x=u.x;
		y=u.y;
		z=u.z;
		
		return *this;
	}
}

// ----------------------------------------------
// --------------       u+v     -----------------

vector3d vector3d::operator+(const vector3d& v)
{
	vector3d tmp;
	
	tmp.x=x+v.x;
	tmp.y=y+v.y;
	tmp.z=z+v.z;
	
	return tmp;
}

// ----------------------------------------------
// --------------       u-v     -----------------

vector3d vector3d::operator-(const vector3d& v)
{
	vector3d tmp;
	
	tmp.x=x-v.x;
	tmp.y=y-v.y;
	tmp.z=z-v.z;
	
	return tmp;
}

// ----------------------------------------------
// -----------     u*v (scalar)     -------------

double vector3d::operator*(const vector3d& v)
{
	double temp=(x*v.x+y*v.y+z*v.z);
//cerr << "scalar *" << temp << endl;
	return temp;
}

// ----------------------------------------------
// -----------          a*u         -------------

vector3d operator*(double a, const vector3d& u)
{
vector3d tmp;
	
	tmp.x=a*u.x;
	tmp.y=a*u.y;
	tmp.z=a*u.z;
	
	return tmp;

}

// ----------------------------------------------
// -----------         u*a          -------------

vector3d operator*(const vector3d& u, double a)
{
	vector3d tmp;
	
	tmp.x=a*u.x;
	tmp.y=a*u.y;
	tmp.z=a*u.z;
	
	return tmp;

}

// ----------------------------------------------
// -----------         u/a          -------------

vector3d operator/(const vector3d& u, double a)
{
	vector3d tmp;
	
	if(a!=0.0)
	{
		tmp.x=u.x/a;
		tmp.y=u.y/a;
		tmp.z=u.z/a;
	}
	else
		cerr << "vector3d::operator/, try to divide by zero.\n";
		
	return tmp;
}

// ----------------------------------------------
// -----------         u/a          -------------

double vector3d::norm(void)
{
	return sqrt(x*x+y*y+z*z);
}


// ----------------------------------------------
// -----------     u^v (cross)      -------------

vector3d vector3d::operator^(const vector3d& v)
{
	vector3d tmp;
	
	tmp.x=y*v.z-v.y*z;
	tmp.y=-x*v.z+v.x*z;
	tmp.z=x*v.y-v.x*y;
	
	return tmp;
}

// ----------------------------------------------
// -----------     u[i]      -------------

double& vector3d::operator[](const int& index)
{
	try
	{
		if (!((index>=1) && (index<=3)))
			throw "Index out of range from class vector3d::operator[]";
			
		switch(index)
		{
			case 1:
				return x;
				break;
			case 2:
				return y;
				break;
			case 3:
				return z;
				break;
		}	
		
	}
	
	catch(const char *message)
	{
		cerr << message << "\n";
		dummy=knan;
		return dummy;
	}
}


// ----------------------------------------------
// -----------       cout << u      -------------

ostream& operator<<(ostream& s, const vector3d &u)
{
//	s << "(" << u.x << "," << u.y << "," << u.z << ")";
	s << u.x << " " << u.y << " " << u.z ;

	return s;
}

// --------------------------------------------------
// ---------------       angle      -----------------

double angle(vector3d u, vector3d v)
{
	double prod=u*v/(u.norm()*v.norm());
	
	if(prod>1.0)
		prod=1.0;
	if(prod<-1.0)
		prod=-1.0;
	
	return acos(prod);
}