vector.h

030_Particle_system using opengl

#ifndef _VECTOR_H
#define _VECTOR_H
#include <math.h>

#ifndef PI180
#define PI180	0.0174532925199432957692369076848861f	// pi / 180
#endif

#ifndef PI
#define PI		3.1415926535897932384626433832795
#endif

#define MIN(x,y)		((x) < (y) ?  (x)      : (y))					// minimum
#define MAX(x,y)		((x) > (y) ?  (x)      : (y))					// maximum
#define ABS(x)			((x) < (0) ? -(x)      : (x))	

class vec2;
class vec4;

class vec  
{
public:
	union
	{
		struct
		{
			float x, y, z;
		};
		float v[3];
	};

	vec(){}
	inline vec( const vec &a);
	inline vec( const vec2 &a);
	inline vec( const vec4 &a);
	inline vec( const float a, const float b, const float c);
	inline void set( const float a, const float b, const float c);
	inline void clear();

	inline vec operator+() const;
	inline vec operator-() const;		// negative

	inline void operator+= ( const vec &a);
	inline void operator-= ( const vec &a);
	inline void operator*= ( const vec &a);
	inline void operator/= ( const vec &a);

	inline void operator+= ( const float f);
	inline void operator-= ( const float f);
	inline void operator*= ( const float f);
	inline void operator/= ( const float f);

	inline vec CROSS( const vec &a) const;
	inline float DOT3( const vec &a) const;
	inline float Length() const;
	inline float Length2() const;		// = Length ^ 2

	inline void RotX( const float angle);
	inline void RotY( const float angle);
	inline void RotZ( const float angle);
	
	inline void Normalize();
	inline void ArbitraryRotate(float theta, vec r);
/*
	inline friend vec operator+ ( const vec &a, const vec &b );
	inline friend vec operator- ( const vec &a, const vec &b );
	inline friend vec operator* ( const vec &a, const vec &b );
	inline friend vec operator/ ( const vec &a, const vec &b );

	inline friend vec operator+ ( const float f, const vec &a );
	inline friend vec operator- ( const float f, const vec &a );
	inline friend vec operator* ( const float f, const vec &a );
//	inline friend vec operator/ ( const float f, const vec &a );
	
	inline friend vec operator+ ( const vec &a, const float f );
	inline friend vec operator- ( const vec &a, const float f );
	inline friend vec operator* ( const vec &a, const float f );
	inline friend vec operator/ ( const vec &a, const float f );

	inline friend int operator< ( const vec &a, const vec &b );
	inline friend int operator<=( const vec &a, const vec &b );
	inline friend int operator> ( const vec &a, const vec &b );
	inline friend int operator>=( const vec &a, const vec &b );
	inline friend int operator==( const vec &a, const vec &b );
	inline friend int operator!=( const vec &a, const vec &b );

	inline friend vec CROSS ( const vec &a, const vec &b );
	inline friend float DOT3  ( const vec &a, const vec &b );
	inline friend vec Normalize( const vec &a );
	inline friend float Length( const vec &a );
	inline friend float Length2( const vec &a );
	inline friend float Distance ( const vec &a, const vec &b);
	inline friend float Distance2 ( const vec &a, const vec &b);
	inline friend vec MakeNormal( const vec &a, const vec &b, const vec &c);
	inline friend vec ReflectedRay( const vec &in, const vec &n);
	inline float PlaneDistance( const vec &normal, const vec &point);
	inline float PlanePointDelta( const vec &normal, float distance, const vec &point);
	inline float PlanePointDistance( const vec &normal, float distance, const vec &point);
	// Collision
	inline vec ClosestPointOnLine( const vec &a, const vec &b, const vec &point);
	inline int SpherePlaneCollision( const vec &center, float radius, const vec &normal, float distance, float *distanceCenterToPlane );
	double AngleBetweenVectors( const vec &a, const vec &b);
	inline int PointInsidePolygon( const vec &point, const vec &a, const vec &b, const vec &c);
	inline int PointInsidePolygon( const vec &point, const vec v[], unsigned int count);
	inline int EdgeSphereCollision( const vec &center, float radius, const vec &a, const vec &b, const vec &c);
	inline int EdgeSphereCollision( const vec &center, float radius, vec v[], unsigned int count);
	inline int SpherePolygonCollision(const vec &center, float radius, vec v[], unsigned int count, const vec &normal);
	inline int SpherePolygonCollision(const vec &center, float radius, const vec &a, const vec &b, const vec &c, const vec &normal);
	inline int LinePlaneCollision( const vec &a, const vec &b, const vec &normal, float distance);
	inline vec LinePlaneIntersection( const vec &a, const vec &b, const vec &normal, float distance);
	inline vec LinePlaneIntersectionDir( const vec &a, const vec &ab, const vec &normal, float distance);
	inline float LinePlaneIntersectionDirParameter( const vec &a, const vec &ab, const vec &normal, float distance);
	inline int LinePolygonCollision( const vec &a, const vec &b, vec v[], unsigned int count, const vec &normal);
	inline int LineSphereIntersection( const vec &p1, const vec &p2, const vec &center, float radius, double *t1, double *t2 );
	inline int LineSphereIntersectionDir( const vec &p1, const vec &p12, const vec &center, float radius, double *t1, double *t2 );
*/
};
inline char isInfPlus( const float &x);
inline char isInfMinus( const float &x);

