/*
        Author: S.N.pattanaik
        Date : 25 March 1991.

	NOTE : Ray Direction must be NORMALISED.
*/

#include <stdio.h>
#include <math.h>
#include "GraphicsGems.h"
#include "data_structure.h"

double sphere_intersection(ray_start,ray_direction,
				sphere_radius, sphere_center)
/*
        This routine returns the parameter of intersection
        of a ray defined as (ray_start+t.ray_direction) and
        a sphere.
NOTE : Assumption of the Ray Direction being NORMALISED has been made.
                                             ----------
*/
Point3 *ray_start;					/* Input */
Vector3 *ray_direction;					/* Input */
double sphere_radius;					/* Input */
Point3 *sphere_center;					/* Input */
{
	int quadraticRoots();
	double t = -1.0,roots[2];
	double dx=ray_start->x-sphere_center->x;	
	double dy=ray_start->y-sphere_center->y;	
	double dz=ray_start->z-sphere_center->z;	
	double B=2.0*(ray_direction->x*dx+
		ray_direction->y*dy+ray_direction->z*dz);
	double C=(dx*dx+dy*dy+dz*dz)-sphere_radius*sphere_radius;
	int nroots = quadraticRoots((double)1.,B,C,roots);
#define check(r,ray_start,ray_direction,lo,hi) 1
	if (nroots)
		set_t(0.,1.); /* Dummy parameters. */
#undef check
/*
	double discriminant=B*B-4.0*C;
	if (discriminant >= 0.0){
		if (discriminant>0.0) discriminant = sqrt(discriminant);
		t=(-B-discriminant)/2.0;
		if ((discriminant>0.0) && (t <= 0.0)){
			double t1;
			t1 = (-B+discriminant)/2.0;
			if (t > t1) t= t1;
		}
	}
*/
	return(t);
}

double ray_sphere(ray_start,ray_direction,
		sphere_radius,
		sphere_center,
		bubble_flag,
		pu,pv)
/*
	This routine finds out the intersection of ray with Sphere.
	If the ray intersects the sphere then the routine returns 
		the ray parameter t
		and
		surface parameters (u,v) of intersection.
	t > 0 and
	(u,v) are in the range 0 to 1 for a valid intersection.

Source: "Essential Ray Tracing Algorithms" by Eric Haines in
        "An Introduction to RayTracing", page 48.

*/	
Point3 *ray_start;					/* Input */
Vector3 *ray_direction;					/* Input */
double sphere_radius;					/* Input */
Point3 *sphere_center;					/* Input */
int bubble_flag;					/* Input */
double *pu,*pv;						/* Output */
{
	double t = sphere_intersection(ray_start,ray_direction,
			sphere_radius,sphere_center);
        if (t > EPSILON){
                Point3 R,N;
		double u,v,phi;

		point_on_line(*ray_start,*ray_direction,t,R);
		N.x = (R.x - sphere_center->x)/sphere_radius;
		N.y = (R.y - sphere_center->y)/sphere_radius;
		N.z = (R.z - sphere_center->z)/sphere_radius;

		if (N.z >= 1.) phi = 0.;
		else if (N.z <= -1.) phi = PI;
		else phi = acos(N.z);
		v = phi / PI;
		if ((v == 0.) || (v == 1.)) u = 0;
		else{
			double r = sin(phi);
			circle_inverse_mapping(r,N.x,N.y,u)
		}
/*
		*pu = (bubble_flag)?(1.-u):u;
*/
		*pu = u;
		*pv = v;
		return(t);
	}
	return(-1.0);
}