#ifndef TRIPLETS_VEC_H_
#define TRIPLETS_VEC_H_

#include <math.h>
#include "util.h"
#include "layout-ptr.h"

using namespace layout_ptr;

#define F_PI 3.1415926f

/** 2D vector with layout support */
template <class Tl, class T> struct vec2_flex {
	// necessary for layouting magic to work
	typedef T layfield_t;
	typedef laygroup_t<Tl, layfield_t> group_t;
	// a provisional typedef, for struct-internal use only
	typedef vec2_flex<id_layout, T> vec2;
	// vector components
	union {T x; group_t dummy0;};
	union {T y; group_t dummy1;};

	// no need to construct for non-id layouts
	static hostdevice__ vec2 from(T x, T y) {
		vec2 a;
		a.x = x;
		a.y = y;
		return a;
	}
	static hostdevice__ vec2 from(T r) {	return from(r, r); }

	// but a need to cast into id layout from non-id layout
	hostdevice__ operator vec2() const { return vec2::from(x, y);	}

	// assignment operators; ordinary assignment must be overridden as well
	hostdevice__ vec2	operator= (const vec2 &a) {
		x = a.x;
		y = a.y;
		return (vec2)*this;
	}
	hostdevice__ vec2	operator+= (const vec2 &a) {
		x += a.x;
		y += a.y;
		return (vec2)*this;
	}
	hostdevice__ vec2	operator-= (const vec2 &a) {
		x -= a.x;
		y -= a.y;
		return (vec2)*this;
	}
	hostdevice__ vec2	operator*= (const vec2 &a) {
		x *= a.x;
		y *= a.y;
		return (vec2)*this;
	}
	hostdevice__ vec2	operator/= (const vec2 &a) {
		x /= a.x;
		y /= a.y;
		return (vec2)*this;
	}
};  // struct vec2_flex<T>

// specialize layflex_t for vec2<T> to provide additional info
template<class Tl, class T>
struct layflex_t<Tl, vec2_flex<id_layout, T> > {
	typedef vec2_flex<Tl, T> flex_t;
};

typedef vec2_flex<id_layout, float> fvec2;
// TODO: add alias for vec2<T> = vec2_flex<id_layout, T> when it's C++11
template <class Tl, class T> struct vec3_flex {
	typedef T layfield_t;
	typedef laygroup_t<Tl, layfield_t> group_t;
	typedef vec3_flex<id_layout, T> vec3;
	static hostdevice__ vec3 from(T x, T y, T z) {
		vec3 a;
		a.x = x;
		a.y = y;
		a.z = z;
		return a;
	}
	static hostdevice__ vec3 from(T r) { return from(r, r, r); }
	hostdevice__ operator vec3() const { return vec3::from(x, y, z); }
	hostdevice__ vec3 operator= (const vec3 &a) {
		x = a.x;
		y = a.y;
		z = a.z;
		return (vec3)*this;
	}
	hostdevice__ vec3 operator+= (const vec3 &a) {
		x += a.x;
		y += a.y;
		z += a.z;
		return (vec3)*this;
	}
	hostdevice__ vec3 operator-= (const vec3 &a) {
		x -= a.x;
		y -= a.y;
		z -= a.z;
		return (vec3)*this;
	}
	hostdevice__ vec3 operator*= (const vec3 &a) {
		x *= a.x;
		y *= a.y;
		z *= a.z;
		return (vec3)*this;
	}
	hostdevice__ vec3 operator/= (const vec3 &a) {
		x /= a.x;
		y /= a.y;
		z /= a.z;
		return (vec3)*this;
	}
	hostdevice__ vec2_flex<id_layout, T> xy() const { 
		return vec2_flex<id_layout, T>::from(x, y); 
	}
 	
	// vector components
	//T x, y, z;
	union { T x; group_t dummy0;};
	union { T y; group_t dummy1;};
	union { T z; group_t dummy2;};
};  // class vec3<T>

// specialize layflex_t for vec3<T> to provide additional info
template<class Tl, class T>
struct layflex_t<Tl, vec3_flex<id_layout, T> > {
	typedef vec3_flex<Tl, T> flex_t;
};

typedef vec3_flex<id_layout, float> fvec3;