inline vec::vec( const vec &a)
{	x = a.x; y = a.y; z = a.z;}
inline vec::vec( const float a, const float b, const float c )
{	x = a; y = b; z = c;}

inline void vec::set( const float a, const float b, const float c)
{	x = a; y = b; z = c;}

inline void vec::clear()
{	x = 0.f; y = 0.f; z = 0.f;}

inline vec vec::operator+() const	{	return *this;}
inline vec vec::operator-() const	{	return vec( -x, -y, -z);}

inline void vec::operator+= ( const vec &a){	x += a.x; y += a.y; z += a.z;}
inline void vec::operator-= ( const vec &a){	x -= a.x; y -= a.y; z -= a.z;}
inline void vec::operator*= ( const vec &a){	x *= a.x; y *= a.y; z *= a.z;}
inline void vec::operator/= ( const vec &a){	x /= a.x; y /= a.y; z /= a.z;}
inline void vec::operator+= ( const float f){	x += f; y += f; z += f; }
inline void vec::operator-= ( const float f){	x -= f; y -= f; z -= f; }
inline void vec::operator*= ( const float f){	x *= f; y *= f; z *= f; }
inline void vec::operator/= ( const float f){	x /= f; y /= f; z /= f; }

inline vec operator+ ( const vec &a, const vec &b ){return vec( a.x+b.x, a.y+b.y, a.z+b.z );}
inline vec operator- ( const vec &a, const vec &b ){return vec( a.x-b.x, a.y-b.y, a.z-b.z );}
inline vec operator* ( const vec &a, const vec &b ){return vec( a.x*b.x, a.y*b.y, a.z*b.z );}
inline vec operator/ ( const vec &a, const vec &b ){return vec( a.x/b.x, a.y/b.y, a.z/b.z );}

inline vec operator+ ( const float f, const vec &a ){return vec( f+a.x, f+a.y, f+a.z );}
inline vec operator- ( const float f, const vec &a ){return vec( f-a.x, f-a.y, f-a.z );}
inline vec operator* ( const float f, const vec &a ){return vec( f*a.x, f*a.y, f*a.z );}
inline vec operator+ ( const vec &a, const float f ){return vec( a.x+f, a.y+f, a.z+f );}
inline vec operator- ( const vec &a, const float f ){return vec( a.x-f, a.y-f, a.z-f );}
inline vec operator* ( const vec &a, const float f ){return vec( a.x*f, a.y*f, a.z*f );}
inline vec operator/ ( const vec &a, const float f ){float f_ = 1.0f/f;return vec( a.x*f_, a.y*f_, a.z*f_ );}

