c Polynomial interpolation routines from Numerical Recipes c For use with dfrmint.f, frmint.f c Ken Van Riper, kvr@lanl.gov subroutine polynomial_int (xa, ya, n, no, x, y, dy) c Taken from Numerical Recipes integer nmax, n, no parameter (nmax=7) real*8 xa(n), ya(n), c(nmax), d(nmax) real*8 x, y, dy, dif, dift, w, den, ho, hp integer ns, i, m ns = 1 dif = abs(x - xa(1)) do i = 1,no dift = abs(x - xa(i)) if (dift .lt. dif) then ns = i dif = dift endif c(i) = ya(i) d(i) = ya(i) enddo y = ya(ns) ns = ns - 1 do m = 1,no-1 do i = 1,no-m ho = xa(i) - x hp = xa(i+m) - x w = c(i+1) - d(i) den = ho - hp if (den .eq. 0) then write (6,'(a)') 'Trouble in polynomial_int' stop endif den = w/den d(i) = hp*den c(i) = ho*den enddo if (2*ns .lt. no-m) then dy = c(ns+1) else dy = d(ns) ns = ns -1 endif y = y + dy enddo return end subroutine TwoDpolyint (etas, betas, ya, m, n, mo, no, \$ eta, beta, y, dy) c Taken from Numerical Recipes c M for eta direction, N for beta direction integer nmax, mmax, m, n, mo, no parameter (nmax=7, mmax=7) real*8 etas(m), betas(n), ya(m,n), yntmp(nmax), ymtmp(mmax), \$ eta, beta, dy, y integer j, k do j = 1,mo do k = 1,no yntmp(k) = ya(j,k) enddo call polynomial_int (betas, yntmp, n, no, beta, ymtmp(j), dy) enddo call polynomial_int (etas, ymtmp, m, mo, eta, y, dy) return end