#ifndef lint
static char	*rcs = "$Header: spline.c,v 1.1 88/01/15 13:05:27 simpson Rel $";
#endif
/*
 * $Log:	spline.c,v $
 * Revision 1.1  88/01/15  13:05:27  simpson
 * initial release
 * 
 * Revision 0.1  87/12/18  11:20:40  simpson
 * beta test
 * 
 */
#include <math.h>

static double U[50];

/* Returns a set of (m+1) spline curve points given a set of (n+1) control
 * points.  (t-1) is the degree of the polynomial used for the B-splines.
 * Cubic polynomials (i.e., t=4) are usually sufficient.  [3] designates
 * (x,y,z).  For a discussion of the spline curve algorithm, see the book
 * ``Computer Graphics'', Donald Hearn and M. Pauline Baker, Prentice-Hall,
 * 1986, pp. 200-2.
 */
void spline(controlpoints, n, curvepoints, m, t)
double  controlpoints[][3];
int	n;
double  curvepoints[][3];
int	m;
int	t;
{
    int 	j, k;
    double	N(), temp, u;

    for (j = 0; j <= n + t; j++) {
	if (j < t)
	    U[j] = 0.0;
	else if (t <= j && j <= n)
	    U[j] = j - t + 1;
	else
	    U[j] = n - t + 2;
    }
    for (j = 0; j <= m; j++) {
	u = ((double)j / m)*(n-t+2-.00000001);
	curvepoints[j][0] = curvepoints[j][1] = curvepoints[j][2] = 0.0;
	for (k = 0; k <= n; k++) {
	    temp = N(k, t, u);
	    curvepoints[j][0] += controlpoints[k][0] * temp;
	    curvepoints[j][1] += controlpoints[k][1] * temp;
	    curvepoints[j][2] += controlpoints[k][2] * temp;
	}
    }
}

static double N(k, t, u)
int	k;
int	t;
double	u;
{
    double firstterm, secondterm;

    if (t == 1)
    	if (U[k] <= u && u < U[k+1])
	    return 1.0;
        else
	    return 0.0;
    if (U[k+t-1]-U[k] < 1.0e-10)	/* Test for zero with real #s */
        firstterm = 0.0;
    else
        firstterm = ((u-U[k])/(U[k+t-1]-U[k]))*N(k,t-1,u);
    if (U[k+t]-U[k+1] < 1.0e-10)
        secondterm = 0.0;
    else
        secondterm = ((U[k+t]-u)/(U[k+t]-U[k+1]))*N(k+1,t-1,u);
    return firstterm + secondterm;
}