#if defined HAVE_CONFIG_H #include "config.h" #endif !!@LICENSE ! ! ============================================================================== ! MODULE interpolation ! ============================================================================== ! Cubic spline interpolation utilities ! ============================================================================== ! Public types in this module: ! spline_t : derived type to hold spline info of a function ! Public parameters, variables and arrays in this module: ! none ! Public procedures provided by this module: ! generate_spline : generate spline interpolation info of a function ! evaluate_spline : evaluate a function at given point(s) using splines ! ============================================================================== ! SUBROUTINE generate_spline( dat, x, y, n, dydx1, dydxn, d2ydx2, stat ) ! Generate data required to interpolate a 1D function with cubic splines ! Required input: ! real(dp) x(n) : independent variable at each point ! real(dp) y(n) : function values at each point ! integer n : number of points ! Optional input: ! real(dp) dydx1 : dy/dx at x(1) (if not present, assumes d2y/dx2=0 at x(1)) ! real(dp) dydxn : dy/dx at x(n) (if not present, assumes d2y/dx2=0 at x(n)) ! Output: ! type(spline_t) dat : spline info of function y(x) ! Optional output: ! real(dp) d2ydx2(n) : d2y/dx2 at mesh points ! integer stat : error status ! Behaviour: ! If dydx1 and/or dydxn are not present, assumes d2y/dx2=0 at x(1) and/or x(n) ! Returns with stat=-1 if mesh is not monotonic, without any other output ! Returns with stat=-2 if n<2, without any other output ! Algorithm: ! Computes d2y/dx2 at each point, to make y(x) and dy/dx continuous. ! The x values are analyzed to check if the mesh is linear, logarithmic, ! or general. This is used by evaluate_spline to find the mesh interval at ! which the evaluation points lie. ! Ref: "Numerical Recipes", W.H. Press et al, Cambridge U.P. ! ============================================================================== ! SUBROUTINE evaluate_spline( dat, x, y ) ! Evaluate function at given point(s) using spline data info ! Input: ! type(spline_t) dat : spline info of function y(x), from generate_spline ! real(dp) x(:) : point(s) at which function must be evaluated ! Output: ! real(dp) y(:) : value of function at given point(s) ! Behaviour: ! Single- and multiple-point evaluators are overloaded with the same name, ! so that x and y may be also scalar values ! Stops with an error message if any point x is out of the interp. interval ! ============================================================================== ! subroutine clean_spline(dat) ! type(spline_t), intent(inout) :: dat ! Behavior: ! Deallocates the array components and resets %n to zero ! end subroutine clean_spline ! ============================================================================== ! SUBROUTINE spline( x, y, n, dydx1, dydxn, d2ydx2 ) ! Included for compatibility with Numerical Recipes interface ! Input: ! real(dp) x(n) ! mesh points ! real(dp) y(n) ! function value at mesh points ! integer n ! number of mesh points ! real(dp) dydx1 ! dy/dx at x(1). If >1.e30, assumes d2y/dx2=0 at x(1) ! real(dp) dydxn ! dy/dx at x(n). If >1.e30, assumes d2y/dx2=0 at x(n) ! Output: ! real(dp) d2ydx2(n) ! d2y/dx2 at mesh points ! ============================================================================== ! SUBROUTINE spline( dx, y, n, dydx1, dydxn, d2ydx2 ) ! Included for compatibility with interface used in siesta ! Input: ! real(dp) dx ! mesh interval of a uniform mesh starting at x=0 ! real(dp) y(n) ! function value at mesh points ! integer n ! number of mesh points ! real(dp) dydx1 ! dy/dx at x(1). If >1.e30, assumes d2y/dx2=0 at x(1) ! real(dp) dydxn ! dy/dx at x(n). If >1.e30, assumes d2y/dx2=0 at x(n) ! Output: ! real(dp) d2ydx2(n) ! d2y/dx2 at mesh points ! ============================================================================== ! SUBROUTINE splint( xi, yi, d2ydx2, n, x, y, dydx ) ! Included for compatibility with Numerical Recipes interface ! Input: ! real(dp) xi(n) ! mesh points ! real(dp) yi(n) ! function value at mesh points ! real(dp) d2ydx2(n) ! second derivative at mesh points ! integer n ! number of mesh points ! real(dp) x ! point at which function is needed ! Output: ! real(dp) y ! function value at point x ! real(dp) dydx ! function derivative at point x ! ============================================================================== ! SUBROUTINE splint( dx, yi, d2ydx2, n, x, y, dydx ) ! Included for compatibility with interface used in siesta ! Input: ! real(dp) dx ! mesh interval of a uniform mesh starting at x=0 ! real(dp) yi(n) ! function value at mesh points ! real(dp) d2ydx2(n) ! second derivative at mesh points ! integer n ! number of mesh points ! real(dp) x ! point at which function is needed ! Output: ! real(dp) y ! function value at point x ! real(dp) dydx ! function derivative at point x ! ============================================================================== ! SUBROUTINE polint(XA,YA,N,X,Y,DYDX) ! Lagrange interpolation ! Input: ! real*8 XA(N) : x values of the function y(x) to interpolate ! real*8 YA(N) : y values of the function y(x) to interpolate ! integer N : Number of data points ! real*8 X : x value at which the interpolation is desired ! Output: ! real*8 Y : interpolated value of y(x) at X ! real*8 DYDX : interpolated derivative dy/dx at X ! Notice: this argument has a different meaning ! in the Numerical Recipes' polint subroutine ! Ref: ! W.H.Press et al, Numerical Recipes, Cambridge U.P. ! ============================================================================== ! Written by J.M.Soler and A.Garcia. May.2015 ! ============================================================================== ! ! This module calls an external subroutine 'die', with interface ! ! interface ! subroutine die(str) ! character(len=*), intent(in) :: str ! end subroutine die ! end interface MODULE gridxc_interpolation implicit none integer, parameter :: dp = selected_real_kind(10,100) PUBLIC :: & spline_t, &! derived type to hold spline info of a function generate_spline, &! generate spline info evaluate_spline, &! evaluate a function at given point(s) clean_spline, &! deallocate the components of a spline_t object spline, &! compatible with Numerical Recipes interface splint, &! compatible with Numerical Recipes interface polint ! Lagrange polynomial interpolation PRIVATE ! nothing is declared public beyond this point ! Derived type to hold spline info of a function type spline_t private character(len=3):: mesh ! mesh type ('lin'|'log'|'gen') real(dp) :: xmin ! x(1)-xtol real(dp) :: xmax ! x(n)+xtol real(dp) :: x1 ! x(1) real(dp) :: dx ! linear mesh interval real(dp) :: a,b ! log-mesh parameters x(k)=b*[exp(a*(k-1))-1] integer :: n ! number of points real(dp),allocatable:: x(:) ! mesh points real(dp),allocatable:: y(:) ! function value at mesh points real(dp),allocatable:: d2ydx2(:) ! 2nd derivative at mesh points end type ! Overloaded spline generators and evaluators interface generate_spline module procedure generate_spline_master ! new interface module procedure generate_spline_x ! Numerical Recipes interface module procedure generate_spline_dx ! older interface end interface generate_spline interface evaluate_spline module procedure evaluate_spline ! evaluate function at single point module procedure evaluate_spline_n ! evaluate function at multiple points module procedure evaluate_spline_x ! Numerical Recipes interface module procedure evaluate_spline_dx ! siesta interface end interface evaluate_spline interface spline module procedure generate_spline_x ! Numerical Recipes interface module procedure generate_spline_dx ! siesta interface end interface spline interface splint module procedure evaluate_spline_x ! Numerical Recipes interface module procedure evaluate_spline_dx ! siesta interface end interface splint ! Internal module parameters real(dp),parameter:: dydxMax = 0.99e30_dp ! max. dy/dx at x(1) and x(n) real(dp),parameter:: tol = 1.e-6_dp ! tolerance for 'within range' condition CONTAINS !------------------------------------------------------------------------------- SUBROUTINE generate_spline_master( dat, x, y, n, dydx1, dydxn, d2ydx2, store, stat ) ! Generate data for cubic spline interpolation of a function implicit none type(spline_t), intent(out):: dat ! spline interpolation info integer, intent(in) :: n ! number of mesh points real(dp), intent(in) :: x(n) ! independent variable at mesh points real(dp), intent(in) :: y(n) ! function value at mesh points real(dp),optional,intent(in) :: dydx1 ! dy/dx at x(1) real(dp),optional,intent(in) :: dydxn ! dy/dx at x(n) real(dp),optional,intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points logical, optional, intent(in):: store ! Store data in the spline object integer, optional,intent(out):: stat ! error status: ! (-1 if nonmonotonic mesh) ! (-2 if n < 2) ! Internal variables and arrays character(len=*),parameter:: myName = 'generate_spline_master ' real(dp):: a, b, dx, dx1, dx2, dxn, dxm, dxp, dy1, dyn, dym, dyp, & p, s, u(n), v(n), xtol, ypp(n) integer :: flag, k character(len=3):: meshType ! Check that mesh is monotonic if (n<2) then flag = -2 ! single-point elseif (n>2 .and. any( (x(2:n-1)-x(1:n-2)) * (x(3:n)-x(2:n-1)) <= 0.0_dp )) then flag = -1 ! non-monotonic mesh else flag = 0 endif if (present(stat)) stat = flag if (flag<0) return ! Find if mesh is linear (a=0) or logarithmic ! x(k)=x(1)+b*[exp(a*(k-1))-1] => (x(k+1)-x(k))/(x(k)-x(k-1))=exp(a) if (n<3) then meshType = 'lin' dx = x(2)-x(1) else a = log( (x(3)-x(2)) / (x(2)-x(1)) ) if (all(abs( (x(3:n)-x(2:n-1))/(x(2:n-1)-x(1:n-2))-exp(a) ) < 1.e-12_dp)) then if (abs(a)<1.e-12_dp) then meshType = 'lin' dx = x(2)-x(1) else ! try x(k) = x(1) + b*[exp(a*(k-1))-1] meshType = 'log' b = (x(2)-x(1)) / (exp(a)-1._dp) endif else meshType = 'gen' endif endif ! Set boundary conditions at end points if (present(dydx1)) then dx1 = x(2)-x(1) dy1 = y(2)-y(1) v(1) = -0.5_dp u(1) = (3.0_dp/dx1)*(dy1/dx1-dydx1) else v(1) = 0._dp u(1) = 0._dp endif if (present(dydxn)) then dxn = x(n)-x(n-1) dyn = y(n)-y(n-1) v(n) = 0.5_dp u(n) = (3.0_dp/dxn)*(dydxn-dyn/dxn) else v(n) = 0.0_dp u(n) = 0.0_dp endif ! Decomposition loop of the tridiagonal equations do k = 2,n-1 dxm = x(k)-x(k-1) dym = y(k)-y(k-1) dxp = x(k+1)-x(k) dyp = y(k+1)-y(k) s = dxm / (dxm+dxp) p = s*v(k-1) + 2.0_dp v(k) = (s-1.0_dp)/p u(k) = ( 6.0_dp*(dyp/dxp-dym/dxm)/(dxm+dxp) - s*u(k-1) )/p enddo ! Backsubstitution loop of the tridiagonal equations ypp(n) = (u(n)-v(n)*u(n-1)) / (v(n)*v(n-1)+1.0_dp) do k = n-1,1,-1 ypp(k) = v(k)*ypp(k+1) + u(k) enddo ! Store output, if requested if (present(d2ydx2)) d2ydx2 = ypp if ( present(store) ) then if ( .not. store ) return end if ! Store spline info in output variable allocate(dat%x(n), dat%y(n), dat%d2ydx2(n)) dx = (x(n)-x(1))/(n-1) xtol = abs(dx)*tol dat%mesh = meshType dat%xmin = min(x(1),x(n)) - xtol dat%xmax = max(x(1),x(n)) + xtol dat%x1 = x(1) dat%dx = dx dat%a = a if (dat%mesh == 'log') dat%b = b dat%n = n dat%x(:) = x(:) dat%y(:) = y(:) dat%d2ydx2(:) = ypp(:) end subroutine generate_spline_master !------------------------------------------------------------------------------- SUBROUTINE evaluate_spline( dat, x, y, dydx ) ! Evaluate function at one point implicit none type(spline_t), intent(in) :: dat ! info for spline interpolation of y(x) real(dp), intent(in) :: x ! point at which function is needed real(dp), intent(out):: y ! function value at point x real(dp),optional,intent(out):: dydx ! function derivative at point x ! Internal variables integer :: kh, kl real(dp):: xh, xl ! Find mesh interval of point x call find_interval( dat, dat%x, x, kl, kh, xl, xh ) ! Find interpolation within given interval call interpolate_interval( xl, xh, dat%y(kl), dat%y(kh), & dat%d2ydx2(kl), dat%d2ydx2(kh), x, y, dydx ) END SUBROUTINE evaluate_spline !------------------------------------------------------------------------------- SUBROUTINE find_interval( dat, xi, x, kl, kh, xl, xh ) ! Find interpolation interval implicit none type(spline_t), intent(in) :: dat ! spline info real(dp), intent(in) :: xi(:) ! mesh points real(dp), intent(in) :: x ! point at which function is needed integer, intent(out):: kl ! lower interval index integer, intent(out):: kh ! higher interval index real(dp), intent(out):: xl ! lower interval mesh point real(dp), intent(out):: xh ! higher interval mesh point ! Internal variables character(len=*),parameter:: myName = 'evaluate_spline/find_interval ' real(dp):: a, b, dx, h, x1, xmin, xmax, xtol integer :: k, n ! Get some parameters xmin = dat%xmin xmax = dat%xmax n = size(xi) ! Check that point is within interpolation range if (x<dat%xmin .or. x>dat%xmax) & call die(myName//'ERROR: x out of range') ! Find mesh interval of point x select case(dat%mesh) case('lin') ! linear mesh x1 = dat%x1 dx = dat%dx kl = 1+floor((x-x1)/dx) kl = min(n-1,max(1,kl)) kh = kl+1 xl = x1+(kl-1)*dx xh = x1+(kh-1)*dx case ('log') ! logarithmic mesh x1 = dat%x1 a = dat%a b = dat%b kl = 1+floor(log(1+(x-x1)/b)/a) kl = min(n-1,max(1,kl)) kh = kl+1 xl = x1+b*(exp((kl-1)*a)-1) xh = x1+b*(exp((kh-1)*a)-1) case('gen') ! general mesh => use bisection kl = 1 kh = n do while(kh-kl>1) k = (kh+kl)/2 if (xi(k)>x) then kh = k else kl = k endif enddo xl = xi(kl) xh = xi(kh) case default call die(myName//'ERROR: unknown mesh type') end select END SUBROUTINE find_interval !------------------------------------------------------------------------------- SUBROUTINE interpolate_interval( xl, xh, yl, yh, d2ydx2l, d2ydx2h, x, y, dydx ) ! Evaluate function at point x, within a given interval implicit none real(dp), intent(in) :: xl ! lower interval point real(dp), intent(in) :: xh ! higher interval point real(dp), intent(in) :: yl ! function value at xl real(dp), intent(in) :: yh ! function value at xh real(dp), intent(in) :: d2ydx2l ! d2y/dx2 at at xl real(dp), intent(in) :: d2ydx2h ! d2y/dx2 at at xh real(dp), intent(in) :: x ! point at which function is needed real(dp), intent(out):: y ! function value at point x real(dp),optional,intent(out):: dydx ! function derivative at point x ! Internal variables real(dp):: A, B, dAdx, dBdx, dx ! Find interpolated value dx = xh-xl A = (xh-x)/dx B = (x-xl)/dx y = A*yl + B*yh + ( (A**3-A)*d2ydx2l + (B**3-B)*d2ydx2h ) * (dx**2)/6.0_dp ! Find interpolated derivative if (present(dydx)) then dAdx = -1._dp/dx dBdx = +1._dp/dx dydx = dAdx*yl + dBdx*yh + ( dAdx*(3*A**2-1._dp)*d2ydx2l + & dBdx*(3*B**2-1._dp)*d2ydx2h ) * (dx**2)/6.0_dp endif END SUBROUTINE interpolate_interval !------------------------------------------------------------------------------- SUBROUTINE evaluate_spline_n( dat, x, y, dydx ) ! Evaluate a function at several points implicit none type(spline_t), intent(in) :: dat ! info for spline interpolation of y(x) real(dp), intent(in) :: x(:) ! point at which function is needed real(dp), intent(out):: y(:) ! function value at point x real(dp),optional,intent(out):: dydx(:) ! function derivative at point x integer:: ix, nx nx = size(x) do ix = 1,nx call evaluate_spline( dat, x(ix), y(ix), dydx(ix) ) enddo END SUBROUTINE evaluate_spline_n !------------------------------------------------------------------------------- SUBROUTINE generate_spline_x( x, y, n, dydx1, dydxn, d2ydx2 ) ! Included for compatibility with an older interface implicit none integer, intent(in) :: n ! number of mesh points real(dp),intent(in) :: x(n) ! mesh of independent variable real(dp),intent(in) :: y(n) ! function value at mesh points real(dp),intent(in) :: dydx1 ! dy/dx at x(1) real(dp),intent(in) :: dydxn ! dy/dx at x(n) real(dp),intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points ! Internal variables integer :: stat type(spline_t):: dat ! Generate spline dat if (dydx1>dydxMax .and. dydxn>dydxMax) then call generate_spline_master( dat, x, y, n, d2ydx2=d2ydx2, store=.false., stat=stat) elseif (dydxn>dydxMax) then call generate_spline_master( dat, x, y, n, dydx1=dydx1, d2ydx2=d2ydx2, store=.false., stat=stat) elseif (dydx1>dydxMax) then call generate_spline_master( dat, x, y, n, dydxn=dydxn, d2ydx2=d2ydx2, store=.false., stat=stat) else call generate_spline_master( dat, x, y, n, dydx1, dydxn, d2ydx2=d2ydx2, store=.false., stat=stat) endif ! If the status is faulty the resulting d2y/dx2 will be ! forcefully set to 0 ! Then the user has to do something differently if ( stat /= 0 ) d2ydx2(:) = 0._dp ! Deallocate arrays in the spline data-structure call clean_spline(dat) END SUBROUTINE generate_spline_x !------------------------------------------------------------------------------- SUBROUTINE generate_spline_dx( dx, y, n, dydx1, dydxn, d2ydx2 ) ! Included for compatibility with an older interface implicit none integer, intent(in) :: n ! number of mesh points real(dp),intent(in) :: dx ! mesh interval real(dp),intent(in) :: y(n) ! function value at mesh points real(dp),intent(in) :: dydx1 ! dy/dx at x(1) real(dp),intent(in) :: dydxn ! dy/dx at x(n) real(dp),intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points ! Internal variables integer :: ix real(dp):: x(n) ! Generate spline data x = (/( (ix-1)*dx, ix=1,n )/) call generate_spline_x( x, y, n, dydx1, dydxn, d2ydx2 ) END SUBROUTINE generate_spline_dx !------------------------------------------------------------------------------- SUBROUTINE evaluate_spline_x( xi, yi, d2ydx2, n, x, y, dydx ) ! Included for compatibility with an older interface implicit none integer, intent(in) :: n ! number of mesh points real(dp), intent(in) :: xi(n) ! mesh of independent variable real(dp), intent(in) :: yi(n) ! function value at mesh points real(dp), intent(in) :: d2ydx2(n) ! function value at point x real(dp), intent(in) :: x ! point at which function is needed real(dp), intent(out):: y ! function value at point x real(dp),optional,intent(out):: dydx ! function derivative at point x integer :: kh, kl real(dp):: xh, xl, xtol type(spline_t):: dat ! Find mesh interval of point x (warning: no check that x is within range) dat%mesh = 'gen' dat%x1 = xi(1) dat%dx = (xi(n)-xi(1))/(n-1) xtol = abs(dat%dx)*tol dat%xmin = min(xi(1),xi(n)) - xtol dat%xmax = max(xi(1),xi(n)) + xtol call find_interval( dat, xi, x, kl, kh, xl, xh ) ! Find interpolation within given interval call interpolate_interval( xl, xh, yi(kl), yi(kh), & d2ydx2(kl), d2ydx2(kh), x, y, dydx ) END SUBROUTINE evaluate_spline_x !------------------------------------------------------------------------------- SUBROUTINE evaluate_spline_dx( dx, yi, d2ydx2, n, x, y, dydx ) ! Included for compatibility with an older interface implicit none integer, intent(in) :: n ! number of mesh points real(dp), intent(in) :: dx ! mesh interval real(dp), intent(in) :: yi(n) ! function value at mesh points real(dp), intent(in) :: d2ydx2(n) ! function value at point x real(dp), intent(in) :: x ! point at which function is needed real(dp), intent(out):: y ! function value at point x real(dp),optional,intent(out):: dydx ! function derivative at point x integer :: kh, kl real(dp):: xh, xl, xmax, xtol ! Check that point is within interpolation range xmax = (n-1)*dx xtol = dx*tol if (x<-xtol .or. x>xmax+xtol) & call die('evaluate_spline ERROR: x out of range') ! Find mesh interval of point x (warning: no check that x is within range) kl = 1+floor(x/dx) kl = min(n-1,max(1,kl)) kh = kl+1 xl = (kl-1)*dx xh = (kh-1)*dx ! Find interpolation within given interval call interpolate_interval( xl, xh, yi(kl), yi(kh), & d2ydx2(kl), d2ydx2(kh), x, y, dydx ) END SUBROUTINE evaluate_spline_dx !------------------------------------------------------------------------------- SUBROUTINE polint( xi, yi, n, x, y, dydx ) ! Lagrange interpolation implicit none integer, intent(in) :: n ! number of mesh points real(dp),intent(in) :: xi(*) ! mesh points real(dp),intent(in) :: yi(*) ! function values at mesh points real(dp),intent(in) :: x ! point at which function is required real(dp),intent(out):: y ! function value at point x real(dp),intent(out):: dydx ! function derivative at point x integer :: i0, m, nm real(dp):: c(n), d(n), e(n), dcdx(n), dddx(n), dedx(n), dx1(n), dx2(n), dxi(n) ! Neville's algorithm, as described in NR c = yi(1:n) d = yi(1:n) dcdx = 0._dp dddx = 0._dp dydx = 0._dp i0 = n/2 y = yi(i0) do m = 1,n-1 ! loop on increasing polynomial order nm = n-m dxi(1:nm) = xi(1:nm)-xi(m+1:n) if (any(dxi(1:nm)==0.0_dp)) & call die('polint: ERROR: two mesh points are equal') dx1(1:nm) = xi(1:nm)-x dx2(1:nm) = xi(m+1:n)-x e(1:nm) = (c(2:nm+1)-d(1:nm))/dxi(1:nm) c(1:nm) = dx1(1:nm)*e(1:nm) d(1:nm) = dx2(1:nm)*e(1:nm) dedx(1:nm) = (dcdx(2:nm+1) - dddx(1:nm)) / dxi(1:nm) dcdx(1:nm) = - e(1:nm) + dx1(1:nm)*dedx(1:nm) dddx(1:nm) = - e(1:nm) + dx2(1:nm)*dedx(1:nm) i0 = i0-1 if (2*i0 < nm) then i0 = i0+1 y = y+c(i0) dydx = dydx+dcdx(i0) else y = y+d(i0) dydx = dydx+dddx(i0) endif end do END SUBROUTINE polint subroutine clean_spline(dat) type(spline_t), intent(inout) :: dat integer :: n n = dat%n if (allocated(dat%x)) then deallocate(dat%x, dat%y, dat%d2ydx2) endif dat%n = 0 end subroutine clean_spline END MODULE gridxc_interpolation