/*
 * Copyright (c) 1992, The Geometry Center
 *                     University of Minnesota 
 *                     1300 South Second Street
 *                     Minneapolis, MN  55454
 *
 * email address: software@geom.umn.edu
 *
 * This software is copyrighted as noted above.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is
 * preserved on all copies.
 *
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the authors, who may or may not act on them as they desire.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 *
 *     The National Science and Technology Research Center for
 *      Computation and Visualization of Geometric Structures
 */

#include <math.h>

void splineval(double *p0, double *p1,double *p2,double *p3,double *pm,double *dd,double u)
{
    double a,b,c,d,e,f;

    int i;

    for (i=0; i<2; ++i)
	{
	a = p0[i] + (p1[i] - p0[i]) * u; 
	b = p1[i] + (p2[i] - p1[i]) * u; 
	c = p2[i] + (p3[i] - p2[i]) * u; 
	d = a + (b - a) * u;
	e = b + (c - b) * u;
	f = d + (e - d) * u;
	dd[i] = e - d; 	/* store the derivative */
	pm[i] = f;
	}
}

void splineval2(patch, pm, du, dv, u, v)
double patch[16][4], pm[4], du[4], dv[4], u,v;
{
	double p0[4], p1[4], p2[4], p3[4], 
		du0[4], du1[4], du2[4], du3[4], ddu[4];
	/* interpolate in the v-constant curves */
	splineval(patch[0], patch[1], patch[2], patch[3], p0, du0, u);
	splineval(patch[4], patch[5], patch[6], patch[7], p1, du1,  u);
	splineval(patch[8], patch[9], patch[10], patch[11], p2, du2, u);
	splineval(patch[12], patch[13], patch[14], patch[15], p3, du3, u);
	/* interpolate the results */
	splineval(p0, p1, p2, p3, pm, dv, v);
	/* and compute the du vector by interpolation too */
	splineval(du0, du1, du2, du3, du, ddu, v);
}
