1 files changed, 140 insertions, 0 deletions
diff --git a/pkg/utilities/nttools/stxtools/savgol.x b/pkg/utilities/nttools/stxtools/savgol.x
new file mode 100644
index 00000000..48958ede
--- /dev/null
+++ b/pkg/utilities/nttools/stxtools/savgol.x
@@ -0,0 +1,140 @@
+# Define some memory management.
+define	A		Memd[a+((($2)-1)*(m+1))+($1)-1]
+define	B		Memd[b+($1)-1]
+
+#---------------------------------------------------------------------------
+.help savgol Jun93 source
+.ih
+NAME
+savgol -- Create a kernel for Savitzky-Golay smoothing.
+.ih
+USAGE
+call savgol (c, np, nl, nr, ld, m)
+.ih
+ARGUMENTS
+.ls c (O: double[np])
+The smoothing kernel in "wrap-around" order.  See discussion for
+details.
+.le
+.ls np (I: int)
+The number of points allocated in the array represented by the "c"
+argument.
+.le
+.ls nl (I: int)
+Size of the kernel "to the left" of the central point.  See discussion
+for more details.
+.le
+.ls nr (I: int)
+Size of the kernel "to the right" of the central point. See discussion
+for more details.
+.le
+.ls ld (I: int)
+Order of the derivative desired.  Should be 0 for smoothing, higher
+for smoothed versions of the specified derivative.
+.le
+.ls m (I: int)
+Order of the smoothing polynomial.  Should be 0 or 1 for standard
+"boxcar" or "moving window" averaging.
+.le
+.ih
+DISCUSSION
+For an introduction to Savitzky-Golay filtering, see:
+
+.nf
+	Press, Teukolsky, Vetterling, & Falnnery, "Numeric Recipies:
+        The Art of Scientifitc Computing, Second Edition", Cambridge,
+	1992.
+.fi
+
+This routine returns the set of Savitzky-Golay smoothing coefficients
+given the size, order of smoothing polynomial, and derivative to
+return.  The coefficients are returned in "wrap-around" order.  Thus,
+if the smoothing coefficients are C[-nl]...C[0]...C[nr], they are
+returned in the array, c[i], as follows:
+
+.nf
+	c[1], c[2],  c[3], ..., c[nl+1],c[nl+2],...,c[np-1],c[np]
+
+	C[0], C[-1], C[-2],..., C[-nl], C[nr],  ...,C[2],   C[1]
+.fi
+
+A code fragment to transform the array c[i] to the orginal order,
+k[i], is: 
+
+.nf
+        do i = 1, nl+1 
+            k[i] = c[nl+2-i]
+        do i = 1, nr
+            k[nl+1+i] = c[np+1-i]
+.fi
+
+Array k[i], is now suitable for routines such as the IRAF VOPS routine
+acnvrd.
+.endhelp
+#---------------------------------------------------------------------------
+procedure savgol (c, np, nl, nr, ld, m)
+
+double	c[np]			# O:  The kernel.
+int	np			# I:  Size of the smoothing kernel.
+int	nl			# I:  Points to the left of center.
+int	nr			# I:  Points to the right of center.
+int	ld			# I:  Order of derivative to return.
+int	m			# I:  Order of the smoothing polynomial.
+
+int	imj, ipj, j, k, kk, mm, ix 
+double	d, fac, sum
+pointer	indx, a, b, sp
+int	shifti()
+
+begin
+	call smark (sp)
+	# Check input parameters.
+	if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m)
+	    call error (1, "savgol: invalid inputs")
+
+	# Allocate memory.
+	call salloc (indx, m+1, TY_INT)
+	call salloc (a, (m+1)**2, TY_DOUBLE)
+	call salloc (b, m+1, TY_DOUBLE)
+
+	# Do it.
+	ipj = shifti (m, 1)
+	do ipj = 0, shifti (m, 1) {
+	    if (ipj != 0)
+		sum = 0.d0
+	    else
+		sum = 1.d0
+	    do k = 1, nr
+		sum = sum + k**ipj
+	    do k = 1, nl
+		sum = sum + (-k)**ipj
+	    mm = min (ipj, 2*m-ipj)
+	    do imj = -mm, mm, 2
+		A(1+(ipj-imj)/2,1+(ipj+imj)/2) = sum
+	}
+	call ludcmd (Memd[a], m+1, m+1, Memi[indx], d, ix)
+	if (ix != OK)
+	    call error (1, "savgol: singular matrix")
+	do j = 1, m+1
+	    B(j) = 0.d0
+	B(ld+1) = 1.d0;
+	call lubksd (Memd[a], m+1, m+1, Memi[indx], Memd[b])
+	do kk = 1, np
+	    c[kk] = 0.d0
+	do k = -nl, nr {
+	    sum = B(1)
+	    fac = 1.d0
+	    do mm = 1, m {
+		fac = fac * k
+		sum = sum + B(mm+1) * fac
+	    }
+	    kk = mod (np - k, np) + 1
+	    c[kk] = sum
+	}
+
+	# That's all folks.
+	call sfree (sp)
+end
+#---------------------------------------------------------------------------
+# End of savgol
+#---------------------------------------------------------------------------