inline int operator< ( const vec &a, const vec &b ){return ((a.x< b.x)&&(a.y< b.y)&&(a.z< b.z));}
inline int operator<=( const vec &a, const vec &b ){return ((a.x<=b.x)&&(a.y<=b.y)&&(a.z<=b.z));}
inline int operator> ( const vec &a, const vec &b ){return ((a.x> b.x)&&(a.y> b.y)&&(a.z> b.z));}
inline int operator>=( const vec &a, const vec &b ){return ((a.x>=b.x)&&(a.y>=b.y)&&(a.z>=b.z));}
inline int operator==( const vec &a, const vec &b ){return ((a.x==b.x)&&(a.y==b.y)&&(a.z==b.z));}
inline int operator!=( const vec &a, const vec &b ){return ((a.x!=b.x)||(a.y!=b.y)||(a.z!=b.z));}

inline vec CROSS ( const vec &a, const vec &b )
{
	return vec( a.y*b.z - a.z*b.y,			//  i   j   k
				a.z*b.x - a.x*b.z,			// a.x a.y a.z
				a.x*b.y - a.y*b.x );		// b.x b.y b.z
}
inline vec vec::CROSS( const vec &a) const
{
	return vec( y*a.z - z*a.y,
				z*a.x - x*a.z,
				x*a.y - y*a.x );
}

inline float DOT3( const vec &a, const vec &b )
{
	return (a.x*b.x + a.y*b.y + a.z*b.z);
}
inline float vec::DOT3( const vec &a) const
{
	return (x*a.x + y*a.y +z*a.z);
}

inline float vec::Length() const
{	
	return (float) sqrt(x*x + y*y + z*z);
}

inline float vec::Length2() const
{	
	return (float) x*x + y*y + z*z;
}

inline void vec::RotX( const float angle)
{
	float s = (float) sin( PI180*angle );
	float c = (float) cos( PI180*angle );
	float Y=y;
	y =  y*c - z*s;
	z =  Y*s + z*c;
}
inline void vec::RotY( const float angle)
{
	float s = (float) sin( PI180*angle );
	float c = (float) cos( PI180*angle );
	float X=x;
	x =  x*c + z*s;
	z = -X*s + z*c;
}
inline void vec::RotZ( const float angle)
{
	float s = (float) sin( PI180*angle );
	float c = (float) cos( PI180*angle );
	float X=x;
	x =  x*c - y*s;
	y =  X*s + y*c;
}
inline void vec::Normalize()
{
	float length;

	length = (float) sqrt(x*x + y*y + z*z);
	if (length != 0) {
		x /= length;
		y /= length;
		z /= length;
	} else {
		x = 0;
		y = 0;
		z = 0;
	}
}

//	Rotate a vector by angle theta around an arbitrary axis r
//	Positive angles are anticlockwise looking down the axis  towards the origin.
inline void vec::ArbitraryRotate(float theta, vec r)	// Rotate a vector by angle theta around an arbitrary axis r
{
	vec p(x,y,z);
	float costheta,sintheta;

	theta*=PI180;
	r.Normalize();
	costheta = (float)cos(theta);
	sintheta = (float)sin(theta);

	x  = (costheta + (1 - costheta) * r.x * r.x) * p.x;
	x += ((1 - costheta) * r.x * r.y - r.z * sintheta) * p.y;
	x += ((1 - costheta) * r.x * r.z + r.y * sintheta) * p.z;

	y  = ((1 - costheta) * r.x * r.y + r.z * sintheta) * p.x;
	y += (costheta + (1 - costheta) * r.y * r.y) * p.y;
	y += ((1 - costheta) * r.y * r.z - r.x * sintheta) * p.z;

	z  = ((1 - costheta) * r.x * r.z - r.y * sintheta) * p.x;
	z += ((1 - costheta) * r.y * r.z + r.x * sintheta) * p.y;
	z += (costheta + (1 - costheta) * r.z * r.z) * p.z;
}

inline vec Normalize( const vec &a){	return a/a.Length();}
inline float Length( const vec &a ){	return a.Length();}
inline float Length2( const vec &a ){	return a.Length2();}
inline float Distance ( const vec &a, const vec &b){	return Length( a-b );}
inline float Distance2 ( const vec &a, const vec &b){	return Length2( a-b );}

