/*
 * *****************************************************************
 *
 * ROUTINE:  bspl_eval()
 *
 * FUNCTION: Evaluate a B-spline basis function
 *           evaluates N^n_i(u) over knot sequence
 *	     u_0,....,u_{l+2n-2}
 *
 * INPUT:  n .............. degree
 *         l .............. no. of domain intervals
 *         i .............. interval in which to evaluate
 *         knot ........... knot sequence
 *         u .............. parameter value at which to evaluate
 *                          u must be in the domain intervals. 
 *                          This is not checked.
 *
 * OUTPUT: bspl_eval ...... basis function value at input u
 *
 * GLOBAL: none 
 * 
 * METHOD: Use recursive definition of B-spline basis functions.
 *         This is a recursive routine.
 *         Special care was taken when a basis function was not
 *         defined to set it to zero.
 *
 * CALLS:  bspl_eval() ..... itself
 *
 * AUTHOR: GF&DCH
 *
 * CHANGES:
 * 
 *
 * *****************************************************************
 */

#include <stdio.h>
#include <math.h>
#include <string.h>


/* extern FILE *infofile; */


double bspl_eval(); 

double bspl_eval(int    n, 
                 int    num_knots, 
                 int    i, 
                 double knot[], 
                 double u,
                 int    last_knot)

{
  int    n1;
  int    i1;
  int    left;
  int    right;
  double value;
  double num1;
  double num2;
  double denom1;
  double denom2;
  double left_val;
  double right_val;
  double zero_eps = 0.00000001;

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

  /*
  fprintf(infofile,"BSPL: n=%d i=%d num_knots=%d last_knot=%d  u=%f \n",
          n,i,num_knots,last_knot,u);
  fflush(infofile);
  */

  if(i > 0) {
    if(u < knot[i-1]) return(0.0);
  }
  if(i+n < num_knots) {
    if(u > knot[i+n]) return(0.0);
  }


  n1 = n-1; 
  i1 = i+1;


  if (n == 0) {
    if(i == last_knot) {
      if(knot[i-1]  <= u && u <= (knot[i]+0.000005)) {
        /* fprintf(infofile,"BSPL: special return 1\n"); */
        return(1.0);
      }
    }
    else  {
      if(knot[i-1]  <= u && u < knot[i]) {
        /* fprintf(infofile,"BSPL: return 1\n"); */
        return(1.0);
      }
    }
    /* fprintf(infofile,"BSPL: return 0\n"); */
    return (0.0);
  }

  left = 1;
  right = 1;

  if(i-1 < 0) left = 0;
  if(i+n >= num_knots) right = 0;

  /* fprintf(infofile,"BSPL: index check left=%d right=%d\n",left,right); */

  if(left) {
    denom1 = knot[i+n-1]-knot[i-1];
    if(denom1 < zero_eps)  left = 0;
  }

  if(right) {
    denom2 = knot[i+n]-knot[i];
    if(denom2 < zero_eps) right = 0;
  }

  /* fprintf(infofile,"BSPL: denom check left=%d right=%d\n",left,right); */


  if(left)  num1 = (u - knot[i-1])/denom1;
  if(right) num2 = ( knot[i+n] - u)/denom2;

    
  if(left)  left_val =  num1 * bspl_eval(n1, num_knots, i, knot, u, last_knot);
  else      left_val = 0.0;

  if(right) right_val = num2 * bspl_eval(n1, num_knots, i1, knot, u, last_knot);
  else      right_val = 0.0;
   
  /* 
  fprintf(infofile,"BSPL: n=%d i=%d  left_val=%f right_val=%f\n",
          n,i,left_val,right_val); 
  */

  value =  left_val + right_val;             
  return(value);	
}