       SUBROUTINE spline(x,y,n,yp1,ypn,y2)
!Given arrays x(1:n) and y(1:n) containing a tabulated function,
! i.e., y i = f(xi), with x1<x2< :::<xN , and given values yp1 and ypn
!for the  rst derivative of the inter- polating function at points 1 and n,
! respectively, this routine returns an array y2(1:n) of length n
!which contains the second derivatives of the interpolating function
!at the tabulated points xi.
! Ifyp1 and/or ypn are equal to 1*10^30 or larger,
! the routine is signaled to set the corresponding boundary
!condition for a natural spline,
!with zero second derivative on that boundary.

!----------------------------------------------------------------------
!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 
!  Centre National de la Recherche Scientifique,
!  Université Claude Bernard Lyon 1. All rights reserved.
!
!  This work is registered with the Agency for the Protection of Programs 
!  as IDDN.FR.001.100009.000.S.P.2014.000.30625
!
!  Authors: P. Fleurat-Lessard, P. Dayal
!  Contact: optnpath@gmail.com
!
! This file is part of "Opt'n Path".
!
!  "Opt'n Path" is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Affero General Public License as
!  published by the Free Software Foundation, either version 3 of the License,
!  or (at your option) any later version.
!
!  "Opt'n Path" is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU Affero General Public License for more details.
!
!  You should have received a copy of the GNU Affero General Public License
!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.
!
! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,
! for commercial licensing opportunities.
!----------------------------------------------------------------------

         Use Vartypes

         IMPLICIT NONE

! Number of points
         INTEGER(KINT) :: n
! Coordinate to spline
         REAL(KREAL) :: x(n),y(n)
! Coefficients to compute
         REAL(KREAL) :: y2(n)
! End-points derivatives
         REAL(KREAL) ::yp1,ypn

       INTEGER(KINT), PARAMETER  :: NMAX=500

       INTEGER(KINT) :: i,k
       REAL(KREAL) :: p,qn,sig,un,u(NMAX)
       LOGICAL Debug

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

  END INTERFACE

       Debug=valid("spline")

       IF (DEBUG) THEN
          WRITE(*,*) "Spline 1D",n
          WRITE(*,*) "x:",(x(i),i=1,n)
          WRITE(*,*) "y:",(y(i),i=1,n)
       END IF

       if (yp1.gt..99e30) then
! The lower boundary condition is set either to be \natural"
        y2(1)=0.
        u(1)=0.
       else
! or else to have a speci ed  rst derivative.
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do i=2,n-1
! This is the decomposition loop of the tridiagonal algorithm.
!y2 and u are used for temporary storage of the decomposed factors.
       sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
       y2(i)=(sig-1.)/p
       u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
                /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
       enddo
       if (ypn.gt..99e30) then
!The upper boundary condition is set either to be \natural"
         qn=0.
         un=0.
        else
!or else to have a speci ed  rst derivative.
         qn=0.5
         un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
         endif
        y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
        do k=n-1,1,-1
         y2(k)=y2(k)*y2(k+1)+u(k)
        enddo

        IF (DEBUG)  WRITE(*,*) "y2:",(y2(i),i=1,n)
         return
        END


      SUBROUTINE splint(x,y,N,xa,ya,y2a)

         Use Vartypes

         IMPLICIT NONE

! Number of points
         INTEGER(KINT) :: n
! Spline coefficients
         REAL(KREAL) :: xa(n),ya(n), y2a(n)
! Y to compute for a given x
         REAL(KREAL) :: x,y


       LOGICAL debug

! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
! function (with the xai's in order), and given the array y2a(1:n),
!which is the output from spline above, and given a value of x,
! this routine returns a cubic-spline interpolated value y.
       INTEGER(KINT) :: k,khi,klo
       REAL(KREAL) :: a,b,h


  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

  END INTERFACE

       Debug=valid("splint")


!           if (debug) THEN
!       WRITE(*,*) "Splint 1D",n
!       WRITE(*,*) "xa:",(xa(i),i=1,n)
!       WRITE(*,*) "ya:",(ya(i),i=1,n)
!       WRITE(*,*) "y2a:",(y2a(i),i=1,n)
!        END IF



          klo=1
!We will  find the right place in the table by means of bisection.
! This is optimal if sequential calls to this routine are at
!random values of x. If sequential calls are in order, and closely
! spaced, one would do better to store previous values of klo and khi
! and test if they remain appropriate on the next call.
          khi=n
 !         WRITE(*,*) xa(klo),xa(khi),n,x
 1        if (khi-klo.gt.1) then
             k=(khi+klo)/2
             if(xa(k).gt.x) then
                khi=k
             else
                klo=k
             endif
          goto 1
          endif
! klo and khi now bracket the input value of x.
          h=xa(khi)-xa(klo)
! The xa's must be distinct.
          if (h.eq.0.) THEN
             WRITE(*,*) 'bad xa input in splint'
             STOP
          END IF
! Cubic spline polynomial is now evaluated.
           a=(xa(khi)-x)/h
           b=(x-xa(klo))/h
           y=a*ya(klo)+b*ya(khi)+ ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi)) &
                *(h**2)/6.

!           WRITE(*,*) "Splint1D x,y",x,y
          return
          END


! SUBROUTINE splinder *************************************************
! *********************************************************************

      SUBROUTINE splinder(x,y,N,xa,ya,y2a)
         Use Vartypes

         IMPLICIT NONE

! Number of points
         INTEGER(KINT) :: n
! Spline coefficients
         REAL(KREAL) :: xa(n),ya(n), y2a(n)
! Y to compute for a given x
         REAL(KREAL) :: x,y


       LOGICAL debug

! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
! function (with the xai's in order), and given the array y2a(1:n),
! which is the output from spline above, and given a value of x,
! this routine returns a cubic-spline interpolated value y first DERIVATIVE.
! this routine returns a cubic-spline interpolated value y.
       INTEGER(KINT) :: k,khi,klo
       REAL(KREAL) :: a,b,h


  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

  END INTERFACE

 
       Debug=valid("splinder")



          klo=1
! We will  find the right place in the table by means of bisection.
! This is optimal if sequential calls to this routine are at
!random values of x. If sequential calls are in order, and closely
! spaced, one would do better to store previous values of klo and khi
! and test if they remain appropriate on the next call.
          khi=n
 !         WRITE(*,*) xa(klo),xa(khi),n,x
 1        if (khi-klo.gt.1) then
             k=(khi+klo)/2
             if(xa(k).gt.x) then
                khi=k
             else
                klo=k
             endif
          goto 1
          endif
! klo and khi now bracket the input value of x.
          h=xa(khi)-xa(klo)
! The xa's must be distinct.
          if (h.eq.0.) THEN
             WRITE(*,*) 'bad xa input in splinder'
             STOP
          END IF
! Cubic spline polynomial is now evaluated.
           a=(xa(khi)-x)/h
           b=(x-xa(klo))/h
! Formula taken from the Numerical Recipies book.
           y=(ya(khi)-ya(klo))/h - (3*a**2-1)/6.*h*y2a(klo)+ &
                (3*b**2-1)/6.*h*y2a(khi)
          return
          END



! SUBROUTINE splintder *************************************************
! *********************************************************************

      SUBROUTINE splintDer(x,y,der,N,xa,ya,y2a)

! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
! function (with the xai's in order), and given the array y2a(1:n),
! which is the output from spline above, and given a value of x,
! this routine returns a cubic-spline interpolated value y.
! and the derivative der.

         Use Vartypes

         IMPLICIT NONE

! Number of points
         INTEGER(KINT) :: n
! Spline coefficients
         REAL(KREAL) :: xa(n),ya(n), y2a(n)
! Y to compute for a given x
         REAL(KREAL) :: x,y, der


       LOGICAL debug

       INTEGER(KINT) :: k,khi,klo
       REAL(KREAL) :: a,b,h

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

  END INTERFACE

       Debug=valid("splintder")


!           if (debug) THEN
!       WRITE(*,*) "SplintDer 1D",n
!       WRITE(*,*) "xa:",(xa(i),i=1,n)
!       WRITE(*,*) "ya:",(ya(i),i=1,n)
!       WRITE(*,*) "y2a:",(y2a(i),i=1,n)
!        END IF



          klo=1
! We will  find the right place in the table by means of bisection.
! This is optimal if sequential calls to this routine are at
! random values of x. If sequential calls are in order, and closely
! spaced, one would do better to store previous values of klo and khi
! and test if they remain appropriate on the next call.
          khi=n
 !         WRITE(*,*) xa(klo),xa(khi),n,x
 1        if (khi-klo.gt.1) then
             k=(khi+klo)/2
             if(xa(k).gt.x) then
                khi=k
             else
                klo=k
             endif
          goto 1
          endif
! klo and khi now bracket the input value of x.
          h=xa(khi)-xa(klo)
! The xa's must be distinct.
          if (h.eq.0.) THEN
             WRITE(*,*) 'bad xa input in splintDer'
             STOP
          END IF

! Cubic spline polynomial is now evaluated.
           a=(xa(khi)-x)/h
           b=(x-xa(klo))/h
           y=a*ya(klo)+b*ya(khi)+ ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi)) &
                *(h**2)/6.