// vec2 methods
template <class T> inline hostdevice__ 
vec2_flex<id_layout, T> operator+
(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(a.x + b.x, a.y + b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator-(const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(-b.x, -b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator-
(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(a.x - b.x, a.y - b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator*
(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(a.x * b.x, a.y * b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator*
(T a, const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(a * b.x, a * b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator*
(const vec2_flex<id_layout, T> &a, T b) {
	return vec2_flex<id_layout, T>::from(a.x * b, a.y * b);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator/
(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return vec2_flex<id_layout, T>::from(a.x / b.x, a.y / b.y);
}
template <class T> inline hostdevice__
vec2_flex<id_layout, T> operator/
(const vec2_flex<id_layout, T> &a, T b) {
	return vec2_flex<id_layout, T>::from(a.x / b, a.y / b);
}
template <class T> inline hostdevice__
T dot(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return a.x * b.x + a.y * b.y;
}
template <class T> inline hostdevice__
T cross(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return a.x * b.y - a.y * b.x;
}
template <class T> inline hostdevice__
T len2(const vec2_flex<id_layout, T> &a) { return dot(a, a); }
template <class T> inline hostdevice__
T dist2(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) 
{ return len2(a - b); }
template <class T> inline hostdevice__
T len(const vec2_flex<id_layout, T> &a) { return sqrtf(len2(a)); }
template <class T> inline hostdevice__
T dist(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) 
{ return len(a - b); }
template <class T> inline hostdevice__ 
T phi(const vec2_flex<id_layout, T> &a) {	return atan2f(a.y, a.x); }
template <class T> inline hostdevice__
T dphi(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	float d = phi(a) - phi(b);
	if(d < -F_PI)
		d += 2 * F_PI;
	else if(d >= F_PI)
		d -= 2 * F_PI;
	return d;
}  // dphi
// this function just returns a number whose sign is the same one as of dphi
template<class T> inline hostdevice__
T dphi_same_sign
(const vec2_flex<id_layout, T> &a, const vec2_flex<id_layout, T> &b) {
	return a.y * b.x - a.x * b.y;
}  // dphi_same_sign
// vector orthogonal to the current one and of same length
template<class T> inline hostdevice__
vec2_flex<id_layout, T> normal_to(const vec2_flex<id_layout, T> &a) { 
	return vec2_flex<id_layout, T>::from(a.y, -a.x); 
}
/** normalized vector */
template<class T> inline hostdevice__
vec2_flex<id_layout, T> normalized(const vec2_flex<id_layout, T> &a) {
	T invlen = 1 / len(a);
	return a * invlen;
}

// vec3 methods
template <class T> inline hostdevice__ 
vec3_flex<id_layout, T> operator+
(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(a.x + b.x, a.y + b.y, a.z + b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator-(const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(-b.x, -b.y, -b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator-
(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(a.x - b.x, a.y - b.y, a.z - b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator*
(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(a.x * b.x, a.y * b.y, a.z * b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator*(T a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(a * b.x, a * b.y, a * b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator*(const vec3_flex<id_layout, T> &a, T b) {
	return vec3_flex<id_layout, T>::from(a.x * b, a.y * b, a.z * b);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator/
(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from(a.x / b.x, a.y / b.y, a.z / b.z);
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> operator/(const vec3_flex<id_layout, T> &a, T b) {
	return vec3_flex<id_layout, T>::from(a.x / b, a.y / b, a.z / b);
}
template <class T> inline hostdevice__
T dot(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return a.x * b.x + a.y * b.y + a.z * b.z;
}
template <class T> inline hostdevice__
vec3_flex<id_layout, T> cross
(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	return vec3_flex<id_layout, T>::from
		(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
}
template <class T> inline hostdevice__
T len2(const vec3_flex<id_layout, T> &a) { return dot(a, a); }
template <class T> inline hostdevice__
T dist2(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) 
{ return len2(a - b); }
template <class T> inline hostdevice__
T len(const vec3_flex<id_layout, T> &a) { return sqrtf(len2(a)); }
template <class T> inline hostdevice__
T dist(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) 
{ return len(a - b); }
template <class T> inline hostdevice__ 
T phi(const vec3_flex<id_layout, T> &a) {	return atan2f(a.y, a.x); }
template <class T> inline hostdevice__
T dphi(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	float d = phi(a) - phi(b);
	if(d < -F_PI)
		d += 2 * F_PI;
	else if(d >= F_PI)
		d -= 2 * F_PI;
	return d;
}  // dphi
// this function returns the cosine of dphi between two vectors
template <class T> inline hostdevice__
T cosdphi(const vec3_flex<id_layout, T> &a, const vec3_flex<id_layout, T> &b) {
	vec2_flex<id_layout, T> a1 = a.xy(), b1 = b.xy();
	return dot(a1, b1) / sqrtf(len2(a1) * len2(b1));
}  // dcosphi
// this function just returns a number whose sign is the same one as of dphi
template<class T> inline hostdevice__
T dphi_same_sign(const vec3_flex<id_layout, T> &a, 
								 const vec3_flex<id_layout, T> &b) {
	return a.y * b.x - a.x * b.y;
}  // dphi_same_sign

#endif