/*
	Test program for C code from An Introduction to NURBS
	by David F. Rogers. Copyright (C) 2000 David F. Rogers,
	All rights reserved.
	
	Name: trbsurf.c
	Purpose: Test the fast rational B-spline surface generator.

	Language: C
	Subroutines called: rbsurf.c sumrbas.c knot.c basis.c
	Book reference:

*/

#include <stdio.h>

/* #include "rbsurf.h"*/

main(){

    int i,i1,i2;
    int npts,mpts;
    int k,l;
    int p1,p2;
    int itest;
    int ibnum;
    int plusminus;

    int ch;
    
    static float b[81];
    static float q[126];
    static float niku[126];
    static float mjlw[126];
    static float rsumij[126];
    static float savrsumij[126];
    static float bold[5];

    float change;

    int stack_probe;
/*  fprintf(stderr, "main:\t\t&stack_probe = %lu\n", (long)&stack_probe); */
    fflush(stderr);

/*  set net and surface parameters */

    npts = 4;     /* use npts = 5 for asymmetrical testing -- see below */
    mpts = 4;
    k = 3;
    l = 3;

    p1 = 5;
    p2 = 5;

    itest = 0;    /* set to zero to get full surface calculation on 1st pass */

/*  printf("k,l,npts,mpts,p1,p2 = %d %d %d %d %d %d \n",k,l,npts,mpts,p1,p2);*/

/*  zero arrays */

for (i = 1; i <= 4*npts; i++){
        b[i] = 0.;
    }

    for (i = 1; i <= 3*p1*p2; i++){
        q[i] = 0.;
    }

    for (i = 1; i <= npts*p1; i++){
        niku[i] = 0.;
    }

    for (i = 1; i <= mpts*p2; i++){
        mjlw[i] = 0.;
    }

    for (i = 1; i <= p1*p2; i++){
        rsumij[i] = 1.;
        savrsumij[i] = 1.;
    }

/*  define polygon net */

    b[1] = -15.;
    b[2] = 0.;
    b[3] = 15.;
    b[4] = 1;
    b[5] = -15.;
    b[6] = 5.;
    b[7] = 5.;
    b[8] = 1;
    b[9] = -15.;
    b[10] = 5.;
    b[11] = -5.;
    b[12] = 1;
    b[13] = -15.;
    b[14] = 0.;
    b[15] = -15.;
    b[16] = 1;

    b[17] = -5.;
    b[18] = 5.;
    b[19] = 15.;
    b[20] = 1;
    b[21] = -5.;
    b[22] = 10.;
    b[23] = 5.;
    b[24] = 1;
    b[25] = -5.;
    b[26] = 10.;
    b[27] = -5.;
    b[28] = 1;
    b[29] = -5.;
    b[30] = 5.;
    b[31] = -15.;
    b[32] = 1;

    b[33] = 5.;
    b[34] = 5.;
    b[35] = 15.;
    b[36] = 1;
    b[37] = 5.;
    b[38] = 10.;
    b[39] = 5.;
    b[40] = 1;
    b[41] = 5.;
    b[42] = 10.;
    b[43] = -5.;
    b[44] = 1;
    b[45] = 5.;
    b[46] = 0.;
    b[47] = -15.;
    b[48] = 1;

    b[49] = 15.;
    b[50] = 0.;
    b[51] = 15.;
    b[52] = 1;
    b[53] = 15.;
    b[54] = 5.;
    b[55] = 5.;
    b[56] = 1;
    b[57] = 15.;
    b[58] = 5.;
    b[59] = -5.;
    b[60] = 1;
    b[61] = 15.;
    b[62] = 0.;
    b[63] = -15.;
    b[64] = 1;

/*  change n1 to 5 and uncomment the b[] assignments below
    to get a 5 x 4 polygon for asymmetrical testing */
/*

    b[65] = 15.;
    b[66] = 0.;
    b[67] = 15.;
    b[68] = 1;
    b[69] = 15.;
    b[70] = -5.;
    b[71] = 5.;
    b[72] = 1;
    b[73] = 15.;
    b[74] = -5.;
    b[75] = -5.;
    b[76] = 1;
    b[77] = 15.;
    b[78] = -5.;
    b[79] = -15.;
    b[80] = 1;
*/
    printf("\nPress Enter to start the timing. You need to use a stop watch.");
    ch=getchar();

/*  To generate single surface uncomment the next line and comment out
    all the statements below the next except the last } labelled end main.
*/

/*  rbsurf(b,k,l,npts,mpts,p1,p2,&itest,ibnum,bold,niku,mjlw,rsumij,savrsumij,q); */

    printf("start\n\n");

/*  generate multiple surfaces for checking incremental calculation and timing */

    plusminus = 1;
    change = 4;
    for (i = 1; i<=30000; i++){ 
/*        printf("i = %d \n",i); */
        rbsurf(b,k,l,npts,mpts,p1,p2,&itest,ibnum,bold,niku,mjlw,rsumij,savrsumij,q);
           ibnum = 21;
           if (i == 1.){              /* assign bold after first pass */
               bold[1] = b[ibnum];
               bold[2] = b[ibnum+1];
               bold[3] = b[ibnum+2];
               bold[4] = b[ibnum+3];
/*             printf("i = %d and bold = %f %f %f \n",i,bold[ibnum],bold[ibnum+1],bold[ibnum+2]); */
           } /* end if */

/*  change the 4 in the for loop above to 10000 or more
    and uncomment the next two lines for timing.
*/

/*
    Uncomment the following two statements to get a valid timing loop
    for a change in the homogeneous weighting factor.
*/

    b[ibnum+3] = b[ibnum+3] + plusminus*change;
    plusminus = -plusminus;
/*    printf("h= %f \n",b[ibnum+3]); */

/*
    Uncomment the following two statements to get a valid timing loop
    for a change in a polygon net vertex.
*/

/*

/*
    b[ibnum] = - b[ibnum];
*/

    } /* end of main for loop */

    printf("\ndone");
    printf("\nPress Enter to list the results of the last calculation.");
    ch=getchar();

/*  print only the last result */

    printf("\nPolygon points\n\n");

    for (i1 = 1; i1 <= 4*npts*mpts; i1=i1+4){
        printf(" %f %f %f %f\n",b[i1],b[i1+1],b[i1+2],b[i1+3]);
    }

    printf("\nPress Enter to see the surface points.");
    ch=getchar();

    printf("\nSurface points\n\n");

    for (i2 = 1; i2 <= 3*p1*p2; i2=i2+3){
        printf("%f, %f, %f \n",q[i2],q[i2+1],q[i2+2]);
    }

} /* end main */