/*  Subroutine to generate a B-spline curve using an uniform open knot vector

    C code for An Introduction to NURBS
    by David F. Rogers. Copyright (C) 2000 David F. Rogers,
    All rights reserved.
    
    Name: dbspline.c
    Language: C
    Subroutines called: knot.c, basis.c, fmtmul.c
    Book reference: Section 3.5, Ex. 3.4, Alg. p. 288

    b[]        = array containing the defining polygon vertices
                  b[1] contains the x-component of the vertex
                  b[2] contains the y-component of the vertex
                  b[3] contains the z-component of the vertex
    d1[]       = array containing the first derivative of the curve
                  d1[1] contains the x-component
                  d1[2] contains the y-component
                  d1[3] contains the z-component
    d2[]       = array containing the second derivative of the curve
                  d2[1] contains the x-component
                  d2[2] contains the y-component
                  d2[3] contains the z-component
    d1nbasis[]  = array containing the first derivative of the basis functions for a single value of t
    d2nbasis[]  = array containing the second derivative of the basis functions for a single value of t
    k           = order of the \bsp basis function
    nbasis      = array containing the basis functions for a single value of t
    nplusc      = number of knot values
    npts        = number of defining polygon vertices
    p[]         = array containing the curve points
                  p[1] contains the x-component of the point
                  p[2] contains the y-component of the point
                  p[3] contains the z-component of the point
    p1          = number of points to be calculated on the curve
    t           = parameter value 0 <= t <= 1
    x[]         = array containing the knot vector
*/

dbspline(npts,k,p1,b,p,d1,d2)

int npts,k,p1;

float b[];
float p[];
float d1[];
float d2[];

{
    int i,j,icount,jcount;
    int i1;
    int x[30];        /* allows for 20 data points with basis function of order 5 */
    int nplusc;

    float step;
    float t;
    float nbasis[20];
    float d1nbasis[20];
    float d2nbasis[20];
    float temp;


    nplusc = npts + k;

/*  zero and redimension the knot vector and the basis arrays */

    for(i = 0; i <= npts; i++){
         nbasis[i] = 0.;
         d1nbasis[i]=0.;
         d2nbasis[i]=0.;
    }

    for(i = 0; i <= nplusc; i++){
         x[i] = 0.;
        }

/* generate the uniform open knot vector */

    knotu(npts,k,x);

/*
    printf("The knot vector is ");
    for (i = 1; i <= nplusc; i++){
        printf(" %d ", x[i]);
    }
    printf("\n");
*/

    icount = 0;

/*    calculate the points on the bspline curve */

    t = k-1;
    step = ((float)(npts-(k-1)))/((float)(p1-1));

    for (i1 = 1; i1<= p1; i1++){

        if ((float)npts - t < 5e-6){
            t = (float)npts;
        }

        dbasisu(k,t,npts,x,nbasis,d1nbasis,d2nbasis);      /* generate the basis function for this value of t */
/*
        printf("t = %f \n",t);
        printf("nbasis = ");
        for (i = 1; i <= npts; i++){
            printf("%f  ",nbasis[i]);
        }
        printf("\n");
*/
/*
        printf("t = %f \n",t);
        printf("d1nbasis = ");
        for (i = 1; i <= npts; i++){
            printf("%f  ",d1nbasis[i]);
        }
        printf("\n");
*/
/*
        printf("t = %f \n",t);
        printf("d2nbasis = ");
        for (i = 1; i <= npts; i++){
            printf("%f  ",d2nbasis[i]);
        }
        printf("\n");
*/

        for (j = 1; j <= 3; j++){      /* generate a point on the curve */
            jcount = j;
            p[icount+j] = 0.;
            d1[icount+j] = 0.;
            d2[icount+j] = 0.;

            for (i = 1; i <= npts; i++){ /* Do local matrix multiplication */
                p[icount + j] = p[icount + j] + nbasis[i]*b[jcount];
                d1[icount + j] = d1[icount + j] + d1nbasis[i]*b[jcount];
                d2[icount + j] = d2[icount + j] + d2nbasis[i]*b[jcount];                
/*
                printf("jcount,nbasis,b,nbasis*b,p = %d %f %f %f %f\n",jcount,nbasis[i],b[jcount],temp,p[icount+j]);
*/
                jcount = jcount + 3;
            }
        }

/*        printf("icount, p %d %f %f %f \n",icount,p[icount+1],p[icount+2],p[icount+3]);*/

        icount = icount + 3;
        t = t + step;
    }
}