!   original from:
!   C. de Boor. A Practical Guide To Splines. Chapter X.
!   REWRITTEN by V.Baturin, 2015
!   for purpose of not-a-knot interpolation with MATLAB
!
!   Values of all (needed) B-splines, nonzero in point X, with order Jhigh
!        simplified in respect of Index
!   knots of B-spline in array T
!
! --Input--------------------------------------------------
!
!   T - knots, length = Left+Jhigh, non-decreasing order
!
!   x - point where B-spline is computed
!
!   Left - integer, T(Left) .LE. T(Left+1)
!
! --Output-------------------------------------------------
!
!   BIATX - array with Jout length, BIATX(i) contains value
!           of polynom order Jout in point x ( B-spline
!           B(Left-Jout+i), Jout, T) on interval ( T(Left),
!           T(Left+1)
!
      SUBROUTINE BSPLVB8v( T, Jhigh, x, LC, BIATX )

      USE mod_numres
!  Variables declaration

      IMPLICIT NONE
      
      INTEGER, INTENT(in) ::  Jhigh
      REAL*8, INTENT(in) :: T(:), x
      INTEGER, INTENT(out) ::  LC
      REAL*8, INTENT(out) :: BIATX(Jhigh)
! locals
      INTEGER Left, ki, kj, jp1
      REAL*8 DeltaL(Jhigh), DeltaR(Jhigh), Saved, Term
      INTEGER lenT, mflag
      
      kj = 1 ! order of curently calculated spline (simplification)
      lenT=SIZE(T)
!  calculation of Left (asymmetric on right edge)
      CALL interv ( T, lenT, x, Left, mflag )    
!
       LC=Left-Jhigh+1        ! LC - index in coef-array. 
!
      BIATX(1) = 1 ! value for spline of first order (constant)
      DO WHILE ( kj < Jhigh )
          jp1 = kj + 1
          DeltaR(kj) = T(Left+kj) - x
          DeltaL(kj) = x - T(Left+1-kj)
          Saved = 0.d0
!
          DO ki = 1, kj
             Term = BIATX(ki) / ( DeltaR(ki) + DeltaL(jp1-ki) )
             BIATX(ki) = Saved + DeltaR(ki)*Term
             Saved = DeltaL(jp1-ki)*Term
          ENDDO
!
         BIATX(jp1) = Saved
         kj = jp1
      ENDDO
      END SUBROUTINE BSPLVB8v