/* Subroutine to generate B-spline basis functions and their derivatives for uniform open knot vectors.

	C code for An Introduction to NURBS
	by David F. Rogers. Copyright (C) 2000 David F. Rogers,
	All rights reserved.
	
	Name: dbasisu.c
	Language: C
	Subroutines called: none
	Book reference: Section 3.10

    b1        = first term of the basis function recursion relation
    b2        = second term of the basis function recursion relation
    c         = order of the B-spline basis function
    d1[]      = array containing the derivative of the basis functions
                d1[1]) contains the derivative of the basis function associated with B1 etc.
    d2[]      = array containing the derivative of the basis functions
                d2[1] contains the derivative of the basis function associated with B1 etc.
    f1        = first term of the first derivative of the basis function recursion relation
    f2        = second term of the first derivative of the basis function recursion relation
    f3        = third term of the first derivative of the basis function recursion relation
    f4        = fourth term of the first derivative of the basis function recursion relation
    npts        = number of defining polygon vertices
    n[]      = array containing the basis functions
                n[1]) contains the basis function associated with B1 etc.
    nplusc    = constant -- npts + c -- maximum knot value
    s1        = first term of the second derivative of the basis function recursion relation
    s2        = second term of the second derivative of the basis function recursion relation
    s3        = third term of the second derivative of the basis function recursion relation
    s4        = fourth term of the second derivative of the basis function recursion relation
    t         = parameter value
    temp[]    = temporary array
    x[]       = knot vector
*/	

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

dbasisu(c,t,npts,x,n,d1,d2)

int c,npts;
int x[];
float t;
float n[],d1[],d2[];

{

	int nplusc;
	int i,k;
	float b1,b2;
	float f1,f2,f3,f4;
	float s1,s2,s3,s4;
	float temp[36];		/* allows for 35 defining polygon vertices */
	float temp1[36];
	float temp2[36];

	nplusc = npts + c;

/*    zero the temporary arrays */

	for (i = 1; i <= nplusc; i++){
		temp[i] = 0.;
        temp1[i] = 0.;
        temp2[i] = 0.;
	}

/* calculate the first order basis functions n[i] */

	for (i = 1; i<= nplusc-1; i++){
    	if (( t >= x[i]) && (t < x[i+1]))
			temp[i] = 1;
	    else
			temp[i] = 0;
	}

	if (t == (float)x[npts+1]){		/*    handle the end specially */
		temp[npts] = 1;   		/*resetting the first order basis functions. */
		temp[npts+1]=0;
	}

/* calculate the higher order basis functions */

	for (k = 2; k <= c; k++){
    	for (i = 1; i <= nplusc-k; i++){
        	if (temp[i] != 0)    /* if the lower order basis function is zero skip the calculation */
           		b1 = ((t-x[i])*temp[i])/(x[i+k-1]-x[i]);
	        else
				b1 = 0;

    	    if (temp[i+1] != 0)     /* if the lower order basis function is zero skip the calculation */
        		b2 = ((x[i+k]-t)*temp[i+1])/(x[i+k]-x[i+1]);
	        else
    			b2 = 0;

/*       calculate first derivative */
       
	        if (temp[i] != 0)       /* if the lower order basis function is zero skip the calculation */
	            f1 = temp[i]/(x[i+k-1]-x[i]);
	        else
	            f1 = 0;

        	if (temp[i+1] != 0)     /* if the lower order basis function is zero skip the calculation */
	            f2 = -temp[i+1]/(x[i+k]-x[i+1]);
	        else
	            f2 = 0;

	        if (temp1[i] != 0)      /* if the lower order basis function is zero skip the calculation */
	            f3 = ((t-x[i])*temp1[i])/(x[i+k-1]-x[i]);
	        else
	            f3 = 0;

	        if (temp1[i+1] != 0)    /* if the lower order basis function is zero skip the calculation */
	            f4 = ((x[i+k]-t)*temp1[i+1])/(x[i+k]-x[i+1]);
	        else
	            f4 = 0;

/*       calculate second derivative */

	        if (temp1[i] != 0)      /* if the lower order basis function is zero skip the calculation */
	            s1 = (2*temp1[i])/(x[i+k-1]-x[i]);
	        else
	            s1 = 0;

	        if (temp1[i+1] != 0)      /* if the lower order basis function is zero skip the calculation */
	            s2 = (-2*temp1[i+1])/(x[i+k]-x[i+1]);
	        else 
	            s2 = 0;

	        if (temp2[i] != 0)       /* if the lower order basis function is zero skip the calculation */
	            s3 = ((t-x[i])*temp2[i])/(x[i+k-1]-x[i]);
	        else
	            s3 = 0;

	        	if (temp2[i+1] != 0)    /* if the lower order basis function is zero skip the calculation */
	            s4 = ((x[i+k]-t)*temp2[i+1])/(x[i+k]-x[i+1]);
	        else
	            s4 = 0;

            temp[i] = b1 + b2;
	        temp1[i] = f1 + f2 + f3 + f4;
	        temp2[i] = s1 + s2 + s3 + s4;
		}
	}

/* put in n array	*/

	for (i = 1; i <= npts; i++) {
    	n[i] = temp[i];
	    d1[i] = temp1[i];
    	d2[i] = temp2[i];
	}
}