/*
QtWagon: a project about 3D objects.
Science and technology promotion license applied. Third party license automatically cascaded.
Zhikai Wang/ www.heteroclinic.net 2013
You can do anything with this file or any file(s) published as part QtWagon project, given this header is kept.
*/
#ifndef __MATHVECTOR
#define __MATHVECTOR
#include <iostream>
#include "mathPoint.h"
// 1. ?! inheritance of assignment, X Y Z, return ? The answer by now (Aug 22) is don't do Point p = Vector();
// 2. ?! mathVector(const mathVector & nmp): mathPoint(nmp).
template <class T>
class mathVector: public mathPoint<T> {
public:
	mathVector(const T ndata[]):mathPoint(ndata) {
		//data[0] = ndata[0]; 	data[1] = ndata[1]; 	data[2] = ndata[2];
	}
	mathVector(T nx = (T)0.0, T ny = (T)0.0, T nz= (T)0.0): mathPoint(nx,ny,nz){
		data[0] = nx; 	data[1] = ny; 	data[2] = nz;
	}
	//// copy constructor
	// our assumption is copy/assign ONLY among same class
	mathVector(const mathVector & nmp): mathPoint(nmp) {
		data[0] = nmp.getx(); 	data[1] = nmp.gety(); 	data[2] = nmp.getz();
	}
	// assignment operator
	mathVector & operator=(const mathVector& nmv) {
		if (this == &nmv) return *this;  // time-saving self-test
		mathPoint::operator = (nmv);
		return *this;                    // for daisy-chaining
	}
	// equality test operator
	bool operator==(const mathVector& nmv) {
		return (mathPoint::operator == (nmv));
	}
	virtual ~ mathVector() {
		//std::cout<<" ~orientation() "<<std::endl;
	}
	//norm 
	//#define Length(v)  sqrt(v.x*v.x+v.y*v.y+v.z*v.z)	// Vector length
	//#define dotProduct(v,w) (v.x*w.x+v.y*w.y+v.z*w.z)
	T firstnorm () const {
		 //fabs(l-r)
		T d0 = fabs(data[0]);
		T d1 = fabs(data[1]);
		T d2 = fabs(data[2]);
		T d3 = ( d0 > d1 )? d0 : d1;
		return (d2>d3)? d2:d3;
	}
	T norm ( ) const{
		//T norm ( ) {
		return sqrt( data[0]* data[0] + data[1]* data[1] + data[2]* data[2]);
	}
	//const mathVector<T> & normalize ( ) { // can't call T n1 =  this->norm();//??
	T  normalize ( ) {
		//T n = norm();//?? Neither 'T norm ( ) const{' nor 'T norm ( ) {'
		T n = sqrt( data[0]* data[0] + data[1]* data[1] + data[2]* data[2]);
		//if (n > 0.0) {
		if (! isAlmost(n,(T)0.0)) {

			data[0] /= n;data[1] /= n;data[2] /= n;
			return (T)1.0; 
		}
		else 
			return (T)0.0;
		//return *this;
	}

	//operators vv: + - cross angleBetweenVW;  pp: -;  
	template <class U>
	friend const mathVector<U> operator + (const mathVector<U> & lhs , const mathVector<U> & rhs );
	template <class U>
	friend const mathVector<U> operator - (const mathVector<U> & lhs , const mathVector<U> & rhs );
	template <class U>
	friend const mathVector<U> operator * (const mathVector<U> & lhs , U l );


	template <class U>
	friend const mathVector<U> cross (const mathVector<U> & V , const mathVector<U> & W );
	template <class U>
	friend U angleBetweenVW (const mathVector<U> & V , const mathVector<U> & W, mathVector<U> & ref );


	template <class U>
	friend U dot (const mathVector<U> & V , const mathVector<U> & W );

};
template <class U>
const mathVector<U> operator * (const mathVector<U> & lhs , U l ) {
	return mathVector<U>(lhs.getx()*l,lhs.gety()*l,lhs.getz()*l);
};

