/*

                              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_makeknots(C(Pfloat) tmin, C(Pfloat) tmax, C(PR_dir *)kk)
PreANSI(Pfloat tmin) 
PreANSI(Pfloat tmax) 
PreANSI(PR_dir *kk)
/*
Produce n knots between tmin and tmax, with tmin repeated k times and tmax
    repeated k  time
*/
{
  Pint nn, i;
  Pfloat dt;
  Pint mm;

  /* Replicate start and end knots*/

  nn = kk->pf_nt - kk->pf_k;
  nrb_allocateknots(kk);

  mm = kk->pf_k;
  for (i = 0; i < mm; i++) {
    kk->pf_kk->knots[i] = tmin;
    kk->pf_kk->knots[i + nn] = tmax;
  }
  dt = (tmax - tmin) / (nn - kk->pf_k + 1);
  for (i = kk->pf_k + 1; i <= nn; i++)   
    kk->pf_kk->knots[i - 1] = tmin + (i - kk->pf_k) * dt;
}

/*---------------------------------------------------------------------------*/
extern void nrb_scaleknots(C(float) new_tmin, C(float) new_tmax, C(PR_dir *)kk)
PreANSI(Pfloat new_tmin) 
PreANSI(Pfloat new_tmax) 
PreANSI(PR_dir *kk)
/*
Scale a knots vector between two new values.
*/
{
PR_dir *new_knots;
float old_1, old_2;
int i;

   /* Set the min and max of the knots
   */
   old_1 = kk->pf_kk->knots[0];
   old_2 = kk->pf_kk->knots[kk->pf_nt-1];

   /* Allocate some space for some new knots
   */
   new_knots->pf_k  = kk->pf_k;
   new_knots->pf_n  = kk->pf_n;
   new_knots->pf_nt = kk->pf_nt;

   nrb_allocateknots(new_knots);

   /* Set the knots
   */
   for(i=0; i< kk->pf_nt; i++)
   {
      new_knots->pf_kk->knots[i] 
      = new_tmin 
      + ((kk->pf_kk->knots[i] - old_1)*(new_tmax - new_tmin))/(old_2 -old_1); 
   }

   nrb_copyknots(new_knots,kk);   
   nrb_deallocateknots(new_knots);

   return;
}
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_unionknots(C(PR_dir *)t1, C(PR_dir *)t2, C(PR_dir *)t3)
PreANSI(PR_dir *t1) 
PreANSI(PR_dir *t2) 
PreANSI(PR_dir *t3)
/*
Create a union of knots sequences t1 and t2 and store in t3
*/
{

  boolean done;
  Pint ntx, i, i1, i2;
  Pfloat t;
  PR_dir kk;
  Pint mm;

  i1 = 1;
  i2 = 1;

  /*Allocate a big one!*/
  kk.pf_kk = NULL;
  kk.pf_nt = t1->pf_nt + t2->pf_nt;
  nrb_allocateknots(&kk);
  ntx = 0;

  done = FALSE;
  while (!done) {
    if (ptk_equal(t1->pf_kk->knots[i1 - 1], t2->pf_kk->knots[i2 - 1])) {
      t = t1->pf_kk->knots[i1 - 1];
      i1++;
      i2++;
    } else {
      if (t1->pf_kk->knots[i1 - 1] < t2->pf_kk->knots[i2 - 1]) {
	t = t1->pf_kk->knots[i1 - 1];
	i1++;
      } else {
	t = t2->pf_kk->knots[i2 - 1];
	i2++;
      }
    }

    ntx++;   /*store the next value*/
    kk.pf_kk->knots[ntx - 1] = t;
    /*finished one vector yet*/
    done = (i1 > t1->pf_nt || i2 > t2->pf_nt);
  }
  /* if both are equal copy one and increment both, otherwise
     take the smaller of the two knots*/

  if (i1 <= t1->pf_nt) {
    mm = t1->pf_nt;
    for (i = i1 - 1; i < mm; i++) {
      ntx++;
      kk.pf_kk->knots[ntx - 1] = t1->pf_kk->knots[i];
    }
  } else {
    mm = t2->pf_nt;
    for (i = i2 - 1; i < mm; i++) {
      ntx++;
      kk.pf_kk->knots[ntx - 1] = t2->pf_kk->knots[i];
    }
  }

  /*t3.pf_K - unaltered!*/

  t3->pf_nt = ntx;
  t3->pf_n = t3->pf_nt - t3->pf_k;

  nrb_allocateknots(t3);

  for (i = 0; i < ntx; i++)
    t3->pf_kk->knots[i] = kk.pf_kk->knots[i];

  nrb_deallocateknots(&kk);

}  /* nrb_unionknots */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_differenceknots(C(PR_dir *)t1, C(PR_dir *)t2, C(PR_dir *)t3)
PreANSI(PR_dir *t1) 
PreANSI(PR_dir *t2) 
PreANSI(PR_dir *t3)
/*
Create a knot sequence t3 which contains those knots in t2 which are NOT in t1,
    i.e. assumes t2.nt > t1.nt
*/
{

  Pint ntx, i1, i2;
  PR_dir kk;
  Pint mm;


  /*Allocate a big one!*/
  kk.pf_kk = NULL;
  kk.pf_nt = P_max(t1->pf_nt, t2->pf_nt);
  nrb_allocateknots(&kk);

  i1 = 1;
  ntx = 0;

  mm = t2->pf_nt;
  for (i2 = 0; i2 < mm; i2++) {
    if (ptk_equal(t1->pf_kk->knots[i1 - 1], t2->pf_kk->knots[i2])) {
      if (i1 < t1->pf_nt)
	i1++;
    } else {
      ntx++;   /*store the next value*/
      kk.pf_kk->knots[ntx - 1] = t2->pf_kk->knots[i2];
    }
  }

  /*t3.pf_K - unaltered!*/

  t3->pf_nt = ntx;
  t3->pf_n = t3->pf_nt - t3->pf_k;

  nrb_allocateknots(t3);

  for (i1 = 0; i1 < ntx; i1++)
    t3->pf_kk->knots[i1] = kk.pf_kk->knots[i1];

  nrb_deallocateknots(&kk);

}  /* nrb_differenceknots */
/*---------------------------------------------------------------------------*/

