Chimie Théorique » Opt'n Path

       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