!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION SPLINE(X,Y,XIN) RESULT(YOUT)
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
!CC   INPUT: X  x-values with X(i) < X(i+1)
!CC          Y  y-values
!CC          XIN  x-value for which y will be returned
!CC --------------------------------------------------------------
!CC   WARNING: will crash if X(i) = X(i+1); exeption not caught
!CC --------------------------------------------------------------
!CC   Function returns y value for a corresponding x value, 
!CC   based on cubic spline.
!CC   Will never oscillates or overshoot. No need to solve matrix.
!CC   Also calculate constants for cubic in case needed (for integration).
!CC --------------------------------------------------------------
!CC   adopted from VBasic routine 
!CC   Developer: C Kruger, Guildford, UK     
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
      REAL(KIND=minuitPrec),DIMENSION(:),INTENT(IN) :: X,Y
      REAL(KIND=minuitPrec),INTENT(IN) :: XIN
      REAL(KIND=minuitPrec) :: YOUT,A,B,C,D,FPI,FPIM1,F2I,F2IM1
      INTEGER :: NPOINTS,LEG,I
      
      NPOINTS = SIZE(X)
      
      IF( NPOINTS <= 0 )THEN
         YOUT = 0.0
         RETURN
      ELSEIF( NPOINTS == 1 )THEN
         YOUT = Y(1)
         RETURN
      ENDIF
      
! 1) Find LineNumber or segment/leg. Linear extrapolate if outside range.
      IF (XIN >= X(NPOINTS)) THEN
!         B=( Y(NPOINTS)-Y(NPOINTS-1) )/( X(NPOINTS)-X(NPOINTS-1) )
!         A= Y(NPOINTS) - B*X(NPOINTS)
!
!         YOUT = A + B*XIN
         YOUT = Y(NPOINTS)
         RETURN
      ELSEIF (XIN <= X(1)) THEN
!         B=( Y(2)-Y(1) )/( X(2)-X(1) )
!         A= Y(2) - B*X(2)

!         YOUT = A + B*XIN
         YOUT= Y(1)
         RETURN
      ELSE
         LEG = NPOINTS
         DO I=2,NPOINTS
            IF (XIN <= X(I)) THEN
               LEG=I
               EXIT
            ENDIF
         ENDDO
      ENDIF
      
! 2) Calc first derivative (slope) for intermediate points
! 3) Reset first derivative (slope) at first and last point
      IF (LEG == 2) THEN
!        First point has 0 2nd derivative  
         FPI   = FPRIME(X(LEG-1:LEG+1),Y(LEG-1:LEG+1))    
         FPIM1 = (3.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) - FPI)/2.0
      ELSEIF (LEG == NPOINTS) THEN
!        Last point has 0 2nd derivative
         FPIM1 = FPRIME(X(LEG-2:LEG),Y(LEG-2:LEG))         
         FPI = (3.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) -  FPIM1)/2.0
      ELSE
         FPIM1 = FPRIME(X(LEG-2:LEG  ),Y(LEG-2:LEG  ))
         FPI   = FPRIME(X(LEG-1:LEG+1),Y(LEG-1:LEG+1))
      ENDIF


! 4) Calc second derivative at points
      F2IM1 = -2.0*( FPI + 2.0*FPIM1 )/( X(LEG)-X(LEG-1) ) &
             + 6.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2
 
      F2I   =  2.0*( 2.0*FPI + FPIM1 )/( X(LEG)-X(LEG-1) ) &
             - 6.0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2


! 5) Calc constants for cubic
      D = ( F2I-F2IM1 )/( 6.0*(X(LEG)-X(LEG-1)) ) 
      C = ( X(LEG)*F2IM1 - X(LEG-1)*F2I )/( 2.0*( X(LEG)-X(LEG-1 )) )
      B = ( Y(LEG)-Y(LEG-1) - C*( X(LEG)**2-X(LEG-1)**2 ) - D*( X(LEG)**3-X(LEG-1)**3 ) ) & 
         / ( X(LEG)-X(LEG-1) )
      A = Y(LEG-1) - B*X(LEG-1) - C*X(LEG-1)**2 - D*X(LEG-1)**3

      YOUT = A + B*XIN + C*XIN**2 + D*XIN**3
      END FUNCTION SPLINE
      FUNCTION FPRIME(X,Y) RESULT(FP)
      REAL(KIND=minuitPrec),DIMENSION(3),INTENT(IN) :: X,Y         
      REAL(KIND=minuitPrec) :: FP,DX1,DX2,DY1,DY2 
      
      DX1 = X(2)-X(1)
      DX2 = X(3)-X(2)
      DY1 = Y(2)-Y(1)
      DY2 = Y(3)-Y(2)
      
      
      IF( DY1==0.0 .OR. DY2==0.0 )THEN
! Prevent divide by 0
         FP = 0.0
      ELSEIF( DY1*DY2 < 0.0 )THEN
! Pos AND neg slope, assume slope = 0 to prevent overshoot
         FP = 0.0         
      ELSE
         FP =  2.0/( DX1/DY1 + DX2/DY2 )
      ENDIF
      END FUNCTION FPRIME