#ifndef lint
static char *sccsid = "@(#)ray.c	1.2 (Steve Hill) 5/24/90";
#endif

/* ray.c
 *
 * Auxiliary maths routines for manipulating rays
 *
 * Routines to create and manipulate rays.  A ray is a point, and
 * a vector that passes through that point.
 * A ray is defined by the parametric equation:
 *
 * (x0 + Vx t, y0 + Vy t, Vz t)
 *
 */

#include <stdio.h>
#include <math.h>

#include "basetype.h"
#include "cartesian.h"
#include "malloc.h"
#include "error.h"
#include "matrix.h"
#include "point.h"
#include "vector.h"
#include "indexlist.h"
#include "ray.h"


/* Ray
 *
 * Create a ray from a point and a vector.
 */

ray_t *
Ray(point, vector)
point_t		*point;
vector_t	*vector;
{
	ray_t	*ray;

	ray = (ray_t *) malloc(sizeof(ray_t));

	if (ray == RayNull)
		FatalError(__FILE__, __LINE__, "Ray: Out of memory");

	ray->point  = point;
	ray->vector = vector;
	ray->list   = IndexListNull;

	return(ray);
}

ray_t *
RaySet(ray, point, vector, list)
ray_t		*ray;
point_t		*point;
vector_t	*vector;
index_list_t	*list;
{
	ray->point  = point;
	ray->vector = vector;
	ray->list   = list;
	return(ray);
}

ray_t *
RayListCopy(ray, point, vector, list)
ray_t		*ray;
point_t		*point;
vector_t	*vector;
index_list_t	*list;
{
	ray->point  = point;
	ray->vector = vector;
	ray->list   = IndexListCopy(list);
	return(ray);
}

void
RayListFree(ray)
ray_t	*ray;
{
	IndexListFree(ray->list);
}


/* RayFree
 *
 * Free a ray (also free the vector and point)
 */

void
RayFree(ray)
ray_t	*ray;
{
	if (ray == RayNull)
		return;

	VectorFree(ray->vector);
	PointFree(ray->point);
	free((char *) ray);
}


/* RayInstantiate
 *
 * Given a parameter and a ray, find the point.
 */

point_t *
RayInstantiate(t, ray)
real_t	t;
ray_t	*ray;
{
	return(Point(
			ray->point->x + t * ray->vector->x,
			ray->point->y + t * ray->vector->y,
			ray->point->z + t * ray->vector->z
		     ));
}


/* RayDistance
 *
 * Calculate parameter of given point along a ray.
 */

real_t
RayDistance(ray, point)
ray_t	*ray;
point_t	*point;
{
	real_t	d;

	if ((d = point->x - ray->point->x) != REAL_ZERO)
	{
		if (ray->vector->x == 0)
			return(REAL_BIG);
		return(d / ray->vector->x);
	}
	if ((d = point->x - ray->point->y) != REAL_ZERO)
	{
		if (ray->vector->y == 0)
			return(REAL_BIG);
		return(d / ray->vector->y);
	}
	if ((d = point->x - ray->point->z) != REAL_ZERO)
	{
		if (ray->vector->z == 0)
			return(REAL_BIG);
		return(d / ray->vector->z);
	}
	return(REAL_ZERO);
}

/* RayPrint
 *
 * Print out a ray to a file.
 */

void
RayPrint(file, ray)
FILE	*file;
ray_t	*ray;
{
	if (ray == RayNull)
	{
		fprintf(file, "Null Ray\n");
		return;
	}

	PointPrint(file, ray->point);
	VectorPrint(file, ray->vector);
}

vector_t *
VectorAddMult(vector1, scale2, vector2, scale3, vector3)
vector_t	*vector1, *vector2, *vector3;
real_t		scale2, scale3;
{
	vector1->x = scale2 * vector2->x + scale3 * vector3->x;
	vector1->y = scale2 * vector2->y + scale3 * vector3->y;
	vector1->z = scale2 * vector2->z + scale3 * vector3->z;

	return(vector1);
}


bool_t
Refract(result, from_index, to_index, incident, normal)
vector_t	*result, *incident, *normal;
real_t		from_index, to_index;
{
	real_t	cos1, cos2, kn, k;
	vector_t	normal1;

	VectorLet(&normal1, normal);
	cos1 = -VectorDot(incident, &normal1);

	if (cos1 < REAL_ZERO)
	{
		VectorNeg(&normal1, &normal1);
		cos1 = -cos1;
	}

	kn = from_index / to_index;
	cos2 = REAL_ONE - kn * kn * (REAL_ONE - cos1 * cos1);
	if (cos2 < REAL_ZERO)
		return(TRUE);
	k = kn * cos1 - sqrt(cos2);
	VectorAddMult(result, kn, incident, k, &normal1);
	return(FALSE);
}