From a31fdc22e4b547e93cf36116fc75c82292d081a0 Mon Sep 17 00:00:00 2001 From: Ole Streicher Date: Fri, 10 Nov 2017 15:41:32 +0100 Subject: [PATCH] Add replacement for spline and polynomial interpolations Spline and polynomial interpolations is used in the `trebin` task. The replacement code (in SPP) is based on the standard methods (described f.e. in Wikipedia) --- pkg/utilities/nttools/doc/trebin.hlp | 4 - pkg/utilities/nttools/trebin/mkpkg | 8 +- pkg/utilities/nttools/trebin/tucspl.f | 52 ------------- pkg/utilities/nttools/trebin/tucspl.x | 99 +++++++++++++++++++++++++ pkg/utilities/nttools/trebin/tuhunt.f | 103 -------------------------- pkg/utilities/nttools/trebin/tuhunt.x | 66 +++++++++++++++++ pkg/utilities/nttools/trebin/tuiep3.f | 71 ------------------ pkg/utilities/nttools/trebin/tuiepn.x | 35 +++++++++ pkg/utilities/nttools/trebin/tuifit.x | 7 +- pkg/utilities/nttools/trebin/tuispl.f | 32 -------- pkg/utilities/nttools/trebin/tuispl.x | 24 ++++++ pkg/utilities/nttools/trebin/tuival.x | 5 +- 12 files changed, 231 insertions(+), 275 deletions(-) delete mode 100644 pkg/utilities/nttools/trebin/tucspl.f create mode 100644 pkg/utilities/nttools/trebin/tucspl.x delete mode 100644 pkg/utilities/nttools/trebin/tuhunt.f create mode 100644 pkg/utilities/nttools/trebin/tuhunt.x delete mode 100644 pkg/utilities/nttools/trebin/tuiep3.f create mode 100644 pkg/utilities/nttools/trebin/tuiepn.x delete mode 100644 pkg/utilities/nttools/trebin/tuispl.f create mode 100644 pkg/utilities/nttools/trebin/tuispl.x diff --git a/pkg/utilities/nttools/doc/trebin.hlp b/pkg/utilities/nttools/doc/trebin.hlp index 293c08ec1..36eb3d181 100644 --- a/pkg/utilities/nttools/doc/trebin.hlp +++ b/pkg/utilities/nttools/doc/trebin.hlp @@ -246,10 +246,6 @@ but it will not modify any scalar column that gives the effective array size. .ih REFERENCES This task was written by Phil Hodge. -The following Numerical Recipes subroutines are used by this task: -TUCSPL, TUHUNT, TUIEP3, and TUISPL. -Numerical Recipes was written by W.H. Press, B.P. Flannery, -S.A. Teukolsky, and W.T. Vetterling. .ih SEE ALSO Type "help tables opt=sys" for a higher-level description of the 'tables' diff --git a/pkg/utilities/nttools/trebin/mkpkg b/pkg/utilities/nttools/trebin/mkpkg index 7b1e83c31..4e85ed549 100644 --- a/pkg/utilities/nttools/trebin/mkpkg +++ b/pkg/utilities/nttools/trebin/mkpkg @@ -11,16 +11,16 @@ libpkg.a: tnamgio.x tnaminit.x trebin.x - tucspl.f + tucspl.x tudcol.x tugcol.x tugetput.x - tuhunt.f - tuiep3.f + tuhunt.x + tuiepn.x tuifit.x trebin.h tuinterp.x tuiset.x trebin.h - tuispl.f + tuispl.x tuival.x trebin.h tutrim.x tuxget.x diff --git a/pkg/utilities/nttools/trebin/tucspl.f b/pkg/utilities/nttools/trebin/tucspl.f deleted file mode 100644 index 511211b89..000000000 --- a/pkg/utilities/nttools/trebin/tucspl.f +++ /dev/null @@ -1,52 +0,0 @@ - subroutine tucspl (xa, ya, n, work, y2) -C -C Compute the second derivative of YA at each point in the array. -C This is the initialization needed in preparation for interpolating -C using cubic splines by the subroutine TUISPL. Input and output -C are all double precision. -C -C This routine was copied with slight modifications from the SPLINE -C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and -C Vetterling. -C -C N i: number of elements in each array -C XA i: array of independent-variable values -C YA i: array of dependent-variable values -C WORK io: scratch array used for work space -C Y2 o: second derivative of YA at each point -C -CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes SPLINE. -C - integer n - double precision xa(n), ya(n), work(n), y2(n) -C-- - integer i - double precision p, sig - -C These values (and y2(n) = 0) are for a "natural" spline. - y2(1) = 0. - work(1) = 0. -C -C This is the decomposition loop of the tridiagonal algorithm. -C Y2 and WORK are used for temporary storage of the decomposed factors. -C - do 10 i = 2, n-1 - sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1)) - p = sig * y2(i-1) + 2. - y2(i) = (sig - 1.) / p - work(i) = (6. * ((ya(i+1) - ya(i)) / (xa(i+1) - xa(i)) - + - (ya(i) - ya(i-1)) / (xa(i) - xa(i-1))) / - + (xa(i+1) - xa(i-1)) - sig * work(i-1)) / p - 10 continue - -C "natural" spline - y2(n) = 0. -C -C This is the backsubstitution loop of the tridiagonal algorithm. -C - do 20 i = n-1, 1, -1 - y2(i) = y2(i) * y2(i+1) + work(i) - 20 continue - - return - end diff --git a/pkg/utilities/nttools/trebin/tucspl.x b/pkg/utilities/nttools/trebin/tucspl.x new file mode 100644 index 000000000..b60efac71 --- /dev/null +++ b/pkg/utilities/nttools/trebin/tucspl.x @@ -0,0 +1,99 @@ +# Setup spline interpolation +# +# Copyright (c) 2017 Ole Streicher + +# Solve a tridiagonal matrix (Thomas algorithm). +# +# Reference: Conte, S.D., and deBoor, C., Elementary Numerical Analysis, +# McGraw-Hill, New York, 1972. +# +procedure solve_tridiagonal(n, a, b, c, d) +int n +double a[n] # subdiagonal +double b[n] # diagonal +double c[n] # superdiagonal +double d[n] # in: right side with constants. out: solution + +int i +pointer cs, ds, sp +begin + call smark(sp) + call salloc(cs, n, TY_DOUBLE) + call salloc(ds, n, TY_DOUBLE) + + # Forward step + Memd[cs] = c[1] / b[1] + do i = 2, n-1 { + Memd[cs+i-1] = c[i] / (b[i] - a[i] * Memd[cs+i-2]) + } + Memd[ds] = d[1] / b[1] + do i = 2, n { + Memd[ds+i-1] = (d[i] - a[i] * Memd[ds+i-2]) / (b[i] - a[i] * Memd[cs+i-2]) + } + + # Back substitution + d[n] = Memd[ds+n-1] + do i = n-1, 1, -1 { + d[i] = Memd[ds+i-1] - Memd[cs+i-1] * d[i+1] + } + + call sfree(sp) +end + + +# Compute the derivative K at each point in the array. This is the +# initialization needed in preparation for interpolating using cubic +# splines by the subroutine TUISPL. Input and output are all double +# precision. +# +# Reference: De Boor, Carl, et al. A practical guide to splines. +# Vol. 27. New York: Springer-Verlag, 1978. +# +# The interface is adopted from IRAFs original tucspl() routine. +# The "work" scratch array is gone, however. +# +procedure tucspl(xa, ya, n, k) + +int n # number of elements in each array +double xa[n] # array of independent-variable values +double ya[n] # array of dependent-variable values +double k[n] # derivative of YA at each point + +int i +double xl, xu +pointer sp, a, b, c + +begin + call smark(sp) + call salloc(a, n, TY_DOUBLE) + call salloc(b, n, TY_DOUBLE) + call salloc(c, n, TY_DOUBLE) + + xu = xa[2] - xa[1] + Memd[a] = 0. + Memd[b] = 2./ xu + Memd[c] = 1. / xu + k[1] = 3 * (ya[2] - ya[1]) / xu**2 + + do i=2, n-1 { + xl = xa[i] - xa[i-1] + xu = xa[i+1] - xa[i] + Memd[a+i-1] = 1. / xl + Memd[b+i-1] = 2. / xl + 2. / xu + Memd[c+i-1] = 1. / xu + k[i] = 3 * ( (ya[i] - ya[i-1])/xl**2 + (ya[i+1] - ya[i])/xu**2 ) + } + + xl = xa[n] - xa[n-1] + Memd[a+n-1] = 1. / xl + Memd[b+n-1] = 2. / xl + Memd[c+n-1] = 0. + k[n] = 3 * (ya[n] - ya[n-1]) / xl**2 + + # Solve triangular matrix + # This replaces the right side with the second derivatives. + call solve_tridiagonal(n, Memd[a], Memd[b], Memd[c], k) + + call sfree(sp) +end + diff --git a/pkg/utilities/nttools/trebin/tuhunt.f b/pkg/utilities/nttools/trebin/tuhunt.f deleted file mode 100644 index 3d8df6480..000000000 --- a/pkg/utilities/nttools/trebin/tuhunt.f +++ /dev/null @@ -1,103 +0,0 @@ - subroutine tuhunt (xa, n, x, klo) -C -C The array XA is searched for an element KLO such that -C xa(klo) <= x <= xa(klo+1) -C If X < XA(1) then KLO is set to zero; if X > XA(N) then KLO is set to N. -C That is, KLO = 0 or N is a flag indicating X is out of bounds. -C -C KLO must be set to an initial guess on input; it will then be replaced -C by the correct value on output. -C -C This routine was copied with some modifications from the HUNT -C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and -C Vetterling. -C -C N i: number of elements in each array -C XA i: array of independent-variable values -C X i: the value to be bracketed by elements in XA -C KLO io: the lower index in XA that brackets X -C -CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes HUNT. -CH Phil Hodge, 21-May-1996 Don't flag endpoints of XA as out of bounds. -C - integer n, klo - double precision xa(n), x -C-- - integer inc, km, khi - logical ascnd - -C Set ASCND, and check for X out of bounds. - - if (xa(n).gt.xa(1)) then - ascnd = .true. - if (x .lt. xa(1)) then - klo = 0 - return - else if (x .gt. xa(n)) then - klo = n - return - end if - else - ascnd = .false. - if (x .gt. xa(1)) then - klo = 0 - return - else if (x .lt. xa(n)) then - klo = n - return - end if - end if - - if ((klo .le. 0) .or. (klo .gt. n)) then - klo = 1 - khi = n - go to 3 - endif - - inc = 1 - if ((x.ge.xa(klo) .and. ascnd) .or. - + (x.lt.xa(klo) .and. .not.ascnd)) then -1 khi = klo + inc - if (khi .gt. n) then - khi = n + 1 - else if ((x.ge.xa(khi) .and. ascnd) .or. - + (x.lt.xa(khi) .and. .not.ascnd)) then - klo = khi - inc = inc + inc - go to 1 - endif - else - khi = klo -2 klo = khi - inc - if (klo .lt. 1) then - klo = 0 - else if ((x.lt.xa(klo) .and. ascnd) .or. - + (x.ge.xa(klo) .and. .not.ascnd)) then - khi = klo - inc = inc + inc - go to 2 - endif - endif - -3 continue -C Before we return, make sure we don't return a value of KLO that -C implies X is out of bounds. We know it isn't because we checked -C at the beginning. - if (khi-klo .eq. 1) then - klo = max (klo, 1) - klo = min (klo, n-1) - return - end if - - km = (khi + klo) / 2 - - if ((x .gt. xa(km) .and. ascnd) .or. - + (x .le. xa(km) .and. .not.ascnd)) then - klo = km - else - khi = km - endif - - go to 3 - - end diff --git a/pkg/utilities/nttools/trebin/tuhunt.x b/pkg/utilities/nttools/trebin/tuhunt.x new file mode 100644 index 000000000..110aa7a54 --- /dev/null +++ b/pkg/utilities/nttools/trebin/tuhunt.x @@ -0,0 +1,66 @@ +# The array XA is searched for an element KLO such that +# xa(klo) <= x <= xa(klo+1) +# If X < XA(1) then KLO is set to zero; if X > XA(N) then KLO is set to N. +# That is, KLO = 0 or N is a flag indicating X is out of bounds. +# +# KLO must be set to an initial guess on input; it will then be replaced +# by the correct value on output. +# +# Copyright (c) 2017 Ole Streicher +# +procedure tuhunt (xa, n, x, klo) +int n # i: number of elements in each array +int klo # i: array of independent-variable values +double xa[n] # i: the value to be bracketed by elements in XA +double x # io: the lower index in XA that brackets X + +int low, up +begin + if (xa[n] > xa[1]) { # increasing + if (xa[1] > x) { + klo = 0 + return + } + if (xa[n] < x) { + klo = n + return + } + klo = min(max(klo, 1), n-1) + low = 1 + up = n-1 + while (true) { + if (x < xa[klo]) { + up = klo-1 + klo = (up + low)/2 + } else if (x > xa[klo+1]) { + low = klo+1 + klo = (up + low)/2 + } else { + return + } + } + } else { # decreasing + if (xa[1] < x) { + klo = 0 + return + } + if (xa[n] > x) { + klo = n + return + } + klo = min(max(klo, 1), n-1) + low = 1 + up = n-1 + while (true) { + if (x > xa[klo]) { + up = klo-1 + klo = (up + low)/2 + } else if (x < xa[klo+1]) { + low = klo+1 + klo = (up + low)/2 + } else { + return + } + } + } +end diff --git a/pkg/utilities/nttools/trebin/tuiep3.f b/pkg/utilities/nttools/trebin/tuiep3.f deleted file mode 100644 index 58f81a99e..000000000 --- a/pkg/utilities/nttools/trebin/tuiep3.f +++ /dev/null @@ -1,71 +0,0 @@ - subroutine tuiep3 (xa, ya, x, y, dy) -C -C Evaluate a cubic polynomial interpolation function at X. -C xa(klo) <= x <= xa(klo+1) -C This routine was copied with slight modifications from the POLINT -C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and -C Vetterling. -C -C XA i: array of four independent-variable values -C YA i: array of four dependent-variable values -C X i: the independent variable -C Y o: the result of interpolations -C DY o: an estimate of the error of interpolation -C -CH Phil Hodge, 18-Apr-1988 Subroutine copied from Numerical Recipes POLINT. -C - double precision xa(4), ya(4), x, y, dy -C-- - integer n -C four-point interpolation - parameter (n = 4) - - double precision c(n), d(n), dif, dift, den - double precision ho, hp, w - integer i, m, ns - - do 10 i = 1, n - if (x .eq. xa(i)) then - y = ya(i) - return - endif - 10 continue - - ns = 1 - dif = abs (x - xa(1)) - - do 20 i = 1, n - dift = abs (x - xa(i)) - if (dift .lt. dif) then - ns = i - dif = dift - endif - c(i) = ya(i) - d(i) = ya(i) - 20 continue - - y = ya(ns) - ns = ns - 1 - - do 40 m = 1, n-1 - - do 30 i = 1, n-m - ho = xa(i) - x - hp = xa(i+m) - x - w = c(i+1) - d(i) - den = w / (ho - hp) - d(i) = hp * den - c(i) = ho * den - 30 continue - - if (2*ns .lt. n-m) then - dy = c(ns+1) - else - dy = d(ns) - ns = ns - 1 - endif - y = y + dy - 40 continue - - return - end diff --git a/pkg/utilities/nttools/trebin/tuiepn.x b/pkg/utilities/nttools/trebin/tuiepn.x new file mode 100644 index 000000000..f53493e46 --- /dev/null +++ b/pkg/utilities/nttools/trebin/tuiepn.x @@ -0,0 +1,35 @@ +# Evaluate a cubic polynomial interpolation function at X. +# xa(klo) <= x <= xa(klo+1) +# +# Copyright (c) 2017 Ole Streicher +# +# References: 1. A. C. Aitken, On interpolation by iteration of +# proportional parts, without the use of differences. +# Proc. Edinburgh Math. Soc, v. 3, 1932, pp. 56-76. +# 2. E. H. Neville, Iterative interpolation. +# Indian Math. Soc, v. 20, 1934, pp. 87-120 +# +procedure tuiepn (n, xa, ya, x, y) + +int n # i: polynomial degree +double xa[n] # i: array of n independent-variable values +double ya[n] # i: array of n dependent-variable values +double x # i: the independent variable +double y # o: the result of interpolation + +int i, j +pointer sp, p + +begin + call smark(sp) + call salloc(p, n, TY_DOUBLE) + + call achtdd(ya, Memd[p], n) + do i=1, n-1 { + do j=1, n-i { + Memd[p+j-1] = ((x - xa[i+j]) * Memd[p+j-1] - (x - xa[j]) * Memd[p+j]) / (xa[j] - xa[i+j]) + } + } + y = Memd[p] + call sfree(sp) +end diff --git a/pkg/utilities/nttools/trebin/tuifit.x b/pkg/utilities/nttools/trebin/tuifit.x index 2f3a2d845..21ba65f3c 100644 --- a/pkg/utilities/nttools/trebin/tuifit.x +++ b/pkg/utilities/nttools/trebin/tuifit.x @@ -22,8 +22,6 @@ double ya[ARB] # o: array of dependent-variable values double y2[ARB] # o: used only by spline interpolation int n # o: size of xa, ya, y2 arrays #-- -pointer sp -pointer work # scratch for spline fit int k # loop index begin @@ -52,10 +50,7 @@ begin case I_SPLINE: if (n >= 4) { - call smark (sp) - call salloc (work, n, TY_DOUBLE) - call tucspl (xa, ya, n, Memd[work], y2) - call sfree (sp) + call tucspl (xa, ya, n, y2) } else { n = 0 } diff --git a/pkg/utilities/nttools/trebin/tuispl.f b/pkg/utilities/nttools/trebin/tuispl.f deleted file mode 100644 index 2f4c68b85..000000000 --- a/pkg/utilities/nttools/trebin/tuispl.f +++ /dev/null @@ -1,32 +0,0 @@ - subroutine tuispl (xa, ya, y2, x, y) -C -C Interpolate at point X using cubic splines. The array Y2 must have -C previously been computed by calling TUCSPL. Note that XA, YA, Y2 -C are two-element subarrays of the arrays with the same names elsewhere. -C Input and output are all double precision. -C This routine was copied with slight modifications from the SPLINT -C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and -C Vetterling. -C -C XA i: pair of independent-variable values -C YA i: pair of dependent-variable values -C Y2 i: second derivatives of YA at each point -C X i: value at which spline is to be computed -C Y o: interpolated value at X -C -CH Phil Hodge, 14-Apr-1988 Subroutine copied from Numerical Recipes SPLINT. -C - double precision xa(2), ya(2), y2(2), x, y -C-- - double precision h, a, b - - h = xa(2) - xa(1) - - a = (xa(2) - x) / h - b = (x - xa(1)) / h - y = a * ya(1) + b * ya(2) + - + ((a**3 - a) * y2(1) + (b**3 - b) * y2(2)) - + * h * h / 6. - - return - end diff --git a/pkg/utilities/nttools/trebin/tuispl.x b/pkg/utilities/nttools/trebin/tuispl.x new file mode 100644 index 000000000..32929e51b --- /dev/null +++ b/pkg/utilities/nttools/trebin/tuispl.x @@ -0,0 +1,24 @@ +# Interpolate at point X using cubic splines. The array K must have +# previously been computed by calling TUCSPL. Note that XA, YA, K +# are two-element subarrays of the arrays with the same names elsewhere. +# Input and output are all double precision. +# +# Copyright (c) 2017 Ole Streicher +# +# Reference: De Boor, Carl, et al. A practical guide to splines. +# Vol. 27. New York: Springer-Verlag, 1978. + +procedure tuispl (xa, ya, k, x, y) +double xa[2] # pair of independent-variable values +double ya[2] # pair of dependent-variable values +double k[2] # derivatives of YA at each point +double x # value at which spline is to be computed +double y # interpolated value at X + +double a, b, t +begin + a = k[1] * (xa[2]-xa[1]) - (ya[2]-ya[1]) + b = -k[2] * (xa[2]-xa[1]) + (ya[2]-ya[1]) + t = (x - xa[1]) / (xa[2] - xa[1]) + y = (1 - t)*ya[1] + t*ya[2] + t*(1-t)*(a*(1-t) + b*t) +end diff --git a/pkg/utilities/nttools/trebin/tuival.x b/pkg/utilities/nttools/trebin/tuival.x index 20354a5d5..0518e0da8 100644 --- a/pkg/utilities/nttools/trebin/tuival.x +++ b/pkg/utilities/nttools/trebin/tuival.x @@ -35,7 +35,6 @@ double y # o: the interpolated value #-- double x # independent variable double x1, x2 # beginning and end of output pixel -double dy # an estimate of the error of interpolation int klo, khi # indexes within xa that bracket x int npts # number of points in interval @@ -97,8 +96,8 @@ begin call tuhunt (xa, n, x, klo) klo = max (klo, 2) klo = min (klo, n-2) - # Pass tuiep3 the four points at klo-1, klo, klo+1, klo+2. - call tuiep3 (xa[klo-1], ya[klo-1], x, y, dy) + # Pass tuiepn the four points at klo-1, klo, klo+1, klo+2. + call tuiepn (4, xa[klo-1], ya[klo-1], x, y) case I_SPLINE: