FUNCTION BSKFUN(X,K,T,NKN,C)
*     Returned is the function value at X of B-spline function of 
*     order K within the knot range.
*     C(*) is the array of coefficients, and T(NKN) is the knot array.
      PARAMETER (KMAX=10) 
      REAL T(NKN),B(KMAX),C(*)
      DOUBLE PRECISION SUM
*     ...
      L=MAX(K,MIN(LEFT(X,T,NKN),NKN-K))   ! index within limits
      CALL BSPLNK(X,K,B,L,T)                ! calculate K B-splines 
      SUM=B(1)*C(L-K+1)
      DO I=2,K
       SUM=SUM+B(I)*C(I+L-K) 
      END DO 
      BSKFUN=SUM
      END