static Pint ilow = 1;

/*function:external*/
extern void interv(C(Pfloat *)xt, C(Pint) lxt, C(Pfloat) x, 
                   C(Pint *)left, C(Pint *)mflag)
PreANSI(Pfloat *xt) 
PreANSI(Pint lxt) 
PreANSI(Pfloat x) 
PreANSI(Pint *left) 
PreANSI(Pint *mflag)
/*
Documentation - to be completed
*/
{
    /* Initialized data */

static Pint istep, middle, ihi;

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* computes left = max( i , 1 .le. i .le. lxt .and. xt(i) .le. x )  . */

/* ******  I N P U T  ****** */
/* xt(real)......a real sequence, of length  lxt , assumed to be nondecreasing*/
/*  lxt(integer)..number of terms in the sequence  xt . */
/* x(real).......the point whose location with respect to the sequence  xt is*/
/*          to be determined. */

/* ******  O U T P U T  ****** */
/*  left, mflag.....both integers, whose value is */

/*   1	  -1	  if		   x .lt.  xt(1) */
/*   i	   0	  if   xt(i)  .le. x .lt. xt(i+1) */
/*  lxt	   1	  if  xt(lxt) .le. x */

/* 	in particular,  mflag = 0 is the 'usual' case.  mflag .ne. 0 */
/* 	indicates that  x  lies outside the halfopen interval */
/* 	xt(1) .le. y .lt. xt(lxt)  . the asymmetric treatment of the */
/* 	interval is due to the decision to make all pp functions cont- */
/* 	inuous from the right. */

/* ******  M E T H O D  ****** */
/*  The program is designed to be efficient in the common situation that */
/*  it is called repeatedly, with  x  taken from an increasing or decrea- */
/*  sing sequence. this will happen, e.g., when a pp function is to be */
/*  graphed. the first guess for  left  is therefore taken to be the val- */
/*  ue returned at the previous call and stored in the  l o c a l  varia- */
/*  ble  ilow . a first check ascertains that  ilow .lt. lxt (this is nec- */
/*  essary since the present call may have nothing to do with the previ- */
/*  ous call). then, if  xt(ilow) .le. x .lt. xt(ilow+1), we set  left = */
/*  ilow  and are done after just three comparisons. */
/*     otherwise, we repeatedly double the difference  istep = ihi - ilow */
/*  while also moving  ilow  and  ihi  in the direction of  x , until */
/* 			xt(ilow) .le. x .lt. xt(ihi) , */
/*  after which we use bisection to get, in addition, ilow+1 = ihi . */
/*  left = ilow  is then returned. */

/* -----------------------------------------------------------------------*/
    /* Parameter adjustments */
    
    --xt;

    /* Function Body */

    ihi = ilow + 1;
    if (ihi < lxt) {
	goto L20;
    }
    if (x >= xt[lxt]) {
	goto L110;
    }
    if (lxt <= 1) {
	goto L90;
    }
    ilow = lxt - 1;
    ihi = lxt;

L20:
    if (x >= xt[ihi]) {
	goto L40;
    }
    if (x >= xt[ilow]) {
	goto L100;
    }

/* 		**** now x .lt. xt(ilow) . decrease  ilow  to capture  x . */

    istep = 1;
L31:
    ihi = ilow;
    ilow = ihi - istep;
    if (ilow <= 1) {
	goto L35;
    }
    if (x >= xt[ilow]) {
	goto L50;
    }
    istep <<= 1;
    goto L31;
/* --------------- */
L35:
    ilow = 1;
    if (x < xt[1]) {
	goto L90;
    }
    goto L50;
/* --------------- */

/*  **** now x .ge. xt(ihi) . increase  ihi  to capture  x . */

L40:
    istep = 1;
L41:
    ilow = ihi;
    ihi = ilow + istep;
    if (ihi >= lxt) {
	goto L45;
    }
    if (x < xt[ihi]) {
	goto L50;
    }
    istep <<= 1;
        goto L41;

L45:
    if (x >= xt[lxt]) {
	goto L110;
    }
    ihi = lxt;

/*  **** now xt(ilow) .le. x .lt. xt(ihi) . narrow the interval. */

L50:
    middle = (ilow + ihi) / 2;
    if (middle == ilow) {
	goto L100;
    }

/* 	note. it is assumed that middle = ilow in case ihi = ilow+1 . */

    if (x < xt[middle]) {
	ihi = middle;
    } else {
	ilow = middle;
    }
    goto L50;
/* --------------- */

/* **** set output and return. */

/* value to left of sequence */

L90:
    *mflag = -1;
    *left = 1;
    return;
/* ------------- */
/* **** Normal return */

L100:
    *mflag = 0;
    *left = ilow;
    return;
/* ------------- */

/* Value to right of sequence */

L110:
    *mflag = 1;
    *left = lxt;
    return;
/* ------------- */
} /* interv_ */

