1!!@LICENSE 2! 3! ============================================================================== 4! MODULE interpolation 5! ============================================================================== 6! Cubic spline interpolation utilities 7! ============================================================================== 8! Public types in this module: 9! spline_t : derived type to hold spline info of a function 10! Public parameters, variables and arrays in this module: 11! none 12! Public procedures provided by this module: 13! generate_spline : generate spline interpolation info of a function 14! evaluate_spline : evaluate a function at given point(s) using splines 15! ============================================================================== 16! SUBROUTINE generate_spline( dat, x, y, n, dydx1, dydxn, d2ydx2, stat ) 17! Generate data required to interpolate a 1D function with cubic splines 18! Required input: 19! real(dp) x(n) : independent variable at each point 20! real(dp) y(n) : function values at each point 21! integer n : number of points 22! Optional input: 23! real(dp) dydx1 : dy/dx at x(1) (if not present, assumes d2y/dx2=0 at x(1)) 24! real(dp) dydxn : dy/dx at x(n) (if not present, assumes d2y/dx2=0 at x(n)) 25! Output: 26! type(spline_t) dat : spline info of function y(x) 27! Optional output: 28! real(dp) d2ydx2(n) : d2y/dx2 at mesh points 29! integer stat : error status 30! Behaviour: 31! If dydx1 and/or dydxn are not present, assumes d2y/dx2=0 at x(1) and/or x(n) 32! Returns with stat=-1 if mesh is not monotonic, without any other output 33! Returns with stat=-2 if n<2, without any other output 34! Algorithm: 35! Computes d2y/dx2 at each point, to make y(x) and dy/dx continuous. 36! The x values are analyzed to check if the mesh is linear, logarithmic, 37! or general. This is used by evaluate_spline to find the mesh interval at 38! which the evaluation points lie. 39! Ref: "Numerical Recipes", W.H. Press et al, Cambridge U.P. 40! ============================================================================== 41! SUBROUTINE evaluate_spline( dat, x, y ) 42! Evaluate function at given point(s) using spline data info 43! Input: 44! type(spline_t) dat : spline info of function y(x), from generate_spline 45! real(dp) x(:) : point(s) at which function must be evaluated 46! Output: 47! real(dp) y(:) : value of function at given point(s) 48! Behaviour: 49! Single- and multiple-point evaluators are overloaded with the same name, 50! so that x and y may be also scalar values 51! Stops with an error message if any point x is out of the interp. interval 52! ============================================================================== 53! subroutine clean_spline(dat) 54! type(spline_t), intent(inout) :: dat 55! Behavior: 56! Deallocates the array components and resets %n to zero 57! end subroutine clean_spline 58! ============================================================================== 59! SUBROUTINE spline( x, y, n, dydx1, dydxn, d2ydx2 ) 60! Included for compatibility with Numerical Recipes interface 61! Input: 62! real(dp) x(n) ! mesh points 63! real(dp) y(n) ! function value at mesh points 64! integer n ! number of mesh points 65! real(dp) dydx1 ! dy/dx at x(1). If >1.e30, assumes d2y/dx2=0 at x(1) 66! real(dp) dydxn ! dy/dx at x(n). If >1.e30, assumes d2y/dx2=0 at x(n) 67! Output: 68! real(dp) d2ydx2(n) ! d2y/dx2 at mesh points 69! ============================================================================== 70! SUBROUTINE spline( dx, y, n, dydx1, dydxn, d2ydx2 ) 71! Included for compatibility with interface used in siesta 72! Input: 73! real(dp) dx ! mesh interval of a uniform mesh starting at x=0 74! real(dp) y(n) ! function value at mesh points 75! integer n ! number of mesh points 76! real(dp) dydx1 ! dy/dx at x(1). If >1.e30, assumes d2y/dx2=0 at x(1) 77! real(dp) dydxn ! dy/dx at x(n). If >1.e30, assumes d2y/dx2=0 at x(n) 78! Output: 79! real(dp) d2ydx2(n) ! d2y/dx2 at mesh points 80! ============================================================================== 81! SUBROUTINE splint( xi, yi, d2ydx2, n, x, y, dydx ) 82! Included for compatibility with Numerical Recipes interface 83! Input: 84! real(dp) xi(n) ! mesh points 85! real(dp) yi(n) ! function value at mesh points 86! real(dp) d2ydx2(n) ! second derivative at mesh points 87! integer n ! number of mesh points 88! real(dp) x ! point at which function is needed 89! Output: 90! real(dp) y ! function value at point x 91! real(dp) dydx ! function derivative at point x 92! ============================================================================== 93! SUBROUTINE splint( dx, yi, d2ydx2, n, x, y, dydx ) 94! Included for compatibility with interface used in siesta 95! Input: 96! real(dp) dx ! mesh interval of a uniform mesh starting at x=0 97! real(dp) yi(n) ! function value at mesh points 98! real(dp) d2ydx2(n) ! second derivative at mesh points 99! integer n ! number of mesh points 100! real(dp) x ! point at which function is needed 101! Output: 102! real(dp) y ! function value at point x 103! real(dp) dydx ! function derivative at point x 104! ============================================================================== 105! SUBROUTINE polint(XA,YA,N,X,Y,DYDX) 106! Lagrange interpolation 107! Input: 108! real*8 XA(N) : x values of the function y(x) to interpolate 109! real*8 YA(N) : y values of the function y(x) to interpolate 110! integer N : Number of data points 111! real*8 X : x value at which the interpolation is desired 112! Output: 113! real*8 Y : interpolated value of y(x) at X 114! real*8 DYDX : interpolated derivative dy/dx at X 115! Notice: this argument has a different meaning 116! in the Numerical Recipes' polint subroutine 117! Ref: 118! W.H.Press et al, Numerical Recipes, Cambridge U.P. 119! ============================================================================== 120! Written by J.M.Soler and A.Garcia. May.2015 121! ============================================================================== 122! 123! This module calls an external subroutine 'die', with interface 124! 125! interface 126! subroutine die(str) 127! character(len=*), intent(in) :: str 128! end subroutine die 129! end interface 130 131MODULE interpolation 132 133 implicit none 134 135 integer, parameter :: dp = selected_real_kind(10,100) 136 137 138PUBLIC :: & 139 spline_t, &! derived type to hold spline info of a function 140 generate_spline, &! generate spline info 141 evaluate_spline, &! evaluate a function at given point(s) 142 clean_spline, &! deallocate the components of a spline_t object 143 spline, &! compatible with Numerical Recipes interface 144 splint, &! compatible with Numerical Recipes interface 145 polint ! Lagrange polynomial interpolation 146 147PRIVATE ! nothing is declared public beyond this point 148 149! Derived type to hold spline info of a function 150type spline_t 151 private 152 character(len=3):: mesh ! mesh type ('lin'|'log'|'gen') 153 real(dp) :: xmin ! x(1)-xtol 154 real(dp) :: xmax ! x(n)+xtol 155 real(dp) :: x1 ! x(1) 156 real(dp) :: dx ! linear mesh interval 157 real(dp) :: a,b ! log-mesh parameters x(k)=b*[exp(a*(k-1))-1] 158 integer :: n ! number of points 159 real(dp),allocatable:: x(:) ! mesh points 160 real(dp),allocatable:: y(:) ! function value at mesh points 161 real(dp),allocatable:: d2ydx2(:) ! 2nd derivative at mesh points 162end type 163 164! Overloaded spline generators and evaluators 165interface generate_spline 166 module procedure generate_spline_master ! new interface 167 module procedure generate_spline_x ! Numerical Recipes interface 168 module procedure generate_spline_dx ! older interface 169end interface generate_spline 170interface evaluate_spline 171 module procedure evaluate_spline ! evaluate function at single point 172 module procedure evaluate_spline_n ! evaluate function at multiple points 173 module procedure evaluate_spline_x ! Numerical Recipes interface 174 module procedure evaluate_spline_dx ! siesta interface 175end interface evaluate_spline 176interface spline 177 module procedure generate_spline_x ! Numerical Recipes interface 178 module procedure generate_spline_dx ! siesta interface 179end interface spline 180interface splint 181 module procedure evaluate_spline_x ! Numerical Recipes interface 182 module procedure evaluate_spline_dx ! siesta interface 183end interface splint 184 185! Internal module parameters 186real(dp),parameter:: dydxMax = 0.99e30_dp ! max. dy/dx at x(1) and x(n) 187real(dp),parameter:: tol = 1.e-6_dp ! tolerance for 'within range' condition 188 189CONTAINS 190 191!------------------------------------------------------------------------------- 192 193SUBROUTINE generate_spline_master( dat, x, y, n, dydx1, dydxn, d2ydx2, store, stat ) 194 195! Generate data for cubic spline interpolation of a function 196 197implicit none 198type(spline_t), intent(out):: dat ! spline interpolation info 199integer, intent(in) :: n ! number of mesh points 200real(dp), intent(in) :: x(n) ! independent variable at mesh points 201real(dp), intent(in) :: y(n) ! function value at mesh points 202real(dp),optional,intent(in) :: dydx1 ! dy/dx at x(1) 203real(dp),optional,intent(in) :: dydxn ! dy/dx at x(n) 204real(dp),optional,intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points 205logical, optional, intent(in):: store ! Store data in the spline object 206integer, optional,intent(out):: stat ! error status: 207 ! (-1 if nonmonotonic mesh) 208 ! (-2 if n < 2) 209 210! Internal variables and arrays 211character(len=*),parameter:: myName = 'generate_spline_master ' 212real(dp):: a, b, dx, dx1, dx2, dxn, dxm, dxp, dy1, dyn, dym, dyp, & 213 p, s, u(n), v(n), xtol, ypp(n) 214integer :: flag, k 215character(len=3):: meshType 216 217! Check that mesh is monotonic 218if (n<2) then 219 flag = -2 ! single-point 220elseif (n>2 .and. any( (x(2:n-1)-x(1:n-2)) * (x(3:n)-x(2:n-1)) <= 0.0_dp )) then 221 flag = -1 ! non-monotonic mesh 222else 223 flag = 0 224endif 225if (present(stat)) stat = flag 226if (flag<0) return 227 228! Find if mesh is linear (a=0) or logarithmic 229! x(k)=x(1)+b*[exp(a*(k-1))-1] => (x(k+1)-x(k))/(x(k)-x(k-1))=exp(a) 230if (n<3) then 231 meshType = 'lin' 232 dx = x(2)-x(1) 233else 234 a = log( (x(3)-x(2)) / (x(2)-x(1)) ) 235 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 236 if (abs(a)<1.e-12_dp) then 237 meshType = 'lin' 238 dx = x(2)-x(1) 239 else ! try x(k) = x(1) + b*[exp(a*(k-1))-1] 240 meshType = 'log' 241 b = (x(2)-x(1)) / (exp(a)-1._dp) 242 endif 243 else 244 meshType = 'gen' 245 endif 246endif 247 248! Set boundary conditions at end points 249if (present(dydx1)) then 250 dx1 = x(2)-x(1) 251 dy1 = y(2)-y(1) 252 v(1) = -0.5_dp 253 u(1) = (3.0_dp/dx1)*(dy1/dx1-dydx1) 254else 255 v(1) = 0._dp 256 u(1) = 0._dp 257endif 258if (present(dydxn)) then 259 dxn = x(n)-x(n-1) 260 dyn = y(n)-y(n-1) 261 v(n) = 0.5_dp 262 u(n) = (3.0_dp/dxn)*(dydxn-dyn/dxn) 263else 264 v(n) = 0.0_dp 265 u(n) = 0.0_dp 266endif 267 268! Decomposition loop of the tridiagonal equations 269do k = 2,n-1 270 dxm = x(k)-x(k-1) 271 dym = y(k)-y(k-1) 272 dxp = x(k+1)-x(k) 273 dyp = y(k+1)-y(k) 274 s = dxm / (dxm+dxp) 275 p = s*v(k-1) + 2.0_dp 276 v(k) = (s-1.0_dp)/p 277 u(k) = ( 6.0_dp*(dyp/dxp-dym/dxm)/(dxm+dxp) - s*u(k-1) )/p 278enddo 279 280! Backsubstitution loop of the tridiagonal equations 281ypp(n) = (u(n)-v(n)*u(n-1)) / (v(n)*v(n-1)+1.0_dp) 282do k = n-1,1,-1 283 ypp(k) = v(k)*ypp(k+1) + u(k) 284enddo 285 286! Store output, if requested 287if (present(d2ydx2)) d2ydx2 = ypp 288 289if ( present(store) ) then 290 if ( .not. store ) return 291end if 292 293! Store spline info in output variable 294allocate(dat%x(n), dat%y(n), dat%d2ydx2(n)) 295 296dx = (x(n)-x(1))/(n-1) 297xtol = abs(dx)*tol 298dat%mesh = meshType 299dat%xmin = min(x(1),x(n)) - xtol 300dat%xmax = max(x(1),x(n)) + xtol 301dat%x1 = x(1) 302dat%dx = dx 303dat%a = a 304if (dat%mesh == 'log') dat%b = b 305dat%n = n 306dat%x(:) = x(:) 307dat%y(:) = y(:) 308dat%d2ydx2(:) = ypp(:) 309 310end subroutine generate_spline_master 311 312!------------------------------------------------------------------------------- 313 314SUBROUTINE evaluate_spline( dat, x, y, dydx ) ! Evaluate function at one point 315 316implicit none 317type(spline_t), intent(in) :: dat ! info for spline interpolation of y(x) 318real(dp), intent(in) :: x ! point at which function is needed 319real(dp), intent(out):: y ! function value at point x 320real(dp),optional,intent(out):: dydx ! function derivative at point x 321 322! Internal variables 323integer :: kh, kl 324real(dp):: xh, xl 325 326! Find mesh interval of point x 327call find_interval( dat, dat%x, x, kl, kh, xl, xh ) 328 329! Find interpolation within given interval 330call interpolate_interval( xl, xh, dat%y(kl), dat%y(kh), & 331 dat%d2ydx2(kl), dat%d2ydx2(kh), x, y, dydx ) 332 333END SUBROUTINE evaluate_spline 334 335!------------------------------------------------------------------------------- 336 337SUBROUTINE find_interval( dat, xi, x, kl, kh, xl, xh ) 338 339! Find interpolation interval 340 341implicit none 342type(spline_t), intent(in) :: dat ! spline info 343real(dp), intent(in) :: xi(:) ! mesh points 344real(dp), intent(in) :: x ! point at which function is needed 345integer, intent(out):: kl ! lower interval index 346integer, intent(out):: kh ! higher interval index 347real(dp), intent(out):: xl ! lower interval mesh point 348real(dp), intent(out):: xh ! higher interval mesh point 349 350! Internal variables 351character(len=*),parameter:: myName = 'evaluate_spline/find_interval ' 352real(dp):: a, b, dx, h, x1, xmin, xmax, xtol 353integer :: k, n 354 355! Get some parameters 356xmin = dat%xmin 357xmax = dat%xmax 358n = size(xi) 359 360! Check that point is within interpolation range 361if (x<dat%xmin .or. x>dat%xmax) & 362 call die(myName//'ERROR: x out of range') 363 364! Find mesh interval of point x 365select case(dat%mesh) 366case('lin') ! linear mesh 367 x1 = dat%x1 368 dx = dat%dx 369 kl = 1+floor((x-x1)/dx) 370 kl = min(n-1,max(1,kl)) 371 kh = kl+1 372 xl = x1+(kl-1)*dx 373 xh = x1+(kh-1)*dx 374case ('log') ! logarithmic mesh 375 x1 = dat%x1 376 a = dat%a 377 b = dat%b 378 kl = 1+floor(log(1+(x-x1)/b)/a) 379 kl = min(n-1,max(1,kl)) 380 kh = kl+1 381 xl = x1+b*(exp((kl-1)*a)-1) 382 xh = x1+b*(exp((kh-1)*a)-1) 383case('gen') ! general mesh => use bisection 384 kl = 1 385 kh = n 386 do while(kh-kl>1) 387 k = (kh+kl)/2 388 if (xi(k)>x) then 389 kh = k 390 else 391 kl = k 392 endif 393 enddo 394 xl = xi(kl) 395 xh = xi(kh) 396case default 397 call die(myName//'ERROR: unknown mesh type') 398end select 399 400END SUBROUTINE find_interval 401 402!------------------------------------------------------------------------------- 403 404SUBROUTINE interpolate_interval( xl, xh, yl, yh, d2ydx2l, d2ydx2h, x, y, dydx ) 405 406! Evaluate function at point x, within a given interval 407 408implicit none 409real(dp), intent(in) :: xl ! lower interval point 410real(dp), intent(in) :: xh ! higher interval point 411real(dp), intent(in) :: yl ! function value at xl 412real(dp), intent(in) :: yh ! function value at xh 413real(dp), intent(in) :: d2ydx2l ! d2y/dx2 at at xl 414real(dp), intent(in) :: d2ydx2h ! d2y/dx2 at at xh 415real(dp), intent(in) :: x ! point at which function is needed 416real(dp), intent(out):: y ! function value at point x 417real(dp),optional,intent(out):: dydx ! function derivative at point x 418 419! Internal variables 420real(dp):: A, B, dAdx, dBdx, dx 421 422! Find interpolated value 423dx = xh-xl 424A = (xh-x)/dx 425B = (x-xl)/dx 426y = A*yl + B*yh + ( (A**3-A)*d2ydx2l + (B**3-B)*d2ydx2h ) * (dx**2)/6.0_dp 427 428! Find interpolated derivative 429if (present(dydx)) then 430 dAdx = -1._dp/dx 431 dBdx = +1._dp/dx 432 dydx = dAdx*yl + dBdx*yh + ( dAdx*(3*A**2-1._dp)*d2ydx2l + & 433 dBdx*(3*B**2-1._dp)*d2ydx2h ) * (dx**2)/6.0_dp 434endif 435 436END SUBROUTINE interpolate_interval 437 438!------------------------------------------------------------------------------- 439 440SUBROUTINE evaluate_spline_n( dat, x, y, dydx ) 441 442! Evaluate a function at several points 443 444implicit none 445type(spline_t), intent(in) :: dat ! info for spline interpolation of y(x) 446real(dp), intent(in) :: x(:) ! point at which function is needed 447real(dp), intent(out):: y(:) ! function value at point x 448real(dp),optional,intent(out):: dydx(:) ! function derivative at point x 449 450integer:: ix, nx 451 452nx = size(x) 453do ix = 1,nx 454 call evaluate_spline( dat, x(ix), y(ix), dydx(ix) ) 455enddo 456 457END SUBROUTINE evaluate_spline_n 458 459!------------------------------------------------------------------------------- 460 461SUBROUTINE generate_spline_x( x, y, n, dydx1, dydxn, d2ydx2 ) 462 463! Included for compatibility with an older interface 464 465implicit none 466integer, intent(in) :: n ! number of mesh points 467real(dp),intent(in) :: x(n) ! mesh of independent variable 468real(dp),intent(in) :: y(n) ! function value at mesh points 469real(dp),intent(in) :: dydx1 ! dy/dx at x(1) 470real(dp),intent(in) :: dydxn ! dy/dx at x(n) 471real(dp),intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points 472 473! Internal variables 474integer :: stat 475type(spline_t):: dat 476 477! Generate spline dat 478if (dydx1>dydxMax .and. dydxn>dydxMax) then 479 call generate_spline_master( dat, x, y, n, d2ydx2=d2ydx2, store=.false., stat=stat) 480elseif (dydxn>dydxMax) then 481 call generate_spline_master( dat, x, y, n, dydx1=dydx1, d2ydx2=d2ydx2, store=.false., stat=stat) 482elseif (dydx1>dydxMax) then 483 call generate_spline_master( dat, x, y, n, dydxn=dydxn, d2ydx2=d2ydx2, store=.false., stat=stat) 484else 485 call generate_spline_master( dat, x, y, n, dydx1, dydxn, d2ydx2=d2ydx2, store=.false., stat=stat) 486endif 487 488! If the status is faulty the resulting d2y/dx2 will be 489! forcefully set to 0 490! Then the user has to do something differently 491if ( stat /= 0 ) d2ydx2(:) = 0._dp 492 493! Deallocate arrays in the spline data-structure 494call clean_spline(dat) 495 496END SUBROUTINE generate_spline_x 497 498!------------------------------------------------------------------------------- 499 500SUBROUTINE generate_spline_dx( dx, y, n, dydx1, dydxn, d2ydx2 ) 501 502! Included for compatibility with an older interface 503 504implicit none 505integer, intent(in) :: n ! number of mesh points 506real(dp),intent(in) :: dx ! mesh interval 507real(dp),intent(in) :: y(n) ! function value at mesh points 508real(dp),intent(in) :: dydx1 ! dy/dx at x(1) 509real(dp),intent(in) :: dydxn ! dy/dx at x(n) 510real(dp),intent(out):: d2ydx2(n) ! d2y/dx2 at mesh points 511 512! Internal variables 513integer :: ix 514real(dp):: x(n) 515 516! Generate spline data 517x = (/( (ix-1)*dx, ix=1,n )/) 518call generate_spline_x( x, y, n, dydx1, dydxn, d2ydx2 ) 519 520END SUBROUTINE generate_spline_dx 521 522!------------------------------------------------------------------------------- 523 524SUBROUTINE evaluate_spline_x( xi, yi, d2ydx2, n, x, y, dydx ) 525 526! Included for compatibility with an older interface 527 528implicit none 529integer, intent(in) :: n ! number of mesh points 530real(dp), intent(in) :: xi(n) ! mesh of independent variable 531real(dp), intent(in) :: yi(n) ! function value at mesh points 532real(dp), intent(in) :: d2ydx2(n) ! function value at point x 533real(dp), intent(in) :: x ! point at which function is needed 534real(dp), intent(out):: y ! function value at point x 535real(dp),optional,intent(out):: dydx ! function derivative at point x 536 537integer :: kh, kl 538real(dp):: xh, xl, xtol 539type(spline_t):: dat 540 541! Find mesh interval of point x (warning: no check that x is within range) 542dat%mesh = 'gen' 543dat%x1 = xi(1) 544dat%dx = (xi(n)-xi(1))/(n-1) 545xtol = abs(dat%dx)*tol 546dat%xmin = min(xi(1),xi(n)) - xtol 547dat%xmax = max(xi(1),xi(n)) + xtol 548call find_interval( dat, xi, x, kl, kh, xl, xh ) 549 550! Find interpolation within given interval 551call interpolate_interval( xl, xh, yi(kl), yi(kh), & 552 d2ydx2(kl), d2ydx2(kh), x, y, dydx ) 553 554END SUBROUTINE evaluate_spline_x 555 556!------------------------------------------------------------------------------- 557 558SUBROUTINE evaluate_spline_dx( dx, yi, d2ydx2, n, x, y, dydx ) 559 560! Included for compatibility with an older interface 561 562implicit none 563integer, intent(in) :: n ! number of mesh points 564real(dp), intent(in) :: dx ! mesh interval 565real(dp), intent(in) :: yi(n) ! function value at mesh points 566real(dp), intent(in) :: d2ydx2(n) ! function value at point x 567real(dp), intent(in) :: x ! point at which function is needed 568real(dp), intent(out):: y ! function value at point x 569real(dp),optional,intent(out):: dydx ! function derivative at point x 570 571integer :: kh, kl 572real(dp):: xh, xl, xmax, xtol 573 574! Check that point is within interpolation range 575xmax = (n-1)*dx 576xtol = dx*tol 577if (x<-xtol .or. x>xmax+xtol) & 578 call die('evaluate_spline ERROR: x out of range') 579 580! Find mesh interval of point x (warning: no check that x is within range) 581kl = 1+floor(x/dx) 582kl = min(n-1,max(1,kl)) 583kh = kl+1 584xl = (kl-1)*dx 585xh = (kh-1)*dx 586 587! Find interpolation within given interval 588call interpolate_interval( xl, xh, yi(kl), yi(kh), & 589 d2ydx2(kl), d2ydx2(kh), x, y, dydx ) 590 591END SUBROUTINE evaluate_spline_dx 592 593!------------------------------------------------------------------------------- 594 595SUBROUTINE polint( xi, yi, n, x, y, dydx ) ! Lagrange interpolation 596 597 implicit none 598 integer, intent(in) :: n ! number of mesh points 599 real(dp),intent(in) :: xi(*) ! mesh points 600 real(dp),intent(in) :: yi(*) ! function values at mesh points 601 real(dp),intent(in) :: x ! point at which function is required 602 real(dp),intent(out):: y ! function value at point x 603 real(dp),intent(out):: dydx ! function derivative at point x 604 605 integer :: i0, m, nm 606 real(dp):: c(n), d(n), e(n), dcdx(n), dddx(n), dedx(n), dx1(n), dx2(n), dxi(n) 607 608! Neville's algorithm, as described in NR 609 c = yi(1:n) 610 d = yi(1:n) 611 dcdx = 0._dp 612 dddx = 0._dp 613 dydx = 0._dp 614 i0 = n/2 615 y = yi(i0) 616 do m = 1,n-1 ! loop on increasing polynomial order 617 nm = n-m 618 dxi(1:nm) = xi(1:nm)-xi(m+1:n) 619 if (any(dxi(1:nm)==0.0_dp)) & 620 call die('polint: ERROR: two mesh points are equal') 621 dx1(1:nm) = xi(1:nm)-x 622 dx2(1:nm) = xi(m+1:n)-x 623 e(1:nm) = (c(2:nm+1)-d(1:nm))/dxi(1:nm) 624 c(1:nm) = dx1(1:nm)*e(1:nm) 625 d(1:nm) = dx2(1:nm)*e(1:nm) 626 dedx(1:nm) = (dcdx(2:nm+1) - dddx(1:nm)) / dxi(1:nm) 627 dcdx(1:nm) = - e(1:nm) + dx1(1:nm)*dedx(1:nm) 628 dddx(1:nm) = - e(1:nm) + dx2(1:nm)*dedx(1:nm) 629 i0 = i0-1 630 if (2*i0 < nm) then 631 i0 = i0+1 632 y = y+c(i0) 633 dydx = dydx+dcdx(i0) 634 else 635 y = y+d(i0) 636 dydx = dydx+dddx(i0) 637 endif 638 end do 639 640END SUBROUTINE polint 641 642subroutine clean_spline(dat) 643type(spline_t), intent(inout) :: dat 644 645integer :: n 646 647n = dat%n 648if (allocated(dat%x)) then 649 deallocate(dat%x, dat%y, dat%d2ydx2) 650endif 651dat%n = 0 652end subroutine clean_spline 653 654END MODULE interpolation 655 656