! Formula taken from the Numerical Recipies book.
           der=(ya(khi)-ya(klo))/h - (3*a**2-1)/6.*h*y2a(klo)+ &
                (3*b**2-1)/6.*h*y2a(khi)
          return
        END SUBROUTINE splintDer



      SUBROUTINE LinearInt(x,y,der,N,xa,ya)

! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
! function (with the xai's in order),  and given a value of x,
! this routine returns a linear interpolated value y and its derivative der

         Use Vartypes

         IMPLICIT NONE

! Number of points
         INTEGER(KINT) :: n
! Spline coefficients
         REAL(KREAL) :: xa(n),ya(n)
! Y to compute for a given x
         REAL(KREAL) :: x,y, der


       LOGICAL debug

       INTEGER(KINT) :: k,khi,klo
       REAL(KREAL) :: a,b,h

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

  END INTERFACE

       Debug=valid("linearint")


!           if (debug) THEN
!       WRITE(*,*) "SplintDer 1D",n
!       WRITE(*,*) "xa:",(xa(i),i=1,n)
!       WRITE(*,*) "ya:",(ya(i),i=1,n)
!       WRITE(*,*) "y2a:",(y2a(i),i=1,n)
!        END IF



          klo=1
! We will  find the right place in the table by means of bisection.
! This is optimal if sequential calls to this routine are at
! random values of x. If sequential calls are in order, and closely
! spaced, one would do better to store previous values of klo and khi
! and test if they remain appropriate on the next call.
          khi=n
 !         WRITE(*,*) xa(klo),xa(khi),n,x
 1        if (khi-klo.gt.1) then
             k=(khi+klo)/2
             if(xa(k).gt.x) then
                khi=k
             else
                klo=k
             endif
          goto 1
          endif
! klo and khi now bracket the input value of x.
          h=xa(khi)-xa(klo)
! The xa's must be distinct.
          if (h.eq.0.) THEN
             WRITE(*,*) 'bad xa input in LinearIntDer'
             STOP
          END IF

! Linear int now evaluated.
           a=(xa(khi)-x)/h
           b=(x-xa(klo))/h
           y=a*ya(klo)+b*ya(khi)
! Formula taken from the Numerical Recipies book.
           der=a+b
          return
        END SUBROUTINE LinearInt
