#include	"GraphicsGems.h"

/* ----	intell - Intersect a ray with an ellipsoid. -------------------	*/
/*									*/
/*									*/
/*	Description:							*/
/*	    Intell determines the intersection of a ray with an		*/
/*	    ellipsoid.							*/
/*									*/
/*	On entry:							*/
/*	    raybase = The coordinate defining the base of the		*/
/*		      intersecting ray.					*/
/*	    raycos  = The direction cosines of the above ray.		*/
/*	    center  = The location of the center of the ellispoid.	*/
/*	    radius  = The radius of a bounding sphere.			*/
/*	    Q       = A 4x4 coefficient matrix representing the quadric	*/
/*		      surface.						*/
/*									*/
/*	On return:							*/
/*	    rin     = The entering distance of the intersection.	*/
/*	    rout    = The leaving  distance of the intersection.	*/
/*									*/
/*	Returns:  True if the ray intersects the ellipsoid.		*/
/*									*/
/* --------------------------------------------------------------------	*/


boolean	intell	(raybase,raycos,center,radius,Q,rin,rout)

	Point3	raybase;		/* Base of the intersection ray	*/
	Vector3	raycos;			/* Direction cosines of the ray	*/
	Vector3	center;			/* Base point for the quadric	*/
	double	radius;			/* Bounding radius		*/
	Matrix4	Q;			/* Quadric coefficient matrix	*/
	double	*rin;			/* Entering distance		*/
	double	*rout;			/* Leaving  distance		*/

{
	boolean	hit;			/* True if ray intersects quad.	*/
	double	qa, qb, qc, qd, qe,	/* Coefficient matrix terms	*/
		qf, qg, qh, qi, qj;
	double	k2, k1, k0;		/* Quadric coefficients		*/
	double	disc;
	Vector3	D;			/* Ray base to ellipsoid base	*/
	double	d;			/* Shortest distance between	*/
					/* ray and ellipsoid		*/


/*	Compute the minimal distance between the ray and the center	*/
/*	of the ellipsoid.						*/

	V3Sub (&center,&raybase,&D);
	d   = V3Length (&D);
	V3Normalize (&D);
	d   = d * sqrt (1. - V3Dot(&D,&raycos));
	hit  = (d <= radius);
  
	if  (!hit) return (hit);	    /* If outside of sphere	*/


	qa = Q.element[0][0];
	qb = Q.element[0][1] + Q.element[1][0];
	qc = Q.element[0][2] + Q.element[2][0];
	qd = Q.element[0][3] + Q.element[3][0];
	qe = Q.element[1][1];
	qf = Q.element[1][2] + Q.element[2][1];
	qg = Q.element[1][3] + Q.element[3][1];
	qh = Q.element[2][2];
	qi = Q.element[2][3] + Q.element[3][2];
	qj = Q.element[3][3];

	k2 = raycos.x * (raycos.x*qa + raycos.y*qb) +
	     raycos.y * (raycos.y*qe + raycos.z*qf) +
	     raycos.z * (raycos.z*qh + raycos.z*qc);
	k1 = raycos.x * ((raybase.x*2.*qa + raybase.y*qb + raybase.z*qc) + qd) +
	     raycos.y * ((raybase.x*qb + raybase.y*2.*qe + raybase.z*qf) + qg) +
	     raycos.z * ((raybase.x*qc + raybase.y*qg + raybase.z*2.*qh) + qi);
	k0 = raybase.x * (raybase.x*qa + raybase.y*qb + raybase.z*qc + qd) +
	     raybase.y * (               raybase.y*qe + raybase.z*qf + qg) +
	     raybase.z * (                              raybase.z*qh + qi) +
	     qj;

	disc = k1*k1 - 4.*k2*k0;
  
	hit  = (disc >= 0.0);

	if  (hit) { 				/* If ray hits quadric	*/
	    disc  =  sqrt(disc);
	    *rin  = (-k1 - disc) / (2.*k2);	/*    entering distance	*/
	    *rout = (-k1 + disc) / (2.*k2);	/*    leaving distance	*/
	}
  
	return (hit);
}