/*

                              DISCLAIMER
                              ==========

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.   

    If the software is modified by someone else and passed on, we, the authors
    want its recipients to know that what they have is not the original, so
    that any problems introduced by others will not reflect on the original
    authors' reputations.
*/                                            

/*---------------------------------------------------------------

	Module Name: Nrb_Procs

	Module Description: Contains the code for non-uniform rational
			    B-Splines

	Authors: W T Hewitt, M Patel, S Larkin

	Last mod date:  2/8/88   Creation date

	Version: 1

	Changes:

		31/01/90Steve Larkin

				Module has been ported to the VaxStation 3540

		20/02/90Steve Larkin

				Appended some routines onto the end.
				These include :-



 ---------------------------------------------------------------*/
#include "nurbh.h"
#include <stdlib.h>
#include <malloc.h>

/*---------------------------------------------------------------------------*/

static Pint        nrb_maxmess = 10;
static Pint        nrb_messcount = 0;
static PR_nurb*    nrb_nrblist= NULL;
static PR_knots*   nrb_knotlist=NULL;
static PR_pts*     nrb_ptslist=NULL;
static PR_stats    nrb_nrbstats={0};
static PR_stats    nrb_ptstats={0};
static PR_stats    nrb_knotstats={0};
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_initialise(C(void))
PreANSI(void)
/*
Initialise the NURB library.
*/
{
  /*Initialise freelist pointers*/
  nrb_nrblist = NULL;
  nrb_knotlist = NULL;
  nrb_ptslist = NULL;

  memset(&nrb_nrbstats, 0, sizeof(PR_stats));
  memset(&nrb_ptstats, 0, sizeof(PR_stats));
  memset(&nrb_knotstats, 0, sizeof(PR_stats));

  nrb_ptstats.pf_minBsize = LONG_MAX;
  nrb_knotstats.pf_minBsize = -LONG_MAX;
  nrb_maxmess = 10;
  nrb_messcount = 0;
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_error(C(char *)str)
PreANSI(char *str)
/*
Output procedure library Error Messages
*/
{
  if (nrb_messcount < nrb_maxmess) {
    fprintf(stderr,"NURB - %s\n", str);
  }

  nrb_messcount++;

  if (nrb_messcount == nrb_maxmess) {
    fprintf(stderr,"%s\n",
     "NRB_Error --- WARNING further error messages will be suppressed");
  }

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_seterrorcount(C(Pint) i)
PreANSI(Pint i)
/*
Set the maximun number of errors that are reported. The default is 10.
*/
{
  nrb_maxmess = P_max(0, i);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern Pint nrb_inqerrorcount(C(Pint *)i)
PreANSI(Pint *i)
/*
Returns the maxium number of errors that are reported. 
*/
{
  *i = nrb_maxmess;
  return (Pint) nrb_maxmess;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern PR_knots* nrb_deallocate(C(PR_knots *)ptr, C(PR_knots *)listhead, 
                                C(PR_stats *)stats)
PreANSI(PR_knots *ptr)
PreANSI(PR_knots *listhead) 
PreANSI(PR_stats *stats)
/*
Deallocate a set of knots.
*/
{
  boolean found;
  PR_knots *tptr, *prev;

  /*Deallocate knots array by adding it to knots freelist*/
  /*List stored in sorted order! - largest first*/

  prev = NULL;
  tptr = listhead;
  found = FALSE;

  stats->pf_dreqs++;

  if (ptr == NULL)
    return (PR_knots*) listhead;

  stats->pf_inuse--;
  stats->pf_Qlen++;
  stats->pf_maxQlen = P_max(stats->pf_maxQlen, stats->pf_Qlen);
  if (stats->pf_minBsize == 0)
    stats->pf_minBsize = LONG_MAX;
  /*Not been initialised minsize must be >0!*/
  stats->pf_minBsize = P_min(stats->pf_minBsize, ptr->pf_size);
  stats->pf_maxBsize = P_max(stats->pf_maxBsize, ptr->pf_size);
  while (tptr != NULL && !found) {
    if (tptr->pf_size <= ptr->pf_size)
      found = TRUE;
    else {  /*Skip to next element*/
      prev = tptr;
      tptr = tptr->pf_next;
    }
  }

  if (!found) {
    if (prev == NULL) {   /*Add as listhead*/
      listhead = ptr;
      listhead->pf_next = NULL;
    } else {  /*Add at enD of list*/
      prev->pf_next = ptr;
      ptr->pf_next = NULL;
    }
  } else  /*Found*/
  {  /*Link into list.*/
    if (prev == NULL) {   /*Add in front of first element of list*/
      ptr->pf_next = listhead;
      listhead = ptr;
    } else {
      ptr->pf_next = prev->pf_next;
      prev->pf_next = ptr->pf_next;
    }
  }

  return (PR_knots *)listhead;
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern PR_knots* nrb_allocate(C(Pint) ntx, C(PR_knots *)ptr, 
                              C(PR_knots **)listhead, 
                              C(PR_stats *)stats)
PreANSI(Pint ntx)
PreANSI(PR_knots *ptr)
PreANSI(PR_knots **listhead) 
PreANSI(PR_stats *stats)
/*
Allocate a set of knots.
*/
{
  boolean found;
  Pint isize;

  PR_knots *prev;

  /*if not nil and size same and return
    otherwise deallocate, allocate correct size and return*/

  /*printf("NRB_allocate, size:%4ld\n",ntx);*/
  
  stats->pf_areqs++;

  if (ptr != NULL) {
    if (ntx != ptr->pf_size){
      *listhead = nrb_deallocate(ptr, *listhead, stats);
      ptr = NULL;
    }
  }

  /*Here ptr = NIL or is the correct size*/

  if (ptr != NULL) {
    stats->pf_nochange++;
    /*printf("NRB_allocate, ptr:%10ld\n",(long) ptr);*/
    return (PR_knots *)ptr;
  }

  stats->pf_inuse++;
  found = FALSE;   /*Search freelist for one of correct size */
  ptr = *listhead;
  prev = NULL;


  while ((ptr != NULL) && (!found) ) {  
    if (ptr->pf_size == ntx)
      found = TRUE;
    else {  /*Skip to next element*/
      prev = ptr;
      ptr = ptr->pf_next;
    }
  }

  if (!found) {   /*Real allocate!*/

    isize = 2 * sizeof(Pint) + ntx * sizeof(Pfloat);
    ptr = NULL;
    ptr = (PR_knots *)malloc((size_t)isize);
    if (ptr == NULL) nrb_error("NRB_allocate --- malloc failed");
    ptr->pf_size = ntx;
    ptr->pf_next = NULL;
    stats->pf_news++;
    /*printf("NRB_allocate, ptr:%10ld\n",(long) ptr);*/
    return (PR_knots *)ptr;
  }
  stats->pf_Qlen--;
  /*UnLink from List. prev = nil means adjust listhead
		     otherwise adjust pointers*/
  if (prev == NULL)
    *listhead = ptr->pf_next;
  else
    prev->pf_next = ptr->pf_next;

/*printf("NRB_allocate, ptr:%10ld\n",(long) ptr);*/
return (PR_knots *)ptr;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_deallocateknots(C(PR_dir *)kk)
PreANSI(PR_dir *kk)
/*
Deallocate a set of knots.
*/
{
  nrb_knotlist =nrb_deallocate(kk->pf_kk, nrb_knotlist, &nrb_knotstats);
  kk->pf_kk = NULL;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_allocateknots(C(PR_dir *)kk)
PreANSI(PR_dir *kk)
/*
Allocate a set of knots.
*/
{
  kk->pf_kk = nrb_allocate(kk->pf_nt, kk->pf_kk, &nrb_knotlist, &nrb_knotstats);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern  void nrb_deallocatepts(C(PR_pts **)ptr)
PreANSI(PR_pts **ptr)
/*
Deallocate points array by adding it to pts freelist
*/
{
nrb_ptslist = (PR_pts*)nrb_deallocate((PR_knots*) *ptr, (PR_knots*)nrb_ptslist,&nrb_ptstats);
*ptr = NULL;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_allocatepts(C(Pint) nx, C(PR_pts **)ptr)
PreANSI(Pint nx) 
PreANSI(PR_pts **ptr)
/*
Allocate a set of points.
*/
{


PR_knots *hack;
hack=(PR_knots *) nrb_ptslist;

*ptr = (PR_pts*) nrb_allocate(nx*4,(PR_knots*)*ptr,&hack,&nrb_ptstats);
nrb_ptslist= (PR_pts*) hack;

/****ORIG
*ptr = (PR_pts*) nrb_allocate(nx*4,(PR_knots*)*ptr,&(PR_knots*)nrb_ptslist,&nrb_ptstats);
*******/

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_clear(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Initialise a NURB to NULL.
*/
{
  nrb->pf_u.pf_kk = NULL;
  nrb->pf_v.pf_kk = NULL;
  nrb->pf_ppp = NULL;
  nrb->pf_next = NULL;
  nrb->pf_prev = NULL;
  nrb->pf_id = 0;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_allocatenurb(C(PR_nurb **)ptr)
PreANSI(PR_nurb **ptr)
/*
Allocate enough space for a pointer to a NURB.
*/
{
  Pint isize;

  nrb_nrbstats.pf_areqs++;
  nrb_nrbstats.pf_inuse++;

  if (nrb_nrblist == NULL) {   
    isize = sizeof(PR_nurb);
    *ptr = (PR_nurb *) malloc(isize);
    if ((*ptr) == NULL) nrb_error("NRB_allocatenurb --- malloc failed");
    (*ptr)->pf_next = NULL;
    (*ptr)->pf_prev = NULL;    
    nrb_nrbstats.pf_news++;
  } else {
    nrb_nrbstats.pf_Qlen--;
    *ptr = nrb_nrblist;
    nrb_nrblist = (*ptr)->pf_next;
  }

  if (*ptr != NULL) nrb_clear(*ptr);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_deallocatenurb(C(PR_nurb **)ptr)
PreANSI(PR_nurb **ptr)
/*
Deallocate a NURB and place it on the freelist.
*/
{
  /*Deallocate points and knots*/

  nrb_deallocateknots(&(*ptr)->pf_u);
  nrb_deallocateknots(&(*ptr)->pf_v);
  nrb_deallocatepts(&(*ptr)->pf_ppp);

  /*Deallocate Nurb - by adding it to the free list*/

  (*ptr)->pf_next = nrb_nrblist;
  nrb_nrblist = *ptr;

  nrb_nrbstats.pf_inuse--;
  nrb_nrbstats.pf_dreqs++;
  nrb_nrbstats.pf_Qlen++;
  nrb_nrbstats.pf_maxQlen = P_max(nrb_nrbstats.pf_maxQlen,
				   nrb_nrbstats.pf_Qlen);
  *ptr=NULL;			     
}
/*----------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_stats(C(FILE *)fp)
PreANSI(FILE *fp)
/*
Print out the overall library statistics to the file pointed to by fp.
*/
{
  fprintf(fp,"%s\n","***************************");
  fprintf(fp,"%s\n\n%s\n","NURB - free list statistics","NURBs - PR_NURB");

  fprintf(fp,"%s%8ld\n","  Currently inuse: ", nrb_nrbstats.pf_inuse);
  fprintf(fp,"%s%4ld\n","  Deallocate requests: ", nrb_nrbstats.pf_dreqs);
  fprintf(fp,"%s%9ld\n","     CurrentQlen: ", nrb_nrbstats.pf_Qlen);
  fprintf(fp,"%s%13ld\n","     MaxQlen: ", nrb_nrbstats.pf_maxQlen);
  fprintf(fp,"%s%7ld\n","  Allocate requests:", nrb_nrbstats.pf_areqs);
  fprintf(fp,"%s%18ld\n","     new:", nrb_nrbstats.pf_news);

  fprintf(fp,"\n%s\n","Knots - PR_knots");

  fprintf(fp,"%s%8ld\n","  Currently inuse: ", nrb_knotstats.pf_inuse);
  fprintf(fp,"%s%4ld\n","  Deallocate requests: ", nrb_knotstats.pf_dreqs);
  fprintf(fp,"%s%9ld\n","     CurrentQlen: ", nrb_knotstats.pf_Qlen);
  fprintf(fp,"%s%13ld\n","     MaxQlen: ", nrb_knotstats.pf_maxQlen);
  fprintf(fp,"%s%10ld\n","     MinBlkSize: ", nrb_knotstats.pf_minBsize);
  fprintf(fp,"%s%10ld\n","     MaxBlkSize: ", nrb_knotstats.pf_maxBsize);
  fprintf(fp,"%s%7ld\n","  Allocate requests:", nrb_knotstats.pf_areqs);
  fprintf(fp,"%s%13ld\n","     Nochange:", nrb_knotstats.pf_nochange);
  fprintf(fp,"%s%18ld\n","     new:", nrb_knotstats.pf_news);

  fprintf(fp,"\n%s\n","Points - PR_points");
  fprintf(fp,"%s%8ld\n","  Currently inuse: ", nrb_ptstats.pf_inuse);
  fprintf(fp,"%s%4ld\n","  Deallocate requests: ", nrb_ptstats.pf_dreqs);
  fprintf(fp,"%s%9ld\n","     CurrentQlen: ", nrb_ptstats.pf_Qlen);
  fprintf(fp,"%s%13ld\n","     MaxQlen: ", nrb_ptstats.pf_maxQlen);
  fprintf(fp,"%s%10ld\n","     MinBlkSize: ", nrb_ptstats.pf_minBsize);
  fprintf(fp,"%s%10ld\n","     MaxBlkSize: ", nrb_ptstats.pf_maxBsize);
  fprintf(fp,"%s%7ld\n","  Allocate requests:", nrb_ptstats.pf_areqs);
  fprintf(fp,"%s%13ld\n","     Nochange:", nrb_ptstats.pf_nochange);
  fprintf(fp,"%s%18ld\n","     new:", nrb_ptstats.pf_news);

  fprintf(fp,"\n%s\n","***************************");
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_init(C(Pint) kx, C(Pint) Ntx, C(Pint) ky, C(Pint) nty, 
                     C(PR_nurb *)nrb)
PreANSI(Pint kx) 
PreANSI(Pint Ntx) 
PreANSI(Pint ky) 
PreANSI(Pint nty) 
PreANSI(PR_nurb *nrb)
/*
Initialise a 'empty' NURB, given the number of knots and degree
*/
{
  /*nrb_init*/
  nrb->pf_u.pf_k = kx;
  nrb->pf_u.pf_nt = Ntx;
  nrb->pf_u.pf_n = Ntx - kx;

  nrb->pf_v.pf_k = P_max(0, ky);
  nrb->pf_v.pf_nt = P_max(1, nty);
  nrb->pf_v.pf_n = P_max(1, nty - ky);

  nrb_allocateknots(&nrb->pf_u);   /*allocate knot array*/
  nrb_allocateknots(&nrb->pf_v);   /*allocate knot array*/

  nrb_allocatepts(nrb->pf_u.pf_n * nrb->pf_v.pf_n, &nrb->pf_ppp);
      /*allocate control point array*/
  
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_copyknots(C(PR_dir *)k1, C(PR_dir *)k2)
PreANSI(PR_dir *k1) 
PreANSI(PR_dir *k2)
/*
Copy a set of knots.
*/
{
  Pint i, mm;

  k2->pf_nt = k1->pf_nt;
  k2->pf_k = k1->pf_k;
  k2->pf_n = k1->pf_n;
  nrb_allocateknots(k2);
  mm = k1->pf_nt;
  for (i = 0; i < mm; i++)
    k2->pf_kk->knots[i] = k1->pf_kk->knots[i];
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_copy(C(PR_nurb *)nrb1, C(PR_nurb *)nrb2)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb2)
/*
Copy a NURB to another NURB.
*/
{
  Pint i, mm;

  *nrb2 = *nrb1;   /*Copy all the fields, overwriting the array addresses*/
  nrb2->pf_ppp = NULL;
  nrb2->pf_u.pf_kk = NULL;
  nrb2->pf_v.pf_kk = NULL;
  nrb2->pf_next = NULL;
  nrb2->pf_prev = NULL;

  nrb_copyknots(&nrb1->pf_u, &nrb2->pf_u);
  nrb_copyknots(&nrb1->pf_v, &nrb2->pf_v);

  nrb_allocatepts(nrb2->pf_u.pf_n * nrb2->pf_v.pf_n, &nrb2->pf_ppp);

  mm = nrb1->pf_u.pf_n * nrb1->pf_v.pf_n;
  /*Now copy the arrays*/
  for (i = 0; i < mm; i++)   
    nrb2->pf_ppp->pts[i] = nrb1->pf_ppp->pts[i];

}
/*---------------------------------------------------------------------------*/
void swap(C(Pint *)a, C(Pint *)b)
PreANSI(Pint *a)
PreANSI(Pint *b)
{
  Pint temp;

  temp = *a;
  *a = *b;
  *b = temp;
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_interchange(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Put the #knots and ctrl pts back
Most of the NURB Surface algorithms `transpose' x and Y. This
    is used to correct such effects when the algorithm is called
    upon a curve
*/
{
  /*nrb_interchange*/
  swap(&nrb->pf_u.pf_n, &nrb->pf_v.pf_n);   /*interchange # ctrl pts*/
  swap(&nrb->pf_u.pf_nt, &nrb->pf_v.pf_nt);   /*Interchange # knots*/
  swap(&nrb->pf_u.pf_k, &nrb->pf_v.pf_k);   /*interchange order*/
  
  return;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern boolean nrb_getpoint4(C(FILE *)fp, C(Ppoint4 *)pt)
PreANSI(FILE *fp) 
PreANSI(Ppoint4 *pt)
/*
Read a homogeneous point from file;
    Returns TRUE if successful
*/
{

   if(fscanf(fp,"%f %f %f %f", &pt->pfhx,&pt->pfhy,&pt->pfhz,&pt->pfhw)!=4)
      return(FALSE);
   else
      return(TRUE);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_putpoint4(C(FILE *)fp, C(Ppoint4 *)pt)
PreANSI(FILE *fp) 
PreANSI(Ppoint4 *pt)
/*
Write out a 4D point to file fd, in ff format, terminated by newline
*/
{
  fprintf(fp,"%16.6g%16.6g%16.6g%16.6g\n",
		      pt->pfhx,pt->pfhy,pt->pfhz,pt->pfhw);
  return;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dumpl(C(FILE *)fp, C(char *)str, C(PR_knots *)listhead)
PreANSI(FILE *fp) 
PreANSI(char *str) 
PreANSI(PR_knots *listhead)
/* 
Dump the points or knots list
*/
{
PR_knots *ptr;


  fprintf(fp,"------------------------\n");
  fprintf(fp,"%s\n",str);

   if (listhead == NULL){
    fprintf(fp,"    This list is empty\n");
    fprintf(fp,"------------------------\n");
    return;
   }

  ptr = listhead;

  fprintf(fp,"Listhead: %d\n",ptr);
  while (ptr != NULL){
   fprintf(fp,"          *pf_next:%d\n",ptr->pf_next);
   fprintf(fp,"           pf_size:%d bytes\n",ptr->pf_size);
   fprintf(fp,"      *knots/pts[]:%d\n",ptr->knots);
   fprintf(fp,"      ********************\n");
   ptr = ptr ->pf_next;
  } 

fprintf(fp,"------------------------\n");
return;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dumplists(C(FILE *)fp)
PreANSI(FILE *fp)
/*
Dump the list of points or knots.
*/
{
PR_nurb *ptr;

  nrb_dumpl(fp,"Dump of points list",(PR_knots *) nrb_ptslist);
  nrb_dumpl(fp,"Dump of knots list",nrb_knotlist);

  fprintf(fp,"------------------------\n");
  fprintf(fp,"Dump of nurb list\n");

   if (nrb_nrblist == NULL){
    fprintf(fp,"    This list is empty\n");
    fprintf(fp,"------------------------\n");
    return;
   }

  ptr = nrb_nrblist;

  fprintf(fp,"Listhead: %d\n",ptr);
  while (ptr != NULL){
   fprintf(fp,"    *pf_next:%d\n",ptr->pf_next);
   fprintf(fp,"    *****************\n");
   ptr = ptr ->pf_next;
  } 
return;
}/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dumpknots(C(FILE *)fp, C(char *)str, C(PR_dir *)kk)
PreANSI(FILE *fp) 
PreANSI(char *str) 
PreANSI(PR_dir *kk)
/*
Dump the list of knots.
*/
{
  Pint i, mm;

  fprintf(fp,"%s%s%s %4d","#Knots(",str,"):", kk->pf_nt);

  mm = kk->pf_nt;
  for (i = 0; i < mm; i++) {
    if ((i % 4) == 0) fprintf(fp,"\n");
    fprintf(fp,"%16.6g",kk->pf_kk->knots[i]);
  }
  fprintf(fp,"\n");
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dumppts(C(FILE *)fp, C(Pint) nx, C(Pint) ny, C(Ppoint4 *)pts, 
                        C(char *)str)
PreANSI(FILE *fp) 
PreANSI(Pint nx) 
PreANSI(Pint ny)
PreANSI(Ppoint4 *pts) 
PreANSI(char *str)
/*
Dump the list of points.
*/
{
  /* Points*/
  Pint i, j,jj;

  fprintf(fp,"%s%s%s %4ld %4ld\n", "#",str,"(X,Y):",nx,ny);
  fprintf(fp,"%s\n",
    " x  y            pfhx            pfhy            pfhz            pfhw");
  for (j = 1; j <= ny; j++) {
    for (i = 1; i <= nx; i++) {
      jj = ni(i,j,nx);
      fprintf(fp,"%2ld %2ld", i, j);
      nrb_putpoint4(fp, &pts[jj - 1]);
    }
  }
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dumpsummary(C(FILE *)fp, C(PR_nurb *)nrb)
PreANSI(FILE *fp) 
PreANSI(PR_nurb *nrb)
/*
 Write out a summary of a NURB Surface to file fd
*/
{
  fprintf(fp,"%s%4ld%4ld\n", "Order(X,Y):",nrb->pf_u.pf_k,nrb->pf_v.pf_k);

  nrb_dumpknots(fp, "X", &nrb->pf_u);
  nrb_dumpknots(fp, "Y", &nrb->pf_v);

  fprintf(fp,"%s %4ld %4ld\n ", "#Control points(X,Y):",nrb->pf_u.pf_n,
		   nrb->pf_v.pf_n, 2);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_dump(C(FILE *)fp, C(PR_nurb *)nrb)
PreANSI(FILE *fp) 
PreANSI(PR_nurb *nrb)
/*
Write out a NURB Surface to file fd.
*/
{
  fprintf(fp,"%s %4ld %4ld\n", "Order(X,Y):",nrb->pf_u.pf_k,nrb->pf_v.pf_k);

  nrb_dumpknots(fp, "X", &nrb->pf_u);
  nrb_dumpknots(fp, "Y", &nrb->pf_v);

  nrb_dumppts(fp, nrb->pf_u.pf_n, nrb->pf_v.pf_n,
		  nrb->pf_ppp->pts, "Control points");
}  /* nrb_dump*/

/*----------------------------------------------------------------------------*/

#define doerror(A)  {eflag =FALSE; nrb_error(A);}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern boolean nrb_read(C(FILE *)fp, C(PR_nurb *)nrb)
PreANSI(FILE *fp) 
PreANSI(PR_nurb *nrb)
/*
returns TRUE if it can sucessfuly read the surface
Reads
    X order
    Y order
    X number of knots (ntx)
    1..ntx knots
    Y number of knots (nty)
    1..nty knots
    nx,ny (numbers of control points)
    nx*ny 4D points
*/
{
  boolean eflag;
  Pint i, mm;

  eflag = TRUE;   /*Assume no problems*/

  if (fscanf(fp,"%d %d",&nrb->pf_u.pf_k,&nrb->pf_v.pf_k) != 2 ) {
    doerror("NRB_READ : expected pf_u.pf_k");
  }
   printf("TWO numbers %d %d\n",nrb->pf_u.pf_k,nrb->pf_v.pf_k);

  if (fscanf(fp,"%d", &nrb->pf_u.pf_nt) != 1) {
    doerror("NRB_READ : expected pf_u.pf_nt");
  }
   printf("NUmber of U knots is %d\n",nrb->pf_u.pf_nt);

  nrb_allocateknots(&nrb->pf_u);
  mm = nrb->pf_u.pf_nt;
  for (i = 0; i < mm; i++) 
   {
      fscanf(fp,"%g",&nrb->pf_u.pf_kk->knots[i]);
      printf("%f\n",nrb->pf_u.pf_kk->knots[i]);   
   }

  if (fscanf(fp,"%ld", &nrb->pf_v.pf_nt) != 1) 
  {
    doerror("NRB_READ : expected pf_v.pf_nt" );
  }
  printf("NUmber of V knots is %d\n",nrb->pf_v.pf_nt);

  nrb_allocateknots(&nrb->pf_v);

  mm = nrb->pf_v.pf_nt;
  for (i = 0; i < mm; i++) {
    if (1 != fscanf(fp,"%g",&nrb->pf_v.pf_kk->knots[i])) {
      doerror("NRB_READ : expected pf_ty" );
    }
  }

  if (1 != fscanf(fp,"%d", &nrb->pf_u.pf_n)) {
    doerror("NRB_READ : expected pf_u.pf_n" );
  }
  if (1 != fscanf(fp,"%d", &nrb->pf_v.pf_n)) {
    doerror("NRB_READ : expected pf_v.pf_n" );
  }

  nrb_allocatepts(nrb->pf_u.pf_n * nrb->pf_v.pf_n, &nrb->pf_ppp);

  mm = nrb->pf_u.pf_n * nrb->pf_v.pf_n;
  for (i = 0; i < mm; i++) {   
    if (!nrb_getpoint4(fp, &nrb->pf_ppp->pts[i]))
      doerror("NRB_READ : expected PRPoint4 pf_ppp^.pts" );
  }

  return (boolean) eflag;
}  /* nrb_read */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_write(C(FILE *)fp, C(PR_nurb *)nrb)
PreANSI(FILE *fp) 
PreANSI(PR_nurb *nrb)
/*
Write a NURB to the file pointed to by fp.
*/
{
  Pint i, mm;

  fprintf(fp,"%4ld %4ld\n", nrb->pf_u.pf_k,nrb->pf_v.pf_k);
  fprintf(fp,"%4ld\n", nrb->pf_u.pf_nt);

  mm = nrb->pf_u.pf_nt;
  for (i = 0; i < mm; i++) {
    fprintf(fp,"%16.6g\n",nrb->pf_u.pf_kk->knots[i]);
  }

  fprintf(fp,"%4ld\n", nrb->pf_v.pf_nt);

  mm = nrb->pf_v.pf_nt;
  for (i = 0; i < mm; i++) {
    fprintf(fp,"%16.6g\n",nrb->pf_v.pf_kk->knots[i]);
  }

  fprintf(fp,"%4ld %4ld\n", nrb->pf_u.pf_n,nrb->pf_v.pf_n);

  mm = nrb->pf_u.pf_n * nrb->pf_v.pf_n;
  for (i = 0; i < mm; i++)   
    nrb_putpoint4(fp, &nrb->pf_ppp->pts[i]);

}  /* nrb_write */



/*+----------------------------------------------------------------------------*/
/*function:external*/
extern Ppoint4 tp_litv4(C(Pfloat) x, C(Pfloat) y, C(Pfloat) z, C(Pfloat) w)
PreANSI(Pfloat x)
PreANSI(Pfloat y)
PreANSI(Pfloat z)
PreANSI(Pfloat w)
/*
Convert 4 reals to a point
*/
{
  Ppoint4 temp;

  temp.pfhx = x;
  temp.pfhy = y;
  temp.pfhz = z;
  temp.pfhw = w;
  return temp;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern Ppoint4 tp_scale4(C(Ppoint4) s1, C(Pfloat) v)
PreANSI(Ppoint4 s1) 
PreANSI(Pfloat v)
/*
Scale a 4D vector by v and return the result
*/
{
  Ppoint4 temp;

  temp.pfhx = s1.pfhx * v;
  temp.pfhy = s1.pfhy * v;
  temp.pfhz = s1.pfhz * v;
  temp.pfhw = s1.pfhw * v;
  return temp;

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern Ppoint4 tp_add4(C(Ppoint4) s1, C(Ppoint4) s2)
PreANSI(Ppoint4 s1) 
PreANSI(Ppoint4 s2)
/*
Add two 4D vectors, and return the result
*/
{
  Ppoint4 temp;

  temp.pfhx = s1.pfhx + s2.pfhx;
  temp.pfhy = s1.pfhy + s2.pfhy;
  temp.pfhz = s1.pfhz + s2.pfhz;
  temp.pfhw = s1.pfhw + s2.pfhw;
  return temp;

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern Ppoint4 tp_sub4(C(Ppoint4) s1, C(Ppoint4) s2)
PreANSI(Ppoint4 s1)
PreANSI(Ppoint4 s2)
/*
subtract two 4D vectors, and return the result
*/
{
  Ppoint4 temp;

  temp.pfhx = s1.pfhx - s2.pfhx;
  temp.pfhy = s1.pfhy - s2.pfhy;
  temp.pfhz = s1.pfhz - s2.pfhz;
  temp.pfhw = s1.pfhw - s2.pfhw;
  return temp;

}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern Ppoint4 tp_linsum4(C(Pfloat) a, C(Ppoint4) s1, 
                          C(Pfloat) b, C(Ppoint4) s2)
PreANSI(Pfloat a) 
PreANSI(Ppoint4 s1) 
PreANSI(Pfloat b) 
PreANSI(Ppoint4 s2)
/*
linear sum of two 4D vectors(a.s1 +b.s2), and return the result
*/
{
  Ppoint4 temp;

  temp.pfhx = (a*(s1.pfhx)) + (b*(s2.pfhx));
  temp.pfhy = (a*(s1.pfhy)) + (b*(s2.pfhy));
  temp.pfhz = (a*(s1.pfhz)) + (b*(s2.pfhz));
  temp.pfhw = (a*(s1.pfhw)) + (b*(s2.pfhw));

  return temp;

}

/*+----------------------------------------------------------------------------*/
/*function:external*/
extern Pfloat tp_dotv4(C(Ppoint4) v, C(Ppoint4) p)
PreANSI(Ppoint4 v)
PreANSI(Ppoint4 p)
/*
Compute the dot product of two position vectors
*/
{
  return (Pfloat) (v.pfhx * p.pfhx + v.pfhy * p.pfhy + v.pfhz * p.pfhz + v.pfhw * p.pfhw);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern Pfloat tp_modv4(C(Ppoint4) P1)
PreANSI(Ppoint4 P1)
/* 
Compute 4D modulus 
*/
{
  return (Pfloat) sqrt(P1.pfhx * P1.pfhx + P1.pfhy * P1.pfhy +
	     P1.pfhz * P1.pfhz + P1.pfhw * P1.pfhw);
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_transpose(C(PR_nurb *)s1)
PreANSI(PR_nurb *s1)
/*
Transpose the NURB
*/
{
  PR_knots *temp;
  PR_pts *tt;
  Pint i, j, mm, mm1;

/*Get space for the new points*/

  memset(&tt, 0, sizeof(PR_pts *));
  nrb_allocatepts(s1->pf_u.pf_n * s1->pf_v.pf_n,&tt);


  mm = s1->pf_u.pf_n;
  for (i = 1; i <= mm; i++) {   /*Copy the points transposing on the way*/
    mm1 = s1->pf_v.pf_n;
    for (j = 1; j <= mm1; j++)
      tt->pts[ni(j, i, s1->pf_v.pf_n) - 1] =
	s1->pf_ppp->pts[ni(i, j, s1->pf_u.pf_n) - 1];
  }

                                          /*xchange knot vectors, swap pointers!*/
  temp = s1->pf_u.pf_kk;
  s1->pf_u.pf_kk = s1->pf_v.pf_kk;
  s1->pf_v.pf_kk = temp;

  swap(&s1->pf_u.pf_n, &s1->pf_v.pf_n);   /*interchange # ctrl pts*/
  swap(&s1->pf_u.pf_nt, &s1->pf_v.pf_nt); /*Interchange # knots*/
  swap(&s1->pf_u.pf_k, &s1->pf_v.pf_k);   /*interchange order*/

  nrb_deallocatepts(&s1->pf_ppp);         /*Get rid of the old points*/
  s1->pf_ppp = tt;                        /*Copy pointer to new ones*/
}
#ifdef VMS
#include "[-.nrb_source]xpose.c"
#else
#include "./xpose.c"
#endif

/*---------------------------------------------------------------------------*/
/*function:external*/
extern Pint nrb_numsubdivs(C(PR_nurb *)nrb, C(PE_dir) xory, C(Pfloat) tol)
PreANSI(PR_nurb *nrb) 
PreANSI(PE_dir xory) 
PreANSI(Pfloat tol)
/*
This routine determines the level of refinement to be applied to the NURB.
*/
{
  Pint i, j;
  Pfloat max_norm, min_w, vec_norm;
  Pint mm, mm1;

  min_w = LONG_MAX;
  max_norm = 0.0;
  if (xory == PVudir) {   /* dealing with x ctrl points */
    mm = nrb->pf_v.pf_n;
    for (j = 1; j <= mm; j++) {
      mm1 = nrb->pf_u.pf_n;
      for (i = 1; i < mm1; i++) {
	vec_norm = tp_modv4(
	    tp_sub4(nrb->pf_ppp->pts[ni(i, j, nrb->pf_u.pf_n) - 1],
		    nrb->pf_ppp->pts[ni(i + 1, j, nrb->pf_u.pf_n) - 1]));
	max_norm = P_max(max_norm, vec_norm);
      }
    }
    mm = nrb->pf_v.pf_n;
    for (j = 1; j <= mm; j++) {
      mm1 = nrb->pf_u.pf_n;
      for (i = 1; i <= mm1; i++) {
	min_w = P_min(min_w,
	    (Pfloat)nrb->pf_ppp->pts[ni(i, j, nrb->pf_u.pf_n) - 1].pfhw);
      }
    }
  } else {  /*Y direction*/
    mm = nrb->pf_u.pf_n;
    for (i = 1; i <= mm; i++) {
      mm1 = nrb->pf_v.pf_n;
      for (j = 1; j < mm1; j++) {
	vec_norm = tp_modv4(
	    tp_sub4(nrb->pf_ppp->pts[ni(i, j, nrb->pf_u.pf_n) - 1],
		    nrb->pf_ppp->pts[ni(i, j + 1, nrb->pf_u.pf_n) - 1]));
	max_norm = P_max(max_norm, vec_norm);
      }
    }
    mm = nrb->pf_u.pf_n;
    for (i = 1; i <= mm; i++) {
      mm1 = nrb->pf_v.pf_n;
      for (j = 1; j <= mm1; j++) {
	min_w = P_min(min_w,
	    (Pfloat)nrb->pf_ppp->pts[ni(i, j, nrb->pf_u.pf_n) - 1].pfhw);
      }
    }
  }


  if (xory == PVudir)
    i = (long)floor(
	  (Pfloat)sqrt(2.0) * (nrb->pf_u.pf_k - 1) * (max_norm / min_w) / tol + 0.5);
  else
    i = (long)floor(
	  (Pfloat)sqrt(2.0) * (nrb->pf_v.pf_k - 1) * (max_norm / min_w) / tol + 0.5);

  if (i > 50)
    i = 50;

  return (Pint)P_max(P_max(i, nrb->pf_u.pf_k * 2 + 1), nrb->pf_v.pf_k * 2 + 1);
      /*Make sure its sensible!*/
  /* with */
}  /* nrb_numsubdivs */
