/*

                              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.
*/                                            

#include "nurbh.h"
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_evaluate(C(PR_nurb *)s1, C(Pfloat) x, C(PR_nurb *)S2)
PreANSI(PR_nurb *s1) 
PreANSI(Pfloat x) 
PreANSI(PR_nurb *S2)
/*
Evaluate NURB s1 at value x. This routine simply call nrb_bvalue to 
evaluate the NURB at the zeroth derivative.
*/
{
  /*nrb_evaluate*/
  nrb_bvalue(s1, x, 0, S2);   /*Evaluate Zero'th Derivative*/
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_du(C(PR_nurb *)s1, C(Pfloat) x, C(PR_pts **)XX)
PreANSI(PR_nurb *s1) 
PreANSI(Pfloat x) 
PreANSI(PR_pts **XX)
/*
Evaluate NURB s1 at value x. This routine simply call nrb_bvalue to evaluate 
the NURB at the first derivative.
*/
{
  PR_nurb *s2;

  /*nrb_du*/
  s2 = NULL;
  nrb_allocatenurb(&s2);
  nrb_bvalue(s1, x, 1, s2);   /*Evaluate 1st Derivative*/
  *XX = s2->pf_ppp;   /*Return Pointer to the row!*/

  s2->pf_ppp = NULL;
  nrb_deallocatenurb(&s2);
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_partials(C(PR_nurb *)s1, C(PR_pts **)du, C(PR_pts **)dv, 
                         C(PR_pts **)duv)
PreANSI(PR_nurb *s1) 
PreANSI(PR_pts **du) 
PreANSI(PR_pts **dv) 
PreANSI(PR_pts **duv)
/*
Calculate the partial derivatives of the NURB in u, v and uv.
*/
{
  Pint i, j;
  Pfloat x;
  PR_pts *dd;
  Pint mm, mm1;
  PR_pts *tpts;
  Pvector3 v1,v2,v3;
  Ppoint3 p1,p2,p3;

  /*calculates first derivates at the pts corresponding to knots K...n*/
  nrb_allocatepts(s1->pf_u.pf_n * s1->pf_v.pf_n, du);
  nrb_allocatepts(s1->pf_u.pf_n * s1->pf_v.pf_n, dv);
  nrb_allocatepts(s1->pf_u.pf_n * s1->pf_v.pf_n, duv);

  dd = NULL;

  mm = s1->pf_u.pf_n;
  for (i = s1->pf_u.pf_k; i <= mm; i++) {
    x = s1->pf_u.pf_kk->knots[i - 1];
    nrb_du(s1, x, &dd);
    mm1 = s1->pf_v.pf_n;
    /*copy Row to correct place*/
    for (j = 1; j <= mm1; j++)
      (*du)->pts[ni(i, j, s1->pf_u.pf_n) - 1] = dd->pts[j - 1];
  }

  nrb_transpose(s1);
  mm = s1->pf_u.pf_n;
  for (i = s1->pf_u.pf_k; i <= mm; i++) {
    x = s1->pf_u.pf_kk->knots[i - 1];
    nrb_du(s1, x, &dd);
    mm1 = s1->pf_v.pf_n;
    /*copy Row to correct place*/
    for (j = 1; j <= mm1; j++)
      (*dv)->pts[ni(j, i, s1->pf_v.pf_n) - 1] = dd->pts[j - 1];
  }

  nrb_transpose(s1);
  nrb_deallocatepts(&dd);
  mm = s1->pf_u.pf_n * s1->pf_u.pf_n;
  /*Now do normals du X dv*/
  for (i = 0; i < mm; i++) {
    tpts = *duv;
    p1 = ptk_pt4topt3(&((*du)->pts[i]));
    p2 = ptk_pt4topt3(&((*dv)->pts[i]));
    v1 = ptk_pt3tovec3(&p1);
    v2 = ptk_pt3tovec3(&p2);
    v3 = ptk_crossv(&v1,&v2);
    p3 = ptk_vec3topt3(&v3);
    tpts->pts[i] = ptk_pt3topt4(&p3);
  }

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_split(C(PR_nurb *)s1, C(Pfloat) x, C(PR_nurb *)s2, 
                      C(PR_nurb *)s3)
PreANSI(PR_nurb *s1) 
PreANSI(Pfloat x) 
PreANSI(PR_nurb *s2) 
PreANSI(PR_nurb *s3)
/*
Split a NURB in two.
*/
{
  Pint i1, ii, msplit, mbreak, mflag, ll, iin, n1, n2, ir;
  PR_dir kk;
  PR_nurb *ss;
  Pint mm;
  PR_dir *tdir;


  /*nrb_split*/
  ss = NULL;
  nrb_allocatenurb(&ss);
  nrb_interv(&s1->pf_u, x, &mbreak, &mflag);

  kk.pf_kk = NULL;
  kk.pf_nt = s1->pf_u.pf_nt + s1->pf_u.pf_k;

  nrb_allocateknots(&kk);
  iin = 0;
  for (ii = 0; ii < mbreak; ii++) {
    iin++;
    kk.pf_kk->knots[iin - 1] = s1->pf_u.pf_kk->knots[ii];
  }

  /*Make sure its multiplicity is pf_u.pf_k*/

  ir = s1->pf_u.pf_k;
  mm = s1->pf_u.pf_k;
  for (ii = 1; ii <= mm; ii++) {
    if (ptk_equal(x, s1->pf_u.pf_kk->knots[mbreak - ii])) {
      ir--;
    }
  }

  for (ii = 1; ii <= ir; ii++) {
    iin++;
    kk.pf_kk->knots[iin - 1] = x;
  }
  msplit = iin;

  mm = s1->pf_u.pf_nt;
  for (ii = mbreak; ii < mm; ii++) {
    iin++;
    kk.pf_kk->knots[iin - 1] = s1->pf_u.pf_kk->knots[ii];
  }

  nrb_clear(ss);
  nrb_osloc(s1, ss);   /*new kx,ntx*/

  /*Initialise S2 - interchnaging x and Y*/

  /*new ky,nty*/
  nrb_init(s1->pf_v.pf_k, iin, s1->pf_u.pf_k, s1->pf_u.pf_nt, s2);
      /*destination*/

  /*Copy s1 Y knots to s2 X knots*/
  nrb_copyknots(&s1->pf_v, &s2->pf_u);

  /*Copy  'half' temporary knots to Y knots*/
  for (ii = 0; ii < msplit; ii++)   /*new kx,ntx*/
    s2->pf_v.pf_kk->knots[ii] = kk.pf_kk->knots[ii];


  /*Initialize S3*/
  /*new ky,nty*/
  nrb_init(s1->pf_v.pf_k, iin, s1->pf_u.pf_k, s1->pf_u.pf_nt, s3);
      /*destination*/

  tdir = &s3->pf_u;

  mm = tdir->pf_nt;
  /*Copy s1 Y knots to s3 X knots*/
  for (ii = 0; ii < mm; ii++)
    tdir->pf_kk->knots[ii] = s1->pf_v.pf_kk->knots[ii];
  /*Copy  'half' temporary knots to Y knots*/
  n2 = msplit - s1->pf_u.pf_k;
  tdir = &s2->pf_v;

  for (ii = n2; ii <= iin; ii++)
    tdir->pf_kk->knots[ii - n2] = kk.pf_kk->knots[ii - 1];
  n1 = mbreak;
  n2 = iin - s1->pf_u.pf_k - mbreak + 1;

  mm = s1->pf_v.pf_n;
  for (ll = 1; ll <= mm; ll++) {
    for (ii = 1; ii <= n1; ii++)
      s2->pf_ppp->pts[ni(ll, ii, s2->pf_u.pf_n) - 1] =
	ss->pf_ppp->pts[ni(ll, ii, ss->pf_u.pf_n) - 1];
  }

  for (ii = 1; ii <= n2; ii++) {
    i1 = ii + mbreak - 1;
    s3->pf_ppp->pts[ni(ll, ii, s3->pf_u.pf_n) - 1] =
      ss->pf_ppp->pts[ni(ll, i1, ss->pf_u.pf_n) - 1];
  }

  nrb_deallocateknots(&kk);
  nrb_deallocatenurb(&ss);

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_join(C(PR_nurb *)nrb1, C(PR_nurb *)nrb2, C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb2) 
PreANSI(PR_nurb *nrb)
/*
Join two NURBs -- This is not currently implemented.
*/
{
  /*join aPint RHedge of 1 and Left hand edge of 2*/
  *nrb = *nrb1;
  *nrb = *nrb2;
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_xform(C(PR_nurb *)nrb1, C(Pmatrix3) xx)
PreANSI(PR_nurb *nrb1) 
PreANSI(Pmatrix3 xx)
/*
Transfrom the control points of nrb1 using the matrix xx.
*/
{
  Pint i, mm;

  mm = nrb1->pf_u.pf_n * nrb1->pf_v.pf_n;
  for (i = 0; i < mm; i++) {
    nrb1->pf_ppp->pts[i] = ptk_transform4(xx, &(nrb1->pf_ppp->pts[i]));
  }
}  /* nrb_xform */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_extrude(C(PR_nurb *)nrb1, C(Ppoint3) w, C(PR_nurb *)nrb, 
                        C(Pint *)error)
PreANSI(PR_nurb *nrb1) 
PreANSI(Ppoint3 w) 
PreANSI(PR_nurb *nrb) 
PreANSI(Pint *error)
/*
Perform a linear extrusion on nrb along the vector w to produce nrb.
*/
{
  /*Do a linear extrusion along the vector w*/
  /*see Curve and Surface Construction using Rational B-splines
    Piegl and Tiller CAD Vol 19 #9 November 1987 pp 485-498
    Note we use 1..nus instead of their 0..nthem i.e. nus = nthem +1

   Assumptions?
   ------------
   */
  Pint ii, i, j;
  Pfloat wcoord;
  PR_pts *tpts;
  Pint mm;

  *error = 0;  /* No error*/

  nrb->pf_u.pf_k  = nrb1->pf_u.pf_k;   /*copy order*/
  nrb->pf_u.pf_n  = nrb1->pf_u.pf_n;   /*# control points*/
  nrb->pf_u.pf_nt = nrb1->pf_u.pf_nt;  /*# knots*/

  nrb->pf_v.pf_k  = 2;   /*linear*/
  nrb->pf_v.pf_n  = 2;   /*two control points*/
  nrb->pf_v.pf_nt = 4;   /*four knots*/

  nrb_allocateknots(&nrb->pf_v);

  nrb->pf_v.pf_kk->knots[0] = 0.0;
  nrb->pf_v.pf_kk->knots[1] = 0.0;
  nrb->pf_v.pf_kk->knots[2] = 1.0;
  nrb->pf_v.pf_kk->knots[3] = 1.0;

  nrb_copyknots(&nrb1->pf_u, &nrb->pf_u);   /*Copy the X knots*/
  nrb_allocatepts(nrb->pf_u.pf_n * nrb->pf_v.pf_n, &nrb->pf_ppp);

  tpts = nrb->pf_ppp;
  mm = nrb->pf_u.pf_n;
  for (i = 0; i < mm; i++) {
    for (j = 1; j <= 2; j++) {
      ii = ni(i + 1, j, nrb->pf_u.pf_n);

      /*change added 11/02/90 Steve Larkin :
	     Have multiplied the term wcoord to the extrusion vector
	     before adding to the curve coordinates. Not really sure
	     this is the correct fix, but it works ! */

      wcoord = nrb1->pf_ppp->pts[i].pfhw;

tpts->pts[ii - 1].pfhx = nrb1->pf_ppp->pts[i].pfhx + (j - 1) * w.pfhx * wcoord;
tpts->pts[ii - 1].pfhy = nrb1->pf_ppp->pts[i].pfhy + (j - 1) * w.pfhy * wcoord;
tpts->pts[ii - 1].pfhz = nrb1->pf_ppp->pts[i].pfhz + (j - 1) * w.pfhz * wcoord;
tpts->pts[ii - 1].pfhw = nrb1->pf_ppp->pts[i].pfhw;
/*pts[ii].pfhw := nrb1.pf_ppp^.pts[i].pfhw + (j-1)*1.0    * wcoord;*/
    }
  }
}
/*------------------------------------------------------------------------------*/

#define addknot(vv) V.pf_kk->knots[(++V.pf_nt) -1]=vv
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_elevateknots(C(PR_dir *)s1, C(PR_dir *)s2)
PreANSI(PR_dir *s1) 
PreANSI(PR_dir *s2)
/*
Routine it raises the degree of each of the knot in k1 by one, but 
ensures that no knot is repeated more than k1.pf_k+1.
*/
{
  PR_dir V;
  Pint i, ii, repeats;
  Pfloat lastknot;
  Pint size, mm;

  /*nrb_elevate_knots*/
  V.pf_kk = NULL;   /*Get some sapce for the knots*/
  V.pf_nt = s1->pf_nt * 2;
  size = V.pf_nt;   /*Save so that can be deallocated properly*/
  nrb_allocateknots(&V);

  lastknot = s1->pf_kk->knots[0];   /*Save the knot*/
  repeats = 1;
  V.pf_nt = 0;   /*No new knots*/

  mm = s1->pf_nt;
  for (i = 1; i < mm; i++) {   /*i loop*/
    if (ptk_equal(lastknot, s1->pf_kk->knots[i])) {
      repeats++;
    } else {
      repeats++;   /*Add the extra knot */
      repeats = P_min(s1->pf_k + 1, repeats);
	  /*Knot multiplicity is limited to kx+1*/
      for (ii = 1; ii <= repeats; ii++) addknot(lastknot);

      repeats = 1;   /*Reset lastknot and repeats*/
      lastknot = s1->pf_kk->knots[i];
    }
  }

  /*Output `last' knot*/
  repeats++;
  repeats = P_min(s1->pf_k + 1, repeats);
  for (ii = 1; ii <= repeats; ii++)
    addknot(lastknot);

  /*copy knot to new vector*/
  s2->pf_nt = V.pf_nt;
  s2->pf_k = s1->pf_k + 1;   /*Increase degree*/
  s2->pf_n = s2->pf_nt - s2->pf_k;

  nrb_allocateknots(s2);
  mm = V.pf_nt;
  for (i = 0; i < mm; i++)
    s2->pf_kk->knots[i] = V.pf_kk->knots[i];

  V.pf_nt = size;
  nrb_deallocateknots(&V);

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_elevate(C(boolean) transpose, C(PR_nurb *)s1, C(PR_nurb *)s2)
PreANSI(boolean transpose) 
PreANSI(PR_nurb *s1) 
PreANSI(PR_nurb *s2)
/*
Raise the degree (by 1) of S1 (a B-spline Surface) and store the result in S2
*/
{
  Pint i, j, kk, k, r, ist, iend, mflag;
  Pfloat y1, y2, yi, yik, uden;
/*  Ppoint4 cc[PC_NrbMax], cw[PC_NrbMax];*/
  PR_pts *tcc,*tcw;
  
  Ppoint4 t1, t2, hzero;
  Pint ll;
  PR_knots *temp;
  Pint mm, mm1, mm2;

  hzero = tp_litv4(0.0, 0.0, 0.0, 0.0);
  
  kk = s1->pf_u.pf_k;
  
  memset(&tcc, 0, sizeof(PR_pts *));
  nrb_allocatepts(kk,&tcc);
  memset(&tcw, 0, sizeof(PR_pts *));
  nrb_allocatepts(kk,&tcw);     /*Get space for the new points*/

  nrb_copy(s1, s2);   /*Initialise the new NURB*/

  nrb_elevateknots(&s1->pf_u, &s2->pf_u);   /*Elevate the Knots*/

  nrb_allocatepts(s2->pf_u.pf_n * s2->pf_v.pf_n, &s2->pf_ppp);
  /*Get more space for the  control points*/
  
  mm = P_max(1, s2->pf_v.pf_n);
  for (ll = 1; ll <= mm; ll++) {   /*Loop over each Y row*/
    mm1 = s2->pf_u.pf_n;
    for (j = 1; j <= mm1; j++) {   /*Along each row*/
      /*1. Find r such that knots(r)= knots'(j),knots(r+1)*/
      nrb_interv(&s1->pf_u, s2->pf_u.pf_kk->knots[j - 1], &r, &mflag);

      /*2. For i:=1 to m do cc(i) := cw(i) := 0*/
      for (i = 0; i < kk; i++) {   /*i loop*/
	tcc->pts[i] = hzero;
	tcw->pts[i] = hzero;
      }
      /*3 for i:= max(1,m+1-r), min(m,n+m-r)
	  do cc(i) := cw(i) := ctrlpts(i+r-m)*/

      mm2 = P_min(kk, kk + s1->pf_u.pf_n - r);
      for (i = P_max(1, kk - r + 1); i <= mm2; i++) {
	tcc->pts[i - 1] = s1->pf_ppp->pts[ni(i + r - kk, ll, s1->pf_u.pf_n) - 1];
	tcw->pts[i - 1] = s1->pf_ppp->pts[ni(i + r - kk, ll, s1->pf_u.pf_n) - 1];
      }
      /*4 for k := m-1 step -1 until 1
	 a. y1 = knots'(j + k )
	 b. y2 = knots'(j + k +1)
	 c. for i:= min(m,n+2m-k-r) step -1
		    until max(m-k+1,1-r+m)
	    yi = knots(r+i-m)
	    yik = knots(r+i+k-m)
	    den = yik - yi
	    cc(i) = ((y2 - yi)) * cc(i) +
		    (yik - y2) * cc(i-1))/den
	    cw(i) = ((y1 - yi)) * cw(i) +
		    (yik - y1) * cw(i-1))/den
		    + cc(i)
      */
      for (k = kk - 1; k >= 1; k--) {   /*k loop*/
	y1 = s2->pf_u.pf_kk->knots[j + k - 1];
	y2 = s2->pf_u.pf_kk->knots[j + k];
	ist = P_min(kk, s1->pf_u.pf_n + kk * 2 - k - r);
	iend = P_max(kk - k + 1, kk - r + 1);

	for (i = ist; i >= iend; i--) {   /*ii loop*/
	  yi = s1->pf_u.pf_kk->knots[r + i - kk - 1];
	  yik = s1->pf_u.pf_kk->knots[r + i + k - kk - 1];
	  uden = 1.0 / (yik - yi);
	  t1 = tp_scale4(tcc->pts[i - 1], (y2 - yi) * uden);
	  t2 = tp_scale4(tcc->pts[i - 2], (yik - y2) * uden);
	  tcc->pts[i - 1] = tp_add4(t1, t2);
	  t1 = tp_scale4(tcw->pts[i - 1], (y1 - yi) * uden);
	  t2 = tp_scale4(tcw->pts[i - 2], (yik - y1) * uden);
	  t1 = tp_add4(t1, t2);
	  tcw->pts[i - 1] = tp_add4(t1, tcc->pts[i - 1]);
	}
      }

      /*5 newctrlpt(j)= cw(m)/m*/
      /* Transpose!*/
      if (transpose)
	s2->pf_ppp->pts[ni(ll, j, s2->pf_v.pf_n) - 1] = tp_scale4(tcw->pts[kk - 1],
								  1.0 / kk);
      else
	s2->pf_ppp->pts[ni(j, ll, s2->pf_u.pf_n) - 1] = tp_scale4(tcw->pts[kk - 1],
								  1.0 / kk);

    }
    /*j and ll loop and with s2*/
  }

  nrb_deallocatepts(&tcc);
  nrb_deallocatepts(&tcw);

  if (!transpose)  /*transpose the knot vectors*/
    return;

  temp = s2->pf_u.pf_kk;   /*swap pointers!*/
  s2->pf_u.pf_kk = s2->pf_v.pf_kk;
  s2->pf_v.pf_kk = temp;
  nrb_interchange(s2);   /*swap fields*/
}  /* nrb_elevate */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_ruled(C(PR_nurb *)nrb1, C(PR_nurb *)nrb2, C(PR_nurb *)nrb, 
                      C(Pint *)error)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb2) 
PreANSI(PR_nurb *nrb) 
PreANSI(Pint *error)
/*
Do a ruled surface between curve1 and curve2.
see Curve and Surface Construction using Rational B-splines
    Piegl and Tiller CAD Vol 19 #9 November 1987 pp 485-498
    Note we use 1..nus instead of their 0..nthem i.e. nus = nthem +1
*/
{

  PR_nurb *nrb11, *nrb22, *yyptr, *xxptr, *nrb111, *nrb222;
  Pint i, mm;
  PR_dir *tdir;
  PR_knots *tknots;

  *error = 0;  /* No error*/

  /*find the higher degree curve, and raise the degree to that*/
  nrb11 = NULL;
  nrb22 = NULL;
  nrb111 = NULL;
  nrb222 = NULL;

  nrb_allocatenurb(&nrb11);
  nrb_allocatenurb(&nrb22);
  nrb_allocatenurb(&nrb111);
  nrb_allocatenurb(&nrb222);

  nrb->pf_u.pf_k = P_max(nrb1->pf_u.pf_k, nrb2->pf_u.pf_k);

  nrb_copy(nrb1, nrb11);
  nrb11->pf_v.pf_n = 1;   /*Force it to be a curve*/
  nrb11->pf_v.pf_nt = 1;
  nrb11->pf_v.pf_k = 0;
  mm = nrb->pf_u.pf_k;
 
  xxptr = nrb1;
  for (i = nrb1->pf_u.pf_k + 1; i <= mm; i++) {
    nrb_elevate(FALSE, xxptr, nrb11);
    yyptr = xxptr;
    xxptr = nrb11;
    nrb11 = yyptr;
  }

  nrb_copy(nrb2, nrb22);
  nrb22->pf_v.pf_n = 1;   /*Force it to be a curve*/
  nrb22->pf_v.pf_nt = 1;
  nrb22->pf_v.pf_k = 0;

  mm = nrb->pf_u.pf_k;
  xxptr = nrb2;

  for (i = nrb2->pf_u.pf_k + 1; i <= mm; i++) {
    nrb_elevate(FALSE, xxptr, nrb22);
    yyptr = xxptr;
    xxptr = nrb22;
    nrb22 = yyptr;
  }

  nrb_unionknots(&nrb11->pf_u, &nrb22->pf_u, &nrb->pf_u);
      /*pf_v.pf_k =0, pf_v.pf_nt = 1*/


  tdir = &nrb->pf_u;

  nrb_init(tdir->pf_k, tdir->pf_nt, 0, 1, nrb111);   /*new kx,ntx,ky,nty*/
  nrb_init(tdir->pf_k, tdir->pf_nt, 0, 1, nrb222);   /*new kx,ntx,ky,nty*/
  mm = nrb->pf_u.pf_nt;
  /*copy common Knot vector*/
  for (i = 0; i < mm; i++) {
    nrb111->pf_u.pf_kk->knots[i] = nrb->pf_u.pf_kk->knots[i];
    nrb222->pf_u.pf_kk->knots[i] = nrb->pf_u.pf_kk->knots[i];
  }

  nrb_osloc(nrb11, nrb111);
  nrb_osloc(nrb22, nrb222);

  nrb_interchange(nrb111);   /*Put the #knots and ctrl pts back*/
  nrb_interchange(nrb222);

  /*Now make up the Surface*/
  nrb->pf_v.pf_k = 2;   /*linear*/
  nrb->pf_v.pf_n = 2;   /*two control points*/
  nrb->pf_v.pf_nt = 4;   /*four knots*/
  nrb->pf_u.pf_n = nrb111->pf_u.pf_n;

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

  nrb_allocateknots(&nrb->pf_v);

  tknots = nrb->pf_v.pf_kk;

  tknots->knots[0] = 0.0;
  tknots->knots[1] = 0.0;
  tknots->knots[2] = 1.0;
  tknots->knots[3] = 1.0;
  mm = nrb->pf_u.pf_n;
  for (i = 1; i <= mm; i++) {
    nrb->pf_ppp->pts[i - 1] = nrb111->pf_ppp->pts[i - 1];
    nrb->pf_ppp->pts[ni(i, 2, nrb->pf_u.pf_n) - 1] = nrb222->pf_ppp->pts[i - 1];
  }

  nrb_deallocatenurb(&nrb11);
  nrb_deallocatenurb(&nrb22);
  nrb_deallocatenurb(&nrb111);
  nrb_deallocatenurb(&nrb222);
}

/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_revolve(C(PR_nurb *)nrb1, C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb)
/*
Do a surface of revolution about the z axis.
    Assumes curve is in x-z plane
see Curve and Surface Construction using Rational B-splines
    Piegl and Tiller CAD Vol 19 #9 November 1987 pp 485-498
    Note we use 1..nus instead of their 0..nthem i.e. nus = nthem +1

   Assumptions?
   ------------
   This will use the (3D) projection of CURVE onto y= 0 plane
*/
{
  Pint i, j;
  Pfloat t, w, radius;
  Pint ii, jj;
  Pfloat ww;
  Pint mm;

  ww = (Pfloat) sqrt(2.0) * 0.5;   /*magic constant*/
  /*Copy X stuff to Y direction*/
  nrb->pf_v.pf_k = nrb1->pf_u.pf_k;   /*copy order*/
  nrb->pf_v.pf_n = nrb1->pf_u.pf_n;   /*# control poitns*/
  nrb->pf_v.pf_nt = nrb1->pf_u.pf_nt;   /*#knots*/

  nrb->pf_u.pf_k = 3;   /*quadratic*/
  nrb->pf_u.pf_n = 9;   /*9 control points*/
  nrb->pf_u.pf_nt = 12;   /*12 knots*/

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

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

  nrb->pf_u.pf_kk->knots[0] = 0.0;
  jj = 2;
  t = 0.0;
  for (ii = 1; ii <= 5; ii++) {
    nrb->pf_u.pf_kk->knots[jj - 1] = t;
    nrb->pf_u.pf_kk->knots[jj] = t;
    t += 0.25;
    jj += 2;
  }
  nrb->pf_u.pf_kk->knots[11] = 1.0;

  i = 1;   /*copy profile*/
  mm = nrb->pf_v.pf_n;
  for (j = 1; j <= mm; j++) {
    nrb->pf_ppp->pts[ni(1, j, nrb->pf_u.pf_n) - 1] = nrb1->pf_ppp->pts[j - 1];
    nrb->pf_ppp->pts[ni(1, j, nrb->pf_u.pf_n) - 1].pfhy = 0.0;
	/*project profile onto y=0 plane*/
  }

  mm = nrb->pf_v.pf_n;
  /*Sweep each into a square*/
  for (j = 1; j <= mm; j++) {   /*tdir*/
    for (i = 2; i <= 9; i++) {   /*Copy Z's*/
      nrb->pf_ppp->pts[ni(i, j, nrb->pf_u.pf_n) - 1].pfhz =
	nrb->pf_ppp->pts[ni(1, j, nrb->pf_u.pf_n) - 1].pfhz;
    }

    radius = nrb->pf_ppp->pts[ni(1, j, nrb->pf_u.pf_n) - 1].pfhx;
    w = nrb->pf_ppp->pts[ni(1, j, nrb->pf_u.pf_n) - 1].pfhw;

    jj = 2;
    for (ii = 1; ii <= 4; ii++) {
      nrb->pf_ppp->pts[ni(jj, j, nrb->pf_u.pf_n) - 1].pfhz *= ww;
      nrb->pf_ppp->pts[ni(jj, j, nrb->pf_u.pf_n) - 1].pfhw = ww * w;
      nrb->pf_ppp->pts[ni(jj + 1, j, nrb->pf_u.pf_n) - 1].pfhw = w;
      jj += 2;
    }

    nrb->pf_ppp->pts[ni(2, j, nrb->pf_u.pf_n) - 1].pfhx = radius * ww;
    nrb->pf_ppp->pts[ni(2, j, nrb->pf_u.pf_n) - 1].pfhy = radius * ww;
    nrb->pf_ppp->pts[ni(8, j, nrb->pf_u.pf_n) - 1].pfhx = radius * ww;
    nrb->pf_ppp->pts[ni(8, j, nrb->pf_u.pf_n) - 1].pfhy = -radius * ww;

    nrb->pf_ppp->pts[ni(3, j, nrb->pf_u.pf_n) - 1].pfhx = 0.0;
    nrb->pf_ppp->pts[ni(3, j, nrb->pf_u.pf_n) - 1].pfhy = radius;
    nrb->pf_ppp->pts[ni(7, j, nrb->pf_u.pf_n) - 1].pfhx = 0.0;
    nrb->pf_ppp->pts[ni(7, j, nrb->pf_u.pf_n) - 1].pfhy = -radius;

    nrb->pf_ppp->pts[ni(4, j, nrb->pf_u.pf_n) - 1].pfhx = -radius * ww;
    nrb->pf_ppp->pts[ni(4, j, nrb->pf_u.pf_n) - 1].pfhy = radius * ww;

    nrb->pf_ppp->pts[ni(6, j, nrb->pf_u.pf_n) - 1].pfhx = -radius * ww;
    nrb->pf_ppp->pts[ni(6, j, nrb->pf_u.pf_n) - 1].pfhy = -radius * ww;

    nrb->pf_ppp->pts[ni(5, j, nrb->pf_u.pf_n) - 1].pfhx = -radius;
    nrb->pf_ppp->pts[ni(5, j, nrb->pf_u.pf_n) - 1].pfhy = 0.0;

    nrb->pf_ppp->pts[ni(9, j, nrb->pf_u.pf_n) - 1].pfhx = radius;
    nrb->pf_ppp->pts[ni(9, j, nrb->pf_u.pf_n) - 1].pfhy = 0.0;
  }

}  /* nrb_revolve */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_circle(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Make a unit circle in the XY plane
*/
{
  Pint i;
  Pfloat t;
  Pint ii, jj;
  Pfloat ww;

  ww = (Pfloat) sqrt(2.0) * 0.5;
  nrb->pf_u.pf_k = 3;   /*quadratic*/
  nrb->pf_u.pf_n = 9;   /*9 control points*/
  nrb->pf_u.pf_nt = 12;   /*12 knots*/

  nrb->pf_v.pf_k = 0;
  nrb->pf_v.pf_n = 1;
  nrb->pf_v.pf_nt = 1;

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

  nrb->pf_u.pf_kk->knots[0] = 0.0;
  jj = 2;
  t = 0.0;
  
  for (ii = 1; ii <= 5; ii++) {
    nrb->pf_u.pf_kk->knots[jj - 1] = t;
    nrb->pf_u.pf_kk->knots[jj] = t;
    t += 0.25;
    jj += 2;
  }

  nrb->pf_u.pf_kk->knots[11] = 1.0;

  nrb->pf_v.pf_kk->knots[0] = 0.0;

  jj = 2;
  nrb->pf_ppp->pts[0].pfhw = 1;
  for (ii = 1; ii <= 4; ii++) {
    nrb->pf_ppp->pts[jj - 1].pfhw = ww;
    nrb->pf_ppp->pts[jj].pfhw = 1.0;
    jj += 2;
  }

  for (i = 0; i <= 8; i++) {
    nrb->pf_ppp->pts[i].pfhz = 0.0;
  }

  nrb->pf_ppp->pts[0].pfhx = 1.0;
  nrb->pf_ppp->pts[0].pfhy = 0.0;
  nrb->pf_ppp->pts[1].pfhx = 1.0 * ww;
  nrb->pf_ppp->pts[1].pfhy = 1.0 * ww;
  nrb->pf_ppp->pts[7].pfhx = 1.0 * ww;
  nrb->pf_ppp->pts[7].pfhy = -1.0 * ww;
  nrb->pf_ppp->pts[2].pfhx = 0.0;
  nrb->pf_ppp->pts[2].pfhy = 1.0;
  nrb->pf_ppp->pts[6].pfhx = 0.0;
  nrb->pf_ppp->pts[6].pfhy = -1.0;
  nrb->pf_ppp->pts[3].pfhx = -1.0 * ww;
  nrb->pf_ppp->pts[3].pfhy = 1.0 * ww;
  nrb->pf_ppp->pts[5].pfhx = -1.0 * ww;
  nrb->pf_ppp->pts[5].pfhy = -1.0 * ww;
  nrb->pf_ppp->pts[4].pfhx = -1.0;
  nrb->pf_ppp->pts[4].pfhy = 0.0;
  nrb->pf_ppp->pts[8].pfhx = 1.0;
  nrb->pf_ppp->pts[8].pfhy = 0.0;   

}  /* nrb_circle */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_arc(C(Pfloat) theta1, C(Pfloat) theta2, C(PR_nurb *)nrb)
PreANSI(Pfloat theta1)
PreANSI(Pfloat theta2) 
PreANSI(PR_nurb *nrb)
/*
Create a circle centred at (0,0,0), radius 1 in XY Plane
*/
{
  PR_nurb *nrb1, *nrb2, *nrb3;
  Pfloat z, t1, t2;

  nrb1 = NULL;
  nrb2 = NULL;
  nrb3 = NULL;
  nrb_allocatenurb(&nrb1);
  nrb_allocatenurb(&nrb2);
  nrb_allocatenurb(&nrb3);

  nrb_circle(nrb1);
  /*find t1 and t2 from theta 1 and 2*/
  /* t1= 0 at (1,0) - horizontal radius, and 0.25 at (0,1) - vertical*/

  z = sin(theta1 * PC_DTOR) / P_max(cos(theta1 * PC_DTOR), (double)ptkcpceps);

  if (ptk_equal(z, 0.0)) {
    t1 = 0.0;
  } else
    t1 = P_max((sqrt(1 * z * z) - 1) / z, (-1 - sqrt(1 * z * z)) / z);

  z = sin(theta2 * PC_DTOR) / P_max(cos(theta2 * PC_DTOR), (double)ptkcpceps);
  if (ptk_equal(z, 0.0)) {
    t2 = 0.0;
  } else
    t2 = P_max((sqrt(1 * z * z) - 1) / z, (-1 - sqrt(1 * z * z)) / z);

  nrb_split(nrb1, t1, nrb2, nrb3);
  nrb_transpose(nrb3);

  nrb_split(nrb3, t2, nrb, nrb2);
  nrb_transpose(nrb);

  nrb_deallocatenurb(&nrb1);
  nrb_deallocatenurb(&nrb2);
  nrb_deallocatenurb(&nrb3);

}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_sphere(C(PR_nurb *)nrb)
PreANSI(PR_nurb *nrb)
/*
Produce a NURB representation of a sphere.
*/
{
  PR_nurb *nrb1;
  Pfloat ww;

  /*creatre an arc of 180 degrees*/
  nrb1 = NULL;
  nrb_allocatenurb(&nrb1);

  ww = (Pfloat) sqrt(2.0) * 0.5;
  nrb1->pf_u.pf_k = 3;   /*quadratic*/
  nrb1->pf_u.pf_n = 5;
  nrb1->pf_u.pf_nt = 8;

  nrb1->pf_v.pf_n = 1;
  nrb1->pf_v.pf_nt = 1;
  nrb1->pf_v.pf_k = 0;

  /*nrb_revolve - allocates these points!*/
  nrb_allocateknots(&nrb1->pf_u);
  nrb_allocateknots(&nrb1->pf_v);
  nrb_allocatepts(nrb1->pf_u.pf_n * nrb1->pf_v.pf_n, &nrb1->pf_ppp);

  nrb1->pf_u.pf_kk->knots[0] = 0.0;
  nrb1->pf_u.pf_kk->knots[1] = 0.0;
  nrb1->pf_u.pf_kk->knots[2] = 0.0;
  nrb1->pf_u.pf_kk->knots[3] = 0.5;
  nrb1->pf_u.pf_kk->knots[4] = 0.5;
  nrb1->pf_u.pf_kk->knots[5] = 1.0;
  nrb1->pf_u.pf_kk->knots[6] = 1.0;
  nrb1->pf_u.pf_kk->knots[7] = 1.0;

  nrb1->pf_v.pf_kk->knots[0] = 0.0;
  nrb1->pf_ppp->pts[0] = tp_litv4(0.0, 0.0, -1.0,1.0);
  nrb1->pf_ppp->pts[1] = tp_litv4(ww, 0.0, -ww,1.0);
  nrb1->pf_ppp->pts[2] = tp_litv4(1.0, 0.0, 0.0,1.0);
  nrb1->pf_ppp->pts[3] = tp_litv4(ww, 0.0, ww,1.0);
  nrb1->pf_ppp->pts[4] = tp_litv4(0.0, 0.0, 1.0,1.0);

  nrb1->pf_ppp->pts[0].pfhw = 1.0;
  nrb1->pf_ppp->pts[1].pfhw = ww;
  nrb1->pf_ppp->pts[2].pfhw = 1.0;
  nrb1->pf_ppp->pts[3].pfhw = ww;
  nrb1->pf_ppp->pts[4].pfhw = 1.0;   


  nrb_revolve(nrb1, nrb);

  nrb_deallocatenurb(&nrb1);
}  /* nrb_sphere */
