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

@@ -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'

pkg/utilities/nttools/trebin/mkpkg

@@ -11,16 +11,16 @@ libpkg.a:
 	tnamgio.x
 	tnaminit.x
 	trebin.x		<error.h> <tbset.h>
-	tucspl.f
+	tucspl.x
 	tudcol.x		<tbset.h>
 	tugcol.x		<error.h> <tbset.h>
 	tugetput.x		<tbset.h>
-	tuhunt.f
-	tuiep3.f
+	tuhunt.x
+	tuiepn.x
 	tuifit.x		trebin.h
 	tuinterp.x		<error.h> <tbset.h>
 	tuiset.x		trebin.h
-	tuispl.f
+	tuispl.x
 	tuival.x		trebin.h
 	tutrim.x
 	tuxget.x		<tbset.h>

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
+

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