vector.h

030_Particle_system using opengl

		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 EdgeSphereCollision( const vec &center, float radius, const vec &a, const vec &b, const vec &c)
{
	if( Distance(center, ClosestPointOnLine( a, b, center)) < radius) return 1;
	if( Distance(center, ClosestPointOnLine( b, c, center)) < radius) return 1;
	if( Distance(center, ClosestPointOnLine( c, a, center)) < radius) return 1;
	return 0;
}

inline int EdgeSphereCollision( const vec &center, float radius, vec v[], unsigned int count)
{
	for(int i = 0; (unsigned int)i<count; i++)
		if( Distance(center, ClosestPointOnLine(v[i], v[(i+1)%count], center)) < radius) return 1;
	return 0;
}

inline int SpherePolygonCollision(const vec &center, float radius, vec v[], unsigned int count, const vec &normal)
{
	float distance = PlaneDistance( normal, v[0] );
	float distanceCenterToPlane;

	if(SpherePlaneCollision( center, radius, normal, distance, &distanceCenterToPlane) == 1) 
	{
		vec centerProjOnPlane = center - distanceCenterToPlane*normal;			// center project on plane

		if( PointInsidePolygon( centerProjOnPlane, v, count ) ||
			EdgeSphereCollision( center, radius, v, count ) )
		{
			return 1;
		}
	}
	return 0;
}

inline int SpherePolygonCollision(const vec &center, float radius, const vec &a, const vec &b, const vec &c, const vec &normal)
{
	float distance = PlaneDistance( normal, a );
	float distanceCenterToPlane;

	if(SpherePlaneCollision( center, radius, normal, distance, &distanceCenterToPlane) == 1) 
	{
		vec centerProjOnPlane = center - distanceCenterToPlane*normal;			// center project on plane

		if( PointInsidePolygon( centerProjOnPlane, a, b, c )  ||
			EdgeSphereCollision( center, radius, a, b, c ) )
		{
			return 1;
		}
	}
	return 0;
}

inline int SpherePolygonCollisionRadius05(const vec &center, float radius, const vec &a, const vec &b, const vec &c, const vec &normal)
{
	float distance = PlaneDistance( normal, a );
	float distanceCenterToPlane;

	if(SpherePlaneCollision( center, radius, normal, distance, &distanceCenterToPlane) == 1) 
	{
		vec centerProjOnPlane = center - distanceCenterToPlane*normal;			// center project on plane

		if( PointInsidePolygon( centerProjOnPlane, a, b, c )  ||
			EdgeSphereCollision( center, 0.5f*radius, a, b, c ) )
		{
			return 1;
		}
	}
	return 0;
}

inline int LinePlaneCollision( const vec &a, const vec &b, const vec &normal, float distance)
{
	float distance_a_plane = DOT3( normal, a) + distance;
	float distance_b_plane = DOT3( normal, b) + distance;
	if( distance_a_plane*distance_b_plane < 0) return 1;	// must have differently sign
	return 0;
}

inline vec LinePlaneIntersectionDir( const vec &a, const vec &ab, const vec &normal, float distance)
{
	vec ab_norm = Normalize( ab);
	float distance_plane_a = DOT3( normal, a) + distance;
	float divisor = DOT3( normal, ab_norm);
	if(divisor==0.f) return a;				// line lie on plane
	float t = distance_plane_a/divisor;
	return a-t*ab_norm;
}

inline float LinePlaneIntersectionDirParameter( const vec &a, const vec &ab, const vec &normal, float distance)
{
	vec ab_norm = Normalize( ab);
	float distance_plane_a = DOT3( normal, a) + distance;
	float divisor = DOT3( normal, ab_norm);
	if(divisor==0.f) return 0.f;				// line lie on plane
	float t = distance_plane_a/divisor;
	return -t/ab.Length();						// for point a t=0, for a+ab t=1.0f 
}

inline vec LinePlaneIntersection( const vec &a, const vec &b, const vec &normal, float distance)
{
	return LinePlaneIntersectionDir( a, b-a, normal, distance );
}

inline int LinePolygonCollision( const vec &a, const vec &b, vec v[], unsigned int count, const vec &normal)
{
	float distance = PlaneDistance( normal, v[0]);
	if( !LinePlaneCollision( a, b, normal, distance) ) return 0;

	vec Intersection = LinePlaneIntersection( a, b, normal, distance);
	if( PointInsidePolygon( Intersection, v, count) ) return 1;
	return 0;
}