inline vec MakeNormal( const vec &a, const vec &b, const vec &c)
{
	return Normalize(CROSS( b - a, c - a));
}

inline vec ReflectedRay( const vec &in, const vec &n)	// in from eye to vertex
{
	return (in - 2.f*n*DOT3( in, n));
}

inline float PlaneDistance( const vec &normal, const vec &point)
{
	return -DOT3( normal, point);	//	D = - Ax + By + Cz;   plane equation is Ax + By + Cz + D = 0
}
inline float PlanePointDelta( const vec &normal, float distance, const vec &point)
{
	// if point is on normal side, return positive value, else return negative value
	return DOT3( normal, point) + distance;			// d = Ax + By + Cz + D
}
inline float PlanePointDistance( const vec &normal, float distance, const vec &point)
{	
	// return always positive value
	return ABS(DOT3( normal, point) + distance);			// distance = ABS(Ax + By + Cz + D)
}

// najblisi bod na usecke a,b k bodu point
inline vec ClosestPointOnLine( const vec &a, const vec &b, const vec &point)
{
	vec ap = point - a;			// vector from a to point
	vec ab = b-a;				// vector from a to b
	float d = ab.Length();		// = Distance( a, b);  // length Line
	ab.Normalize();				// normalized direction vector from a to b
	float t = DOT3( ab, ap);	// ab is unit vector, t is distance from a to point projected on line ab
	if( t<= 0) return a;		// point projected on line ab is out of line, closest point on line is a
	if( t>= d) return b;		// point projected on line ab is out of line, closest point on line is b
	return a+t*ab;				// point away from a on t length in direction ab
}
inline int SpherePlaneCollision( const vec &center, float radius, const vec &normal, float distance, float *distanceCenterToPlane )
{
	*distanceCenterToPlane = DOT3( normal, center) + distance;			// d = Ax + By + Cz + D
	if( fabs(*distanceCenterToPlane) < radius) return 1;	// sphere intersected plane
//	if( *distanceCenterToPlane<= -radius) return 2;			// sphere behind plane
//	else return 0;											// sphere is front plane
	return 0;												// sphere is front or behind plane
}

inline double AngleBetweenVectors( const vec &a, const vec &b)
{							
	double angle = acos( DOT3( a, b) / ( a.Length()*b.Length() ) );
//	if(_isnan(angle))return 0;		// Here we make sure that the angle is not a -1.#IND0000000 number, which means indefinate
	if( isInfPlus( (float)angle)||isInfMinus( (float)angle) )return 0;
	return( angle );
}

// return, t=0 a, t=1 b
inline vec Interpolate( vec a, const vec &b, float t)
{
	if(a==b)return a;
	float angle = (float)AngleBetweenVectors( a, b);
	a.ArbitraryRotate( angle*t/PI180, CROSS( a, b));
	return a;
}

inline int PointInsidePolygon( const vec &point, const vec &a, const vec &b, const vec &c)
{
	vec pa = a-point;
	vec pb = b-point;
	vec pc = c-point;
	double angle  = AngleBetweenVectors( pa, pb);
	angle += AngleBetweenVectors( pb, pc);
	angle += AngleBetweenVectors( pc, pa);

	const double MATCH_FACTOR = 0.99;			// Used to cover up the error in floating point
	if(angle >= (MATCH_FACTOR * (2.0 * PI)) )	// If the angle is greater than 2 PI, (360 degrees)
		return 1;							// The point is inside of the polygon
	return 0;								// If you get here, it obviously wasn't inside the polygon, so Return FALSE
}

inline int PointInsidePolygon( const vec &point, const vec v[], unsigned int count)
{
	double angle = 0.0;
	
	for (int i = 0; (unsigned int)i<count; i++)				// Go in a circle to each vertex and get the angle between
		angle += AngleBetweenVectors( v[i] - point, v[(i+1)%count] - point);
											
	const double MATCH_FACTOR = 0.99;			// Used to cover up the error in floating point
	if(angle >= (MATCH_FACTOR * (2.0 * PI)) )	// If the angle is greater than 2 PI, (360 degrees)