/*
 * rotspline.c
 *
 * Copyright (C) 1991, Gerald Iles.
 * All rights reserved.
 *
 * This software may be freely copied, modified, and redistributed
 * provided that this copyright notice is preserved on all copies.
 *
 * You may not distribute this software, in whole or in part, as part of
 * any commercial product without the express consent of the authors.
 *
 * There is no warranty or other guarantee of fitness of this software
 * for any purpose.  It is provided solely "as is".
 *
 * $Id$
 *
 * $Log$
 */
#include "geom.h"
#include "rotspline.h"

static Methods *iRotsplineMethods = NULL;
static char rotsplineName[] = "rotspline";

unsigned long RotsplineTests, RotsplineHits;

Rotspline *
RotsplineCreate(bot, x0, d0, top, xw, dw, coeffs)
Vector *bot, *top;
Float x0, xw, d0, dw, *coeffs;
{
	Rotspline *rotspline;
	Float w;
	Vector axis, tmp;

	VecSub(*top, *bot, &axis);
	w = VecNormalize(&axis);
	if (w < EPSILON) {
		RLerror(RL_WARN, "Degenerate rotspline.\n");
		return (Rotspline *)NULL;
	}

	rotspline = (Rotspline *)Malloc(sizeof(Rotspline));

	rotspline->w = w;

	if (coeffs == NULL)
	{
		rotspline->k = (Float *)Calloc(4, sizeof(Float));
		rotspline->k[3] = 2.0 * (x0 * x0 - xw * xw + w * x0 * d0 + w * xw * dw) / (w * w * w);
		rotspline->k[2] = (3.0 * (xw * xw - x0 * x0) - 4.0*w*x0*d0 - 2.0*w*xw*dw) / (w * w);
		rotspline->k[1] = 2.0 * x0 * d0;
		rotspline->k[0] = x0 * x0;
	}
	else
		rotspline->k = coeffs;
	CoordSysTransform(bot, &axis, 1.0, 1.0, &rotspline->trans);

	return rotspline;
}

Methods *
RotsplineMethods()
{
	if (iRotsplineMethods == (Methods *)NULL) {
		iRotsplineMethods = MethodsCreate();
		iRotsplineMethods->name = RotsplineName;
		iRotsplineMethods->create = (GeomCreateFunc *)RotsplineCreate;
		iRotsplineMethods->methods = RotsplineMethods;
		iRotsplineMethods->intersect = RotsplineIntersect;
		iRotsplineMethods->normal = RotsplineNormal;
		iRotsplineMethods->uv = RotsplineUV;
		iRotsplineMethods->bounds = RotsplineBounds;
		iRotsplineMethods->stats = RotsplineStats;
		iRotsplineMethods->checkbounds = TRUE;
		iRotsplineMethods->closed = FALSE;
	}
	return iRotsplineMethods;
}

int
RotsplineIntersect(rotspline, ray, mindist, maxdist)
Rotspline *rotspline;
Ray *ray;
Float mindist, *maxdist;
{
	double c[4], s[3], pp, pd, pz, dz, tmp, dist;
	Ray newray;
	Vector ndir, npos;
	int num, i;

	RotsplineTests++;

	newray = *ray;
	(void) RayTransform(&newray, &rotspline->trans.itrans);
	ndir = newray.dir;
	npos = newray.pos;

	pp = npos.x*npos.x + npos.y*npos.y +npos.z*npos.z;
	pd = npos.x*ndir.x + npos.y*ndir.y +npos.z*ndir.z;
	pz = npos.z;
	dz = ndir.z;

	c[3] = rotspline->k[3] * dz * dz * dz;
	c[2] = dz * dz *(3.0 * rotspline->k[3] * pz + rotspline->k[2] + 1.0) - 1.0;
	c[1] = dz *(pz *(3.0 * rotspline->k[3] * pz + 2.0 * rotspline->k[2] + 2.0) + rotspline->k[1]) - 2.0 * pd;
	c[0] = pz *(pz *(rotspline->k[3] * pz + rotspline->k[2] +1.0) + rotspline->k[1]) + rotspline->k[0] - pp;

	if (fabs(c[3]) > EPSILON)
		num = SolveCubic(c, s);
	else if (fabs(c[2]) > EPSILON)
		num = SolveQuadric(c, s);
	else if (fabs(c[1]) > EPSILON)
	{
		s[0] = -c[0]/c[1];
		num = 1;
	}
	else
		num = 0;
		
	if (num==0) return FALSE;
	dist = 0.0;
	for(i=0;i<num;i++)
	{
		if (s[i] > mindist)
		{
			tmp = pz + dz * s[i];
			if (tmp > 0.0 && tmp < rotspline->w)
			{
				if (dist == 0.0)
					dist = s[i];
				else if (s[i] < dist)
					dist = s[i];
			}
		}
	}
	if (dist > mindist && dist < *maxdist) {
		*maxdist = dist;
		RotsplineHits++;
		return TRUE;
	}
	return FALSE;
}