inline int LineSphereIntersectionDir( const vec &p1, const vec &p12, const vec &center, float radius, double *t1, double *t2 )
{
//	vec	p12 = p2-p1;
	double a = DOT3(p12, p12);		// a is always positive
	double b = 2.0*DOT3( p12, p1 - center);
	double c = DOT3( center, center) + DOT3( p1, p1) - 2.0*DOT3( center, p1) - radius*radius;
	double diskriminant = b*b - 4.0*a*c;

	if(diskriminant<0 || a==0){	*t1 = 0;*t2 = 0;return 0;}
	if(diskriminant==0)
	{
		*t1 = ( -b/(2.0*a) );
		*t2 = *t1;
		return 1;
	}
	double s_diskriminant = sqrt(diskriminant);
	*t1 = ( (-b + s_diskriminant)/(2.0*a) );
	*t2 = ( (-b - s_diskriminant)/(2.0*a) );
	// because a is always positive, t1 > t2
	return 2;
}

inline int LineSphereIntersection( const vec &p1, const vec &p2, const vec &center, float radius, double *t1, double *t2 )
{
	return LineSphereIntersectionDir( p1, p2-p1, center, radius, t1, t2);
}
/*  angle1
   in\  A n-nomala
      \ |     n1 - absolute refraction index of material 1
       V|                        sin(angle1)   v1   c/n1   n2        1.0
   -----+------  nr = n2/n1      ----------- = -- = ---- = -- = nr = ---
         \                       sin(angle2)   v2   c/n2   n1        nri
         |    n2 - absolute refraction index of material 2
         \                       n1*sin(angle1) = n2*sin(angle2)
   angle2 | T - transmission vector
in - vector from eye to vertex, must be normalized
n - normal vector
nr - relative refraction index, water nr = 1.33
nri - nri = 1.f/nr
  T = nri*in - nri*(in*n)*n - sqrt( 1.0-nri^2*(1.0-(in,n)^2) )*n   */
// *i - call with nri = 1/relative refraction index
// *n - vector in is Normalized
inline vec RefractedRayin(const vec &in, const vec &n, float nri)	// in from eye to vertex
{
//	float nri=1.f/nr;
	float dot_in_n = DOT3(in,n);
	float a = (float)sqrt( 1.f-nri*nri*(1.f-dot_in_n*dot_in_n) );
	return nri*in - (a+nri*dot_in_n)*n;
}
inline vec RefractedRayi( vec &in, const vec &n, float nri)		// normalize in
{	return RefractedRayin( Normalize(in), n, nri);}
inline vec RefractedRay( vec &in, const vec &n, float nr)
{	return RefractedRayin( Normalize(in), n, 1.f/nr);}

// work also if DOT3( in, n) > 0.f
inline vec RefractedRayic( vec &in, const vec &n, float nri)	// in from eye to vertex
{
	if( DOT3( in, n) > 0.f)
		return RefractedRayi( in, -n, 1.f/nri);
	else 
		return RefractedRayi( in, n, nri);
}

inline float FresnelF0( float nr)
{
	float nri=1.f/nr;
	float f0,f0_;
	f0 = 1-nri;		//      (1-nri)^2
	f0 *= f0;		// f0 = ----------
	f0_ = 1+nri;	//      (1+nri)^2
	f0_ *= f0_;
	f0 /= f0_;
	return f0;
}
// Fresnel
inline float Fresnelin( const vec &in, const vec &n, float nri, float f0)
{
	float dot = DOT3( n, -in);
	float dot2 = dot*dot;
	dot *= dot2*dot2;		// dot = (n*v)^5
	return f0+(1.f-f0)*dot;
}
inline float Fresneli( const vec &in, const vec &n, float nri, float f0)
{	return Fresnelin( Normalize(in), n, nri, f0);}
inline float Fresnel( const vec &in, const vec &n, float nr, float f0)
{	return Fresnelin( Normalize(in), n, 1.f/nr, f0);}

inline float Fresnelin_1( const vec &in, const vec &n, float nri, float f0)
{
	float dot = DOT3( n, -in);	// dot = (n*v)^1
	return f0+(1.f-f0)*dot;
}
inline float Fresnelin_3( const vec &in, const vec &n, float nri, float f0)
{
	float dot = DOT3( n, -in);
	dot *= dot*dot;				// dot = (n*v)^3
	return f0+(1.f-f0)*dot;
}

inline char isInfPlus( const float &x)
{
	float a;
	*(unsigned int*)&a = (unsigned int)0x7f800000;
	return ( x==a);
//	return (*(unsigned int*) &x==0x7f800000);
}