#if DIFFERENT
/*---------------------------------------------------------------------------*/
static Pint ilo = 0;

/**********/
void chop(C(Pfloat) x, C(float *)xt, C(Pint *)left, C(Pint *)mflag, C(Pint) ihi)
PreANSI(Pfloat x)
PreANSI(Pfloat *xt)
PreANSI(Pint   *left)
PreANSI(Pint   *mflag)
PreANSI(Pint   ihi)
{
Pint	middle,
	chopped = 0;

	while (!chopped)
	{
	    middle = (ilo + ihi)/2;
	    if (middle == ilo)
	    {
/*100*/		*mflag = 0;
		*left = ilo + 1;		/* return fortran format */
		chopped = 1;
	    }
	    else
	    {
		if (x < xt[middle])
		{
		    ihi = middle;
		}
		else
		{
		    ilo = middle;
		}
	    }
	}	/* end of while loop */
}	/* end of subroutine chop */
/**********************************/

/****************************************************************/
/*function:external*/
extern void neil_interv(C(float *)xt, C(Pint) lxt, C(float) x, C(Pint *)left, 
                        C(Pint *)mflag)
PreANSI(float *xt) 
PreANSI(Pint lxt) 
PreANSI(float x) 
PreANSI(Pint *left)
PreANSI(Pint *mflag)
/*
Documentation - to be completed
*/
{
Pint	istep,
	finished	= 0,
	ihi;

	--lxt;		/* change array index to C format */

	/* first try the simple case */
	ihi = ilo + 1;
	if (ihi >= lxt)
	{
	    if (x >= xt[lxt])
	    {
/*110*/		*mflag = 1;
		*left = lxt + 1;		/* return fortran format */
		return;
	    }
	    else
	    {
		if (lxt <= 0)
		{
/*90*/		    *mflag = -1;
		    *left = 1;			/* return fortran format */
		    return;
		}
		else
		{
		    ilo = lxt - 1;
		    ihi = lxt;
		}
	    }
	}

/*20*/	if (x >= xt[ihi])
	{
/*40*/	    istep = 1;
	    while (!finished)
	    {
/*41*/		ilo = ihi;
		ihi += istep;
		if (ihi >= lxt)
		{
/*45*/		    if (x >= xt[lxt])
		    {
/*110*/			*mflag = 1;
			*left = lxt + 1;	/* return fortran format */
			finished = 1;
		    }
		    else
		    {
			ihi = lxt;
			chop(x, xt, left, mflag, ihi);
			finished = 1;
		    }
		}
		else
		{
		    if (x < xt[ihi])
		    {
			chop(x, xt, left, mflag, ihi);
			finished = 1;
		    }
		    else
		    {
			istep *= 2;
		    }
		}
	    }	/* end of while loop */
	}
	else
	{
	    if (x >= xt[ilo])
	    {
/*100*/		*mflag = 0;
		*left = ilo + 1;		/* return fortran format */
	    }
	    else
	    {
		istep = 1;
		while (!finished)
		{
/*31*/		    ihi = ilo;
		    ilo -= istep;
		    if (ilo <= 0)
		    {
/*35*/			ilo = 0;
			if (x < xt[0])
			{
/*90*/			    *mflag = -1;
			    *left = 1;		/* return fortran format */
			    finished = 1;
			}
			else
			{
			    chop(x, xt, left, mflag, ihi);
			    finished = 1;
			}
		    }
		    else
		    {
			if (x >= xt[ilo])
			{
			    chop(x, xt, left, mflag, ihi);
			    finished = 1;
			}
			else
			{
			    istep *= 2;
			}
		    }
		}	/* end of while loop */
	    }
	}
}	/* end of subroutine interv */
/********************************************/
#endif
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_interv(C(PR_dir *)knots, C(Pfloat) X, C(Pint *)left, 
                       C(Pint *)mflag)
