SUBROUTINE polint(xa, ya, x, y, dy)
USE nrtype; USE nrutil, ONLY : assert_eq, iminloc, nrerror
IMPLICIT NONE
REAL(sp), DIMENSION(:), INTENT(in) :: xa, ya
REAL(sp), INTENT(in) :: x
REAL(sp), INTENT(out) :: y, dy
! Given arrays xa and ya of length N, and given a value x, this routine
! returns a value y, and an error estimate dy. If P(x) is the polynomial
! of degree N - 1 such that P(x a_i) = y a_i, i = 1, ..., N, then the
! returned value y = P(x).
INTEGER(i4b) :: m, n, ns
REAL(sp), DIMENSION(SIZE(xa)) :: c, d, den, ho
n = assert_eq(SIZE(xa), SIZE(ya), 'polint')
c = ya ! Initialize the tableau of c's and d's
d = ya
ho = xa - x
ns = iminloc(ABS(x - xa)) ! Find index ns of closest table entry
y = ya(ns) ! Initial approximation to y.
ns = ns - 1
DO m = 1, n - 1 ! For each column of the tableau
den(1:n-m) = ho(1:n-m) - ho(1+m:n) ! we loop over c's and d's and
IF (ANY(den(1:n-m) == 0.0)) & ! update them
CALL nrerror('polint: calculation failure')
! This error can occur only if two input xa's are (to within roundoff)
! identical.
den(1:n - m) = (c(2:n-m+1) - d(1:n-m))/den(1:n-m)
d(1:n-m) = ho(1+m:n) * den(1:n-m) ! Here c's and d's get updated
c(1:n-m) = ho(1:n-m) * den(1:n-m)
IF (2 * ns < n-m) THEN ! After each column in the tableau is
dy=c(ns+1) ! completed decide, which correction
ELSE ! c or d we add to y. We take the
dy=d(ns) ! straightest line through the tableau
ns=ns-1 ! to its apex. The partial approximations
END IF ! are thus centred on x. The last dy
y = y+dy ! is the measure of error.
END DO
END SUBROUTINE polint
SUBROUTINE ratint(xa, ya, x, y, dy)
USE nrtype; USE nrutil, ONLY : assert_eq, iminloc, nrerror
IMPLICIT NONE
REAL(sp), DIMENSION(:), INTENT(in) :: xa, ya
REAL(sp), INTENT(in) :: x
REAL(sp), INTENT(out) :: y, dy
! Given arrays xa and ya of length N, and given a value of x, this routine
! returns a value of y and an accuracy estimate dy. The value returned is
! that of the diagonal rational function, evaluated at x, that passes
! through the N points (xa_i, ya_i), i = 1... N.
INTEGER(i4b) :: m, n, ns
REAL(sp), DIMENSION(SIZE(xa)) :: c, d, dd, h, t
REAL(sp), PARAMETER :: tiny=1.0e-25_sp
n = assert_eq(SIZE(xa), SIZE(ya), 'ratint')
h = xa - x
ns = iminloc(ABS(h))
y = ya(ns)
IF (x == xa(ns)) THEN
dy = 0.0
RETURN
END IF
c = ya
d = ya + tiny ! The tiny is needed to prevent 0/0
ns = ns - 1
DO m=1, n-1
t(1:n-m) = (xa(1:n-m)-x) * d(1:n-m)/h(1+m:n) ! h will never be 0
dd(1:n-m) = t(1:n-m) - c(2:n-m+1)
IF (ANY(dd(1:n-m) == 0.0)) & ! interpolating function
CALL nrerror('failure in ratint') ! has a pole here
dd(1:n-m) = (c(2:n-m+1) - d(1:n-m))/dd(1:n-m)
d(1:n-m) = c(2:n-m+1) * dd(1:n-m)
c(1:n-m) = t(1:n-m) * dd(1:n-m)
IF(2*ns < n-m) THEN
dy = c(ns+1)
ELSE
dy = d(ns)
ns = ns-1
END IF
y=y+dy
END DO
END SUBROUTINE ratint
SUBROUTINE pzextr(iest, xest, yest, yz, dy)
USE nrtype; USE nrutil, ONLY : assert_eq, nrerror
IMPLICIT NONE
INTEGER(i4b), ikntent(in) :: iest
REAL(sp), INTENT(in) :: xest
REAL(sp), DIMENSION(:), INTENT(in) :: yest
REAL(sp), DIMENSION(:), INTENT(out) :: yz, dy
! Use polynomial extrapolation to evaluate N functions at x = 0 by fitting
! a polynomial to a sequence of estimates with progressively smaller
! values x = xest, and corresponding function vectors yest. This call
! is number iest in the sequence of calls. Extrapolated function values
! are output as yz, and their estimated error is output as dy. yest, yz,
! and dy are arrays of length N.
INTEGER(i4b), PARAMETER :: iest_max = 16
INTEGER(i4b) :: j, nv
INTEGER(i4b), SAVE :: nvold = -1
REAL(sp) :: delta, f1, f2
REAL(sp), DIMENSION(SIZE(yz)) :: d, tmp, q
REAL(sp), DIMENSION(iest_max), SAVE :: x
REAL(sp), DIMENSION(:,:), ALLOCATABLE, SAVE :: qcol
nv = assert_eq(SIZE(yz), SIZE(yest), SIZE(dy), 'pzextr')
IF (iest > iest_max) &
CALL nrerror('pzextr: probable misuse, too much extrapolation')
IF(nv /= nvold) THEN ! set up internal storage
IF(ALLOCATED(qcol)) DEALLOCATE(qcol)
ALLOCATE(qcol(nv, iest_max))
nvold=nv
END IF
x(iest) = xest ! save current independent variable
dy(:) = yest(:)
yz(:) = yest(:)
IF (iest == 1) THEN ! store first estimate in first column
qcol(:,1)=yest(:)
ELSE
d(:) = yest(:)
DO j=1, iest-1
delta=1.0_sp/(x(iest-j)-xest)
f1=xest*delta
f2=x(iest-j)*delta
q(:) = qcol(:,j) ! propagate tableau 1 diagonal more
qcol(:,j)=dy(:)
tmp(:)=d(:)-q(:)
dy(:)=f1*tmp(:)
d(:)=f2*tmp(:)
yz(:)=yz(:)+dy(:)
END DO
qcol(:,iest)=dy(:)
END IF
END SUBROUTINE pzextr
SUBROUTINE rzextr(iest, xest, yest, yz, dy)
USE nrtype; USE nrutil, ONLY : assert_eq, nrerror
IMPLICIT NONE
INTEGER(i4b), INTENT(in) :: iest
REAL(sp), INTENT(in) :: xest
REAL(sp), DIMENSION(:), INTENT(in) :: yest
REAL(sp), DIMENSION(:), INTENT(out) :: yz, dy
! Exact substitute for pzextr, but uses diagonal rational function
! extrapolation instead of polynomial extrapolation.
INTEGER(i4b), PARAMETER :: iest_max = 16
INTEGER(i4b) :: k, nv
INTEGER(i4b), SAVE :: nvold = -1
REAL(sp), DIMENSION(SIZE(yz)) :: yy, v, c, b, b1, ddy
REAL(sp), DIMENSION(:,:), ALLOCATABLE, SAVE :: d
REAL(sp), DIMENSION(iest_max), SAVE :: fx, x
nv = assert_eq(SIZE(yz), SIZE(dy), SIZE(yest), 'rzextr')
IF(iest > iest_max) &
CALL nrerror('rzextr: probable misuse, too much extrapolation')
IF (nv /= nvold) THEN
IF(ALLOCATED(d)) DEALLOCATE(d)
ALLOCATE(d(nv, iest_max))
nvold=nv
END IF
x(iest)=xest ! save current independent variable
IF(iest == 1) THEN
yz=yest
d(:,1)=yest
dy=yest
ELSE
fx(2:iest)=x(iest-1:1:-1)/xest
yy=yest ! evaluate next diagonal in tableau
v=d(1:nv, 1)
c=yy
d(1:nv, 1)=yy
DO k=2, iest
b1=fx(k)*v
b=b1-c
WHERE(b/=0.0)
b=(c-v)/b
ddy=c*b
c=b1*b
ELSEWHERE ! care needed to avoid division by 0
ddy=v
END WHERE
IF (k/=iest) v=d(1:nv, k)
d(1:nv, k) = ddy
yy=yy+ddy
END DO
dy=ddy
yz=yy
END IF
END SUBROUTINE rzextr