inline float Fresnel2in( const vec &in, const vec &n, float nri)
{
	float k,g;
	
	k = DOT3( in, n);
	g = (float)sqrt(nri*nri+k*k-1);

	float gmk = g-k;
	float gpk = g+k;

	float a = k*gpk-1;
	float b = k*gmk+1;
	return (gmk*gmk/(2*gpk*gpk)) *(1+(a*a)/(b*b));
}
inline float Fresnel2( const vec &in, const vec &n, float nr)
{	return Fresnel2in( Normalize(in), n, 1.f/nr);}



inline char isInfMinus( const float &x)
{
	float a;
	*(unsigned int*)&a = (unsigned int)0xff800000;
	return ( x==a);
//	return (*(unsigned int*) &x==0xff800000);
}


class vec2  
{
public:
	union
	{
		struct
		{
			float x, y;
		};
		float v[2];
	};
	void set( float x_, float y_){ x=x_; y=y_;}
};

class vec4 
{
public:
	union
	{
		struct
		{
			float x, y, z, w;
		};
		float v[4];
	};

	vec4(){}
	inline vec4( const vec4 &a);
	inline vec4( const vec &in);
	inline vec4( const float a, const float b, const float c, const float d = 1.f);
	inline void set( const float a, const float b, const float c, const float d = 1.f);
	inline void clear();
};

inline vec4::vec4( const vec4 &a)
{	x = a.x; y = a.y; z = a.z; w = a.w;}
inline vec4::vec4( const vec &in)
{	x = in.x; y = in.y; z = in.z; w = 1.f;}
inline vec4::vec4( const float a, const float b, const float c, const float d)
{	x = a; y = b; z = c; w = d;}
inline void vec4::set( const float a, const float b, const float c, const float d)
{	x = a; y = b; z = c; w = d;}
inline void vec4::clear()
{	x = 0; y = 0; z = 0; w = 1;}

inline vec::vec( const vec4 &a)
{	x = a.x; y = a.y; z = a.z;}

inline vec::vec( const vec2 &a)
{	x = a.x; y = a.y; z = 0.f;}


inline vec Bspline( const vec &p0, const vec &p1, const vec &p2, const vec &p3, float t)
{
	float t2=t*t;		// t^2
	float c3=t2*t;		// t^3
	float c0,c1,c2;

	c0 = 1.f - t;	c0 *= c0*c0;		// c0 = (1-t)^3
	c1 = 3.f*c3 - 6.f*t2 + 4;			// c1 = 3*t^3 - 6*t^2 + 4
	c2 =-3.f*c3 + 3.f*t2 + 3.f*t + 1.f;	// c2 =-3*t^3 + 3*t^2 + 3*t + 1
//	c3 = t*t*t;
//	vec v = c0*p0 + c1*p1 + c2*p2 + c3*p3;
//	v *= 0.166666666666666666666666666666667f;	// 1/6
	return 0.166666666666666666666666666666667f*(c0*p0 + c1*p1 + c2*p2 + c3*p3);
}

inline vec BsplineDeriv( const vec &p0, const vec &p1, const vec &p2, const vec &p3, float t)
{
	float t2=t*t;		// t^2
	float c0,c1,c2,c3;
//	(1-t)^3 = (1-t)*(1-t-t+t^2) = (1-t)*(1-2*t+t^2) = (1-2*t+t^2 -t +2*t^2 - t^3) = (1-3*t+3*t^2-t^3)
	c0 =-3.f*t2 + 6.f*t - 3;			// c0 =-3*t^2 + 6*t - 3
	c1 = 9.f*t2 -12.f*t;				// c1 = 9*t^2 -12*t
	c2 =-9.f*t2 + 6.f*t + 3;			// c2 =-9*t^2 + 6*t + 3
	c3 = 3.f*t*t;						// c3 = 3*t^2
	return 0.166666666666666666666666666666667f*(c0*p0 + c1*p1 + c2*p2 + c3*p3);
}

inline vec HermiteCurve( const vec &p0, const vec &p1, const vec &t0, const vec &t1, float t)
{
	float t2 = t*t;
	float t3 = t2*t;
	float h1 =  2*t3 - 3*t2 + 1;
	float h2 = -2*t3 + 3*t2;
	float h3 =    t3 - 2*t2 + t;
	float h4 =    t3 -   t2;
	return h1*p0 + h2*p1 + h3*t0 + h4*t1;
}

inline vec HermiteCurveDeriv( const vec &p0, const vec &p1, const vec &t0, const vec &t1, float t)
{
	float t2 = t*t;
	float h1 =  6*t2 - 6*t;
	float h2 = -6*t2 + 6*t;
	float h3 =  3*t2 - 4*t + 1;
	float h4 =  3*t2 - 2*t;
	return h1*p0 + h2*p1 + h3*t0 + h4*t1;
}

#endif