PreANSI(PR_dir *knots) 
PreANSI(Pfloat X) 
PreANSI(Pint *left) 
PreANSI(Pint *mflag)
/*
Documentation - to be completed
*/
{ 
  Pint nleft, nflag;

  interv(knots->pf_kk->knots, knots->pf_nt, X, left, mflag);

#if DIFFERENT
  neil_interv(knots->pf_kk->knots, knots->pf_nt, X, &nleft, &nflag);
  if ((nleft != *left) || (*mflag != nflag)){
    fprintf(stderr,"Interv and NEIL Interv DISAGREE\n");
    fprintf(stderr,"Interv, left:%d, flag:%d\n",*left,*mflag);
    fprintf(stderr,"Neil Interv, left:%d, flag:%d\n",nleft,nflag);
  }
#endif

  if (*mflag == 1) {   
    nrb_error("NRB_interv - overran");
    *left = knots->pf_n;
    return;
  }
  if (*mflag == -1) {
    nrb_error("NRB_interv - underran");
    *left = knots->pf_k;
  }
}

/*----------------------------------------------------------------------*/
/*function:external*/
extern void nrb_subdivc(C(PR_nurb *)nrb1, C(PR_nurb *)nrb2, C(Pint) mu, 
                        C(Pint) j)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb2) 