int
RotsplineNormal(rotspline, pos, nrm, gnrm)
Rotspline *rotspline;
Vector *pos, *nrm, *gnrm;
{
	Vector npos;
	Float du, r;

	npos = *pos;
	PointTransform(&npos, &rotspline->trans.itrans);
	
	du = rotspline->k[3] * 3.0 * npos.z * npos.z
	   + rotspline->k[2] * 2.0 * npos.z
	   + rotspline->k[1];

	r = sqrt(npos.x * npos.x + npos.y * npos.y);

	nrm->x = r * npos.x;
	nrm->y = r * npos.y;
	nrm->z = -du / 2.0;

	NormalTransform(nrm, &rotspline->trans.itrans);
	*gnrm = *nrm;
	return FALSE;
}

void
RotsplineUV(rotspline, pos, norm, uv, dpdu, dpdv)
Rotspline *rotspline;
Vector *pos, *norm, *dpdu, *dpdv;
Vec2d *uv;
{
	Vector npos;
	Float val, d;

	npos = *pos;
	PointTransform(&npos, &rotspline->trans.itrans);

	uv->v = npos.z / rotspline->w;
	d = sqrt(npos.x * npos.x + npos.y * npos.y);
	if (d < EPSILON)
		uv->u = 0.;
	else
	{
		val = npos.x / d;
		if (val > 1.)
			uv->u = 0.;
		else if (val < -1.)
			uv->u = 0.5;
		else
		{
			uv->u = acos(val) / TWOPI;
			if (npos.y < 0.)
				uv->u = 1. - uv->u;
		}
	}

	if (dpdu)
	{
		dpdu->x = -npos.y;
		dpdu->y = npos.x;
		dpdu->z = 0.;
		dpdv->x = rotspline->k[3] * 3.0 * npos.z * npos.z
				+ rotspline->k[2] * 2.0 * npos.z
				+ rotspline->k[1];
		dpdv->y = 0.;
		dpdv->z = (d < EPSILON)?0.:d;
		VecTransform(dpdu, &rotspline->trans.trans);
		VecTransform(dpdv, &rotspline->trans.trans);
		(void)VecNormalize(dpdu);
		(void)VecNormalize(dpdv);
	}
}

void
RotsplineBounds(rotspline, bounds)
Rotspline *rotspline;
Float bounds[2][3];
{
	Float r0, rw, rmax, d, z, rz;

	r0 = (rotspline->k[0] < EPSILON ? 0.0 : sqrt(rotspline->k[0]));
	rw = EvalSpline(rotspline->k, rotspline->w);
	rw = (rw < EPSILON ? 0.0 : sqrt(rw));
	rmax = max(r0, rw);
	d = 4.*rotspline->k[2]*rotspline->k[2] - 12.*rotspline->k[3]*rotspline->k[1];
	if (d > 0.)
	{
		z = (-2. * rotspline->k[2] + sqrt(d))/(6. * rotspline->k[3]);
		if (z > 0. && z < rotspline->w)
		{
			rz = EvalSpline(rotspline->k, z);
			rz = (rz < EPSILON ? 0.0 : sqrt(rz));
			rmax = max(rmax, rz);
		}
		z = (-2. * rotspline->k[2] - sqrt(d))/(6. * rotspline->k[3]);
		if (z > 0. && z < rotspline->w)
		{
			rz = EvalSpline(rotspline->k, z);
			rz = (rz < EPSILON ? 0.0 : sqrt(rz));
			rmax = max(rmax, rz);
		}
	}
	bounds[LOW][X] = bounds[LOW][Y] = -rmax;
	bounds[HIGH][X] = bounds[HIGH][Y] = rmax;
	bounds[LOW][Z] = 0.;
	bounds[HIGH][Z] = rotspline->w;
	BoundsTransform(&rotspline->trans.trans, bounds);
}

char *
RotsplineName()
{
	return rotsplineName;
}

void
RotsplineStats(tests, hits)
unsigned long *tests, *hits;
{
	*tests = RotsplineTests;
	*hits = RotsplineHits;
}