template <class U>
const mathVector<U> operator + (const mathVector<U> & lhs , const mathVector<U> & rhs ) {
	return mathVector<U> (lhs.data[0] + rhs.data[0],lhs.data[1] + rhs.data[1],lhs.data[2] + rhs.data[2]);
};
template <class U>
const mathVector<U> operator - (const mathVector<U> & lhs , const mathVector<U> & rhs ){
	return mathVector<U> (lhs.data[0] - rhs.data[0],lhs.data[1] - rhs.data[1],lhs.data[2] - rhs.data[2]);
};
template <class U>
U dot (const mathVector<U> & lhs , const mathVector<U> & rhs ) {
	return (lhs.data[0] * rhs.data[0]+ lhs.data[1] * rhs.data[1]+lhs.data[2] * rhs.data[2]);
};
// nomalization and anomalies is not this func' responsibility
template <class U>
const mathVector<U> cross (const mathVector<U> & V , const mathVector<U> & W ) {
	//crossVector->x = (v->y*w->z - v->z*w->y) / crossVectorLength;
	//crossVector->y = (v->z*w->x - v->x*w->z) / crossVectorLength;
	//crossVector->z = (v->x*w->y - v->y*w->x) / crossVectorLength;
	return mathVector<U> (
		(V.data[1]*W.data[2] - V.data[2]*W.data[1]) ,
		(V.data[2]*W.data[0] - V.data[0]*W.data[2]) ,
		(V.data[0]*W.data[1] - V.data[1]*W.data[0]) 
		);
};
template <class U>
U angleBetweenVW2 (const mathVector<U> & V , const mathVector<U> & W, mathVector<U> & ref ) {
	// it is required that ref is orthorgonal to v1 and v2 plane, if v1 and v2 are not colinear
	U l1 = V.norm(), l2 = W.norm();
	if (isAlmost<U>(l1, (U)0.0) || isAlmost<U>(l2, (U)0.0) ) {
		std::cout<<"Error! In angleBetweenVW2, l1 or l2 is 0.0."<<std::endl;
		//exit(-1);
		//return -1.0;
	}
	U tmpr = dot(V,W)/(l1*l2);
	if (isAlmost<U>(tmpr,(U)1.0))
		tmpr = (U)1.0;
	if (isAlmost<U>(tmpr,(U)-1.0))
		tmpr = (U)-1.0;
	if (isAlmost<U>(tmpr,(U)0.0))
		tmpr = (U)0.0;
	U theta = (     ((U)acos( tmpr)) )* ((U)180.0) / ((U)PI);
	//if ( 

	mathVector<U> tmpv1 = cross(V,W);
	tmpv1.normalize();
	mathVector<U> tmpv2 = ref;
	tmpv2.normalize(); 

	U tmpr2 = (ref - tmpv1).norm();

	if ( isAlmost<U>(tmpr2,(U)0.0) )
		return theta;
	else if ( isAlmost<U>(tmpr2,(U)2.0) ) 
		return ((U)-1.0)*theta;
	else {
		std::cout<<"angleBetweenVW2 (const mathVector<U> & V , const mathVector<U> & W, mathVector<U> & ref )  Error! tmpr2 is neither 0.0 or 2.0 ."<<std::endl;
		//exit(-1);

		return (U)-2.0;
	}


}
// anomalies is not this func' responsibility
template <class U>
U angleBetweenVW (const mathVector<U> & V , const mathVector<U> & W, mathVector<U> & ref ) {
	ref = cross(V,W);
	// |ref| should be done after calling this func
	// The angle returned is always greater than 0, V-W-ref right hand rule
	// -1.0 means anomalies
	// return degree
	U l1 = V.norm(), l2 = W.norm();
	if (isAlmost<U>(l1, (U)0.0) || isAlmost<U>(l2, (U)0.0) ) 
		return -1.0;
	U tmpr = dot(V,W)/(l1*l2);
	if (isAlmost<U>(tmpr,(U)1.0))
		tmpr = (U)1.0;
	if (isAlmost<U>(tmpr,(U)-1.0))
		tmpr = (U)-1.0;

	if (isAlmost<U>(tmpr,(U)0.0))
		tmpr = (U)0.0;
	U theta = acos( tmpr)* ((U)180.0) / ((U)PI);
	return theta;
};
template <class U>
const mathVector<U> operator - (const mathPoint<U> & lhs , const mathPoint<U> & rhs ) {
	return mathVector<U> (lhs.getx() - rhs.getx(),lhs.gety() - rhs.gety(),lhs.getz() - rhs.getz());

};
template <class T>
void drawVector(const mathVector<T> &v,
	T scale = 1.0, mathPoint<T> position = mathPoint<T>((T)0.0,(T)0.0,(T)0.0),
	mathVector<T> color  = mathVector<T>((T)0.5,(T)0.5,(T)0.5)
	) {
		glPushAttrib(GL_ALL_ATTRIB_BITS);
		glPushMatrix();
		glDisable(GL_LIGHTING);

#ifdef __USE_DOUBLE
		glTranslated(position.getx(),position.gety(),position.getz());
		glColor3dv(color.getDataPointer());
#else
		glTranslatef(position.getx(),position.gety(),position.getz());
		glColor3fv(color.getDataPointer());
#endif

		
		glBegin(GL_LINES);
#ifdef __USE_DOUBLE
		glVertex3d((T)0.0,(T)0.0,(T)0.0);
		glVertex3d(v.getx()*scale,v.gety()*scale,v.getz()*scale);
#else
		glVertex3f((T)0.0,(T)0.0,(T)0.0);
		glVertex3f(v.getx()*scale,v.gety()*scale,v.getz()*scale);
#endif
		glEnd();

		glPopMatrix();
		//glEnable(GL_LIGHTING);
		glPopAttrib();

};
template <class T>
void drawVectorAxis(T scale = (T)1.0, mathPoint<T> position = mathPoint<T>((T)0.0,(T)0.0,(T)0.0),
	mathVector<T> colorx  = mathVector<T>((T)1.0,(T)0.0,(T)0.0),
	mathVector<T> colory  = mathVector<T>((T)0.0,(T)1.0,(T)0.0),
	mathVector<T> colorz  = mathVector<T>((T)0.0,(T)0.0,(T)1.0)
	) {
		glPushAttrib(GL_ALL_ATTRIB_BITS);
		glPushMatrix();
		glDisable(GL_LIGHTING);
#ifdef __USE_DOUBLE
		glTranslated(position.getx(),position.gety(),position.getz());
		glColor3dv(colorx.getDataPointer());
#else
		glTranslatef(position.getx(),position.gety(),position.getz());
		glColor3fv(colorx.getDataPointer());
#endif


		glBegin(GL_LINES);
#ifdef __USE_DOUBLE
		glVertex3d((T)0.0,(T)0.0,(T)0.0);
		glVertex3d(scale*((T)1.0),(T)0.0,(T)0.0);
#else
		glVertex3f((T)0.0,(T)0.0,(T)0.0);
		glVertex3f(scale*((T)1.0),(T)0.0,(T)0.0);
#endif

		glEnd();

#ifdef __USE_DOUBLE
		glColor3dv(colory.getDataPointer());
#else
		glColor3fv(colory.getDataPointer());
#endif

		glBegin(GL_LINES);
#ifdef __USE_DOUBLE
		glVertex3d((T)0.0,(T)0.0,(T)0.0);
		glVertex3d((T)0.0,scale*((T)1.0),(T)0.0);
#else
		glVertex3f((T)0.0,(T)0.0,(T)0.0);
		glVertex3f((T)0.0,scale*((T)1.0),(T)0.0);
#endif
		glEnd();

#ifdef __USE_DOUBLE
		glColor3dv(colorz.getDataPointer());;
#else
		glColor3fv(colorz.getDataPointer());
#endif
		
		glBegin(GL_LINES);
#ifdef __USE_DOUBLE
		glVertex3d((T)0.0,(T)0.0,(T)0.0);
		glVertex3d((T)0.0,(T)0.0,scale*((T)1.0));
#else
		glVertex3f((T)0.0,(T)0.0,(T)0.0);
		glVertex3f((T)0.0,(T)0.0,scale*((T)1.0));
#endif

		glEnd();
		glPopMatrix();
		glPopAttrib();
}

#endif