PreANSI(Pint mu) 
PreANSI(Pint j)
/*
Calculates the control points corresponding to new knots.
*/
{
  Pint ii, ir, kk,jj;
  Pfloat tmul, t1, t2;
  PR_pts *tt;
  Pint ttsize;
/*  Ppoint4 temp[PC_NrbKmax][PC_NrbMax];*/
  Pint mm, mm1;
  Pint i1,i2,i3;
  
  mm  = nrb1->pf_v.pf_n;
  mm1 = nrb1->pf_u.pf_k;

  memset(&tt, 0, sizeof(PR_pts *));
  ttsize = mm1+1;
  nrb_allocatepts(ttsize * (mu + 1),&tt);      /*Get space for the new points*/

  for (kk = 1; kk <= mm; kk++) {   /*kk loop*/

    for (ii = mu - nrb1->pf_u.pf_k + 1; ii <= mu; ii++)
      tt->pts[ni(0,ii - 1,ttsize)]
            = nrb1->pf_ppp->pts[ni(ii, kk, nrb1->pf_u.pf_n) - 1];

    for (ir = 1; ir < mm1; ir++) {
      for (ii = mu - nrb1->pf_u.pf_k + ir + 1; ii <= mu; ii++) {
	t1 = nrb2->pf_u.pf_kk->knots[j + nrb1->pf_u.pf_k - ir - 1] -
	     nrb1->pf_u.pf_kk->knots[ii - 1];
	t2 = nrb1->pf_u.pf_kk->knots[ii + nrb1->pf_u.pf_k - ir - 1] -
	     nrb2->pf_u.pf_kk->knots[j + nrb1->pf_u.pf_k - ir - 1];
	tmul = 1.0 / (t1 + t2);
        i1 = ni(ir    ,ii - 1,ttsize);
        i2 = ni(ir - 1,ii - 1,ttsize);
        i3 = ni(ir - 1,ii - 2,ttsize);
        tt->pts[i1].pfhx = (t1 * tt->pts[i2].pfhx + t2 * tt->pts[i3].pfhx) * tmul;
	tt->pts[i1].pfhy = (t1 * tt->pts[i2].pfhy + t2 * tt->pts[i3].pfhy) * tmul;
	tt->pts[i1].pfhz = (t1 * tt->pts[i2].pfhz + t2 * tt->pts[i3].pfhz) * tmul;
	tt->pts[i1].pfhw = (t1 * tt->pts[i2].pfhw + t2 * tt->pts[i3].pfhw) * tmul;
      }  /* ii */
    }

    jj = ni(kk, j, nrb2->pf_v.pf_n);
    nrb2->pf_ppp->pts[jj - 1] = tt->pts[ni(nrb1->pf_u.pf_k - 1,mu - 1,ttsize)];

  }
  nrb_deallocatepts(&tt);
  
}  /* nrb_subdivc */
/*-----------------------------------------------------------------------*/
/*function:external*/
extern void nrb_osloc(C(PR_nurb *)nrb1, C(PR_nurb *)nrb2)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_nurb *nrb2)
/*
Performs refinement on nrb1 and places the new NURB in nrb2.
*/
{
  Pint jj, mu, mflag;
  PR_knots *temp;
  Pint mm;

  temp = NULL;

  /*Make sure there is sufficient space*/

  jj = nrb2->pf_u.pf_nt * nrb2->pf_v.pf_nt;

  nrb_allocatepts(jj, &nrb2->pf_ppp);

  mm = nrb2->pf_u.pf_n;
  for (jj = 1; jj <= mm; jj++) {
    nrb_interv(&nrb1->pf_u, nrb2->pf_u.pf_kk->knots[jj - 1], &mu, &mflag);
    nrb_subdivc(nrb1, nrb2, mu, jj);
  }

  /*transpose the knot vectors - swap pointers!*/
  temp = nrb2->pf_u.pf_kk;
  nrb2->pf_u.pf_kk = nrb2->pf_v.pf_kk;
  nrb2->pf_v.pf_kk = temp;

  nrb_interchange(nrb2);   /*swap fields*/
  
}  /* nrb_osloc */
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_boehmc(C(PR_nurb *)nrb1, C(PR_dir *)tnew, C(PR_nurb *)nrb2)
PreANSI(PR_nurb *nrb1) 
PreANSI(PR_dir *tnew) 
PreANSI(PR_nurb *nrb2)
/*
Inserts the knots in tnew into each row (U direction) of
nrb1 returning the answer transposed
in nrb2. It uses Boehm's algorithm, and nrb1 is unmodified.
*/
{
  Pint ii, jj, j, l, oldcpt, newcpt, oldkpt, newkpt, nnew, nold;
  Ppoint4 t1, t2;
  Pfloat beta;
  Pint mm, mm1;
  PR_dir *tdir;
  PR_pts *tpts1;
  Pint mm2;


  /*Copy much of nrb1 to nrb2 transposing on the way*/
  /*Copy knots*/

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

  /*nrb_copyknots(&nrb1->pf_u, &nrb2->pf_v); But must leave space for others!*/

  nrb2->pf_v.pf_nt = nrb1->pf_u.pf_nt + tnew->pf_nt;
  nrb2->pf_v.pf_k  = nrb1->pf_u.pf_k;
  nrb2->pf_v.pf_n = nrb2->pf_v.pf_nt - nrb2->pf_v.pf_k;
  nrb_allocateknots(&nrb2->pf_v);
  mm = nrb2->pf_v.pf_nt;
  for (mm1 = 0; mm1 < mm; mm1++)
    nrb2->pf_v.pf_kk->knots[mm1] = nrb1->pf_u.pf_kk->knots[mm1];

  nnew = nrb2->pf_u.pf_n;
  nold = nrb1->pf_u.pf_n;

  /*Copy control points, transpose and leave space for new points*/

  mm = (nrb2->pf_u.pf_nt+tnew->pf_nt) * nrb2->pf_v.pf_nt;
  nrb_allocatepts(mm, &nrb2->pf_ppp);

  mm = nrb1->pf_u.pf_n;
  for (ii = 1; ii <= mm; ii++) {
    mm1 = nrb1->pf_v.pf_n;
    for (jj = 1; jj <= mm1; jj++)
      nrb2->pf_ppp->pts[ni(jj, ii, nnew) - 1] =
	nrb1->pf_ppp->pts[ni(ii, jj, nold) - 1];
  }

  tdir = &nrb2->pf_v;

  /*Now work with Y knots as though they are X!*/

  oldcpt = nrb1->pf_u.pf_n;
  oldkpt = nrb1->pf_u.pf_nt - 1;

  newcpt = oldcpt + tnew->pf_nt;
  newkpt = oldkpt + tnew->pf_nt;

  /*Add extra copy of last knot!*/

  tdir->pf_kk->knots[tdir->pf_nt - 1] = tdir->pf_kk->knots[nrb1->pf_u.pf_nt - 1];

  for (jj = tnew->pf_nt - 1; jj >= 0; jj--) {   /*jj*/
    mm1 = nrb1->pf_v.pf_nt;
    for (j = 1; j <= mm1; j++) {
      tpts1 = nrb2->pf_ppp;
      tpts1->pts[ni(j, newcpt, nnew) - 1] = tpts1->pts[ni(j, oldcpt, nnew) - 1];
    }
    while (tnew->pf_kk->knots[jj] <= tdir->pf_kk->knots[oldkpt - 1] &&
	   oldkpt > tdir->pf_k) {
      tdir->pf_kk->knots[newkpt - 1] = tdir->pf_kk->knots[oldkpt - 1];
      newkpt--;
      newcpt--;
      oldkpt--;
      oldcpt--;
      mm1 = nrb1->pf_v.pf_nt;
      for (j = 1; j <= mm1; j++) {
	tpts1 = nrb2->pf_ppp;
	tpts1->pts[ni(j, newcpt, nnew) - 1] = tpts1->pts[ni(j, oldcpt, nnew) - 1];
      }
    }
    mm1 = tdir->pf_k;
    for (l = 1; l < mm1; l++) {
      beta = tdir->pf_kk->knots[newkpt + l - 1] - tnew->pf_kk->knots[jj];
      if (beta <= 0) {
	mm2 = nrb1->pf_v.pf_nt;
	for (j = 1; j <= mm2; j++) {
	  tpts1 = nrb2->pf_ppp;
	  tpts1->pts[ni(j, newcpt + l - 1, nnew) - 1] =
	    tpts1->pts[ni(j, newcpt + l, nnew) - 1];
	}
      } else {
	beta /= tdir->pf_kk->knots[newkpt + l - 1] -
		tdir->pf_kk->knots[oldcpt + l - 1];
	mm2 = nrb1->pf_v.pf_nt;
	for (j = 1; j <= mm2; j++) {
	  tpts1 = nrb2->pf_ppp;
	  t1 = tp_scale4(tpts1->pts[ni(j, newcpt + l - 1, nnew) - 1], beta);
	  t2 = tp_scale4(tpts1->pts[ni(j, newcpt + l, nnew) - 1], 1.0 - beta);
	  tpts1->pts[ni(j, newcpt + l - 1, nnew) - 1] = tp_add4(t1, t2);
	}
      }
    }
    tdir->pf_kk->knots[newkpt - 1] = tnew->pf_kk->knots[jj];
    newkpt--;
    newcpt--;
  }
  /*with nrb2*/
}  /* nrb_boehmc*/
/*---------------------------------------------------------------------------*/
/*function:external*/
extern void nrb_bvalue(C(PR_nurb *)s1, C(Pfloat) x, C(Pint) jderiv, 
                       C(PR_nurb *)s2)
