Pfloat bvalue(Pfloat *t, Pfloat *bcoef, Pint n, Pint k, Pfloat x)
{

  /*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, ilo, kmj, jj, nin;
  Pfloat t1;
  Pfloat  p1, p2;
  Pfloat aj[20], dr[20],dl[20];
  Pint mm,jcmin,jcmax,nmi;
  Pint mm1, mm2;

   
  /*nrb_bvalue*/
  /*printf("Bvalue: n:%3ld, k:%3ld, x:%12.3f\nKnots:",n,k,x);
  for (j=1;j<=n+k; j++)printf("%8.3f",t[j-1]);
  printf("\nBcoef:");
  for (j=1;j<=n; j++)printf("%8.3f",bcoef[j-1]);
  printf("\n");
  */
  /*
     *** 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
  */
  interv(t,n+k,x, &i, &mflag);

  km1 = k - 1;

  if (km1 <= 0) return (Pfloat) bcoef[i-1];

  /*km1 <=0  => k>1 - not linear splines*/
  /*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 - k;
    jcmin = 1;
    if (imk < 0) {
      jcmin = i - imk;
      for (j = 1; j <= i; j++)   dl[j - 1] = x - t[i - j];
      for (j = i; j <= km1; j++){
        aj[k-j-1] = 0.0;
        dl[j-1] = dl[i-1];
      }
    } else
    for (j=1; j<= km1; j++)  dl[j] = x - t[i - j];

    jcmax = k;
    nmi = n - i;
    if (nmi < 0){
      jcmax = k + nmi;
       for (j = 1; j <= jcmax; j++) dr[j-1] = t[i+j-1] -x;
       for ( j= jcmax; j<=km1; j++){
       	 aj[j] =0.0;
       	 dr[j-1] = dr[jcmax-1];
       }
    }
    else for (j=1; j<=km1;j++) dr[j-1] =t[i+j-1] -x;

    for (jc =jcmin;jc<=jcmax;jc++)aj[jc-1] = bcoef[imk+jc-1];
    
    /* *** 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).*/

      for (j = 1; j <= km1; j++) {
	kmj =  - j;
	ilo = kmj;
	for (jj = 1; jj <= kmj; jj++) {
	  aj[jj-1] =
          (aj[jj]*dl[ilo-1] +aj[jj-1]*dr[jj-1])/(dl[ilo-1]+dr[jj-1]);
	  ilo--;
	}
      }
     return (Pfloat) aj[0];

}

