/*
This file contains intersection of ray with Canonical Cylinder, Cone and Ring.
1)
	Intersection of a ray with canonical Cylinder.
	This canonical cylinder has axis coinciding
	with +ve Z axis with one base on origin with
	the Height H and the Radius R.

	The equation :
		x^2 + y^2 - R^2 = 0

	"ray_cylinder" computes the intersection and
	returns rayparameter t of intesection
		where t > 0 for valid intersection,
	and surface parameter (u,v) of the intersection point
		where v is the height parameter
		and   u is along the circle in counterclock direction.
2.
        Ray  Canonical cone intersection.
        Canonical cone has :
                apex at origin with axis aligned with Z.
                base diameter and height 1.
                So the slant angle is 45 degrees.
        In this Cononical form the surface equation :
                x^2 + y^2 - z^2 = 0
                i.e. the radius of the crosssection circle  is equal
                     to its z.
        "ray_cone" computes the intersection and returns parameter t of
        intersection
                where t > 0 for valid intersection,
        and surface parameter (u,v) of the intersection point
                where v is the height of the canonical cone withing the
                limit and u is along the crosssectional cirlce in
		counterclockwise direction.
3.
	Ray Cononical Ring intersection.
	Canonical ring is defined by the ring with
	center at origin and Z-axis as normal.
	Its further defined by its inner radius(R0), outer radius(R1).
	0 <= R0 < R1.
------
sumant
*/
#include <stdio.h>
#include <math.h>
#include "GraphicsGems.h"
#include "data_structure.h"

static int check(t,ray_start,ray_direction,lo,hi)
double t;
Point3 *ray_start;
Vector3 *ray_direction;
double lo,hi;
{
	double zposn=ray_start->z + ray_direction->z*t;
	return(zposn >= lo && zposn <= hi);
}

double ray_cylinder(ray_start,ray_direction,R,H,tubeflag,pu,pv)
/*
	Tubeflag = 0 for the cylinder whose out surface visible.
		   1 for the cylinder whose inside surface visible.
	The intersection algorithm is same except in the surface
		parameterisation.
*/
Point3 *ray_start;                                      /* Input */
Vector3 *ray_direction;                                 /* Input */
double R,H;						/* Input */
int tubeflag;						/* Input */
double *pu,*pv;						/* Output */
{
	int quadraticRoots();
	int nroots;
	double t = -1.,roots[2];
	double A,B,C;
	C  = ray_start->x * ray_start->x + ray_start->y * ray_start->y
		- R * R;
	B = 2.*(ray_start->x * ray_direction->x
		+ ray_start->y * ray_direction->y);
	A = ray_direction->x * ray_direction->x 
		+ ray_direction->y * ray_direction->y;
	if (nroots = quadraticRoots(A,B,C,roots))
		set_t(0.,H);
	if (t > EPSILON){ /* Compute Surface Parameter */
		Point3 p;
		double u;
		point_on_line(*ray_start,*ray_direction,t,p);
		*pv = p.z/H;
		circle_inverse_mapping(R,p.x,p.y,u);
/*
		*pu = (tubeflag)?(1.-u):u;
*/
		*pu = u;
		return(t);
	}
	return(-1.);
}

double ray_cone(ray_start,ray_direction,distance,cupflag,pu,pv)
Point3 *ray_start;                                      /* Input */
Vector3 *ray_direction;                                 /* Input */
double distance;                                        /* Input */
int cupflag;                                            /* Input */
double *pu,*pv;                                         /* Output */
{
	int quadraticRoots();
	int nroots;
	double t = -1.,roots[2];
	double A,B,C;
	C = ray_start->x * ray_start->x + ray_start->y * ray_start->y
		- ray_start->z * ray_start->z;
	B = 2.*(ray_start->x * ray_direction->x
		+ ray_start->y * ray_direction->y
		- ray_start->z * ray_direction->z);
	A = ray_direction->x * ray_direction->x
		+ ray_direction->y * ray_direction->y
		- ray_direction->z * ray_direction->z;
	if (nroots = quadraticRoots(A,B,C,roots))
		set_t(distance,1.);
        if (t > EPSILON){ /* Compute Surface Parameter */
                Point3 p;
                double u;
                point_on_line(*ray_start,*ray_direction,t,p);
                *pv = (p.z-distance)/(1.-distance);
                circle_inverse_mapping(p.z,p.x,p.y,u);
/*
                *pu = (cupflag)?(1.-u):u;
*/
		*pu = u;
		return(t);
        }
	return(-1.);
}

double ray_ring(ray_start,ray_direction,R0,R1,pu,pv)
Point3 *ray_start;                                      /* Input */
Vector3 *ray_direction;                                 /* Input */
double R0,R1;						/* Input */
double *pu,*pv;						/* Output */
{
	double t = -1.;
	if (fabs(ray_direction->z) > EPSILON){ /* Possible intersection */
		Point3 p;
		double radius;
		double value = -ray_start->z/ray_direction->z;
		if (value > EPSILON){
			point_on_line(*ray_start,*ray_direction,value,p);
			radius = V3SquaredLength(&p);
			if (radius >= R0*R0 && radius <= R1*R1){
				t = value;
				radius = sqrt(radius);
				*pv = (radius - R0)/(R1-R0);
				circle_inverse_mapping(radius,p.x,p.y,*pu);
				return(t);
			}
		}
	} 
	return(-1.);
}