PreANSI(PR_nurb *s1) 
PreANSI(Pfloat x) 
PreANSI(Pint jderiv) 
PreANSI(PR_nurb *s2)
/*
NURB evaluation rountine taken from deBoors book. Given a NURB evaluate the NURB at somevalue x at some derivative jderiv.
*/
{

  /*from deBoors Book -REAL FUNCTION BVALUE ( T, BCOEF, N, K, X, JDERIV )
  c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  c
  C CALLS  INTERV
  c
  c calculates value at  x  of  jderiv-th derivative of spline from b-repr.
  c  the spline is taken to be continuous from the right.
  c
  C******  I N P U T  ******
  c  t, bcoef, n, k......forms the b-representation of the spline  f  to
  cbe evaluated. specifically,
  c  t.......knot sequence, of length  n+k, assumed nondecreasing.
  c  bcoef...b-coefficient sequence, of length  n .
  c  n.......length of  bcoef  and dimension of spline(k,t),
  ca s s u m e d positive .
  c  k.....order of the spline .
  c
  c  W A R N I N G . . .  the restriction k .le. kmax (=20)  is imposed
  carbitrarily by the dimension statement for  aj, dl, dr  below,
  cbut is NOWHERE CHECKED FOR.
  c
  c  x.....the point at which to evaluate.
  c  jderiv.....integer giving the order of the derivative to be evaluated
  ca s s u m e d to be zero or positive.
  c
  c****** o u t p u t ******
  c  bvalue.....the value of the (jderiv)-th derivative of f at x.
  c
  c****** m e t h o d ******
  cthe nontrivial knot interval (t(i),t(i+1)) containing x is lo-
  c  cated with the aid of interv.  the k b-coeffs of f relevant for
  c  this interval are then obtained from bcoef (or taken to be zero if
  c  not explicitly available) and are then differenced jderiv times to
  c  obtain the b-coeffs of (d**jderiv)f  relevant for that interval.
  c  precisely, with j = jderiv, we have from x.(12) of the text that
  c
  c (d**j)f = sum ( bcoef(.,j)*b(.,k-j,t) )
  c
  c  where
  c         /bcoef(.),, j .eq. 0
  c  /
  c   bcoef(.,j) = /bcoef(.,j-1) - bcoef(.-1,j-1)
  c /----------------------------- , j .gt. 0
  c                /   (t(k.+-j) - t(.))/(k-j)
  c
  c     then, we use repeatedly the fact that
  c
  c    sum ( a(.)*b(.,m,t)(x) ) = sum (a(.,x)*b(.,m-1,t)(x) )
  c  with
  c(x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
  c    a(.,x) =   ---------------------------------------
  c(x - t(.))      + (t(.+m-1) - x)
  c
  c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
  c  of order  1, and the coefficient for b(i,1,t)(x)  must then be the
  c  desired number (d**j)f(x). (see x.(17)-(19) of text).
  c
  c------------------------------------------------------------------
  */
  Pint i, mflag, km1, ll, imk, j, jc, ilow, kmj, jj, nin;
  Pfloat t1;
  Ppoint4 p1, p2;
  PR_pts *aj;
  PR_dir dl, dr;
  Pint mm;
  PR_dir *tdir;
  Pint mm1, mm2;
  Pint ajdim1;
  /*nrb_bvalue*/
  /*
     *** find  i  s.t.  1 .le. i .lt. n+k and t(i) .lt. t(i+1) and
          t(i) .le. x .lt. t(i+1) . if no such i can be found, x lies
          outside the support of the spline f and svalue = 0.
          (the asymmerty in this choice of i makes f rightcontinuous)
   WTH version of interv finds min value pf_K and max value pf_n
  */
  nrb_interv(&s1->pf_u, x, &i, &mflag);

  km1 = s1->pf_u.pf_k - 1;

  nrb_allocatepts(s1->pf_v.pf_n, &s2->pf_ppp);
  if (km1 <= 0) {
    mm = s1->pf_v.pf_n;
    for (ll = 1; ll <= mm; ll++)
      s2->pf_ppp->pts[ll - 1] = s1->pf_ppp->pts[ni(i, ll, s1->pf_u.pf_n) - 1];
  } else  /*km1 > 0  => k>1 - at least linear splines*/
  {   /*km1 <=0  => k>1 - not linear splines*/
    /*allocate points*/
    dl.pf_kk = NULL;
    dl.pf_nt = s1->pf_u.pf_nt;
    dl.pf_k = s1->pf_u.pf_k;   /*Used for consistency only*/
    dl.pf_n = s1->pf_u.pf_n;   /*Used for consistency only*/

    dr = dl;

    nin = s1->pf_u.pf_n;

    nrb_allocateknots(&dl);
    nrb_allocateknots(&dr);
    aj = NULL;
    ajdim1= P_max(s1->pf_u.pf_n,s1->pf_u.pf_k);
    nrb_allocatepts(s1->pf_v.pf_n * ajdim1, &aj);

    tdir = &s1->pf_u;

    /*c
    c   *** store the k b-spline coefficients relevant for the knot interval
    c(t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
    cdr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
    cfrom input to zero. set any t.s not obtainable equal to t(1) or
    cto t(n+k) appropriately.
    c
    */
    imk = i - tdir->pf_k;
    for (j = 1; j <= km1; j++)
      dl.pf_kk->knots[j - 1] = x - tdir->pf_kk->knots[i - j];
    for (j = 0; j < km1; j++)
      dr.pf_kk->knots[j] = tdir->pf_kk->knots[i + j] - x;

    mm = tdir->pf_k;
    for (jc = 1; jc <= mm; jc++) {
      mm1 = s1->pf_v.pf_n;
      for (ll = 1; ll <= mm1; ll++)
	aj->pts[ni(jc, ll, ajdim1) - 1] = s1->pf_ppp->pts[ni(imk + jc, ll, nin) - 1];
    }

    /*   *** difference the coefficients jderiv times.*/
    if (jderiv != 0) {   
      for (j = 1; j <= jderiv; j++) {
	kmj = tdir->pf_k - j;
	ilow = kmj;
	for (jj = 1; jj <= kmj; jj++) {
	  t1 = dl.pf_kk->knots[ilow - 1] + dr.pf_kk->knots[jj - 1];
	  mm2 = s1->pf_v.pf_n;
	  for (ll = 1; ll <= mm2; ll++) {
	    p1 = aj->pts[ni(jj + 1, ll, ajdim1) - 1];
	    p2 = aj->pts[ni(jj, ll, ajdim1) - 1];
	    p1 = tp_scale4(p1, kmj);
	    p2 = tp_scale4(p2, kmj / t1);
	    aj->pts[ni(jj, ll, ajdim1) - 1] = tp_sub4(p1, p2);
	  }
	  ilow--;
	}
      }
    }
    /* *** compute value at  x  in (t(i),t(i+1)) of jderiv-th derivative,
    cgiven its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).*/
    if (jderiv != km1) {
      for (j = jderiv + 1; j <= km1; j++) {
	kmj = s1->pf_u.pf_k - j;
	ilow = kmj;
	for (jj = 1; jj <= kmj; jj++) {
	  t1 = dl.pf_kk->knots[ilow - 1] + dr.pf_kk->knots[jj - 1];
	  mm2 = s1->pf_v.pf_n;
	  for (ll = 1; ll <= mm2; ll++) {
	    p1 = aj->pts[ni(jj + 1, ll, ajdim1) - 1];
	    p2 = aj->pts[ni(jj, ll, ajdim1) - 1];
	    p1 = tp_scale4(p1, dl.pf_kk->knots[ilow - 1] / t1);
	    p2 = tp_scale4(p2, dr.pf_kk->knots[jj - 1] / t1);
	    aj->pts[ni(jj, ll, ajdim1) - 1] = tp_add4(p1, p2);
	  }
	  ilow--;
	}
      }
    }

    mm = s1->pf_v.pf_n;
    for (ll = 1; ll <= mm; ll++)
      s2->pf_ppp->pts[ll - 1] = aj->pts[ni(1, ll, ajdim1) - 1];

    nrb_deallocatepts(&aj);
    nrb_deallocateknots(&dl);
    nrb_deallocateknots(&dr);
  }

  nrb_interchange(s2);

  /*Now force Y's to be 1D!*/
  s2->pf_v.pf_k = 0;
  s2->pf_v.pf_n = 1;
  s2->pf_v.pf_nt = 1;
  nrb_allocateknots(&s2->pf_v);
  s2->pf_v.pf_kk->knots[0] = 0.0;

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