1 files changed, 81 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/mssmooth.x b/noao/twodspec/multispec/mssmooth.x
new file mode 100644
index 00000000..be7e01ca
--- /dev/null
+++ b/noao/twodspec/multispec/mssmooth.x
@@ -0,0 +1,81 @@
+include	<math/curfit.h>
+
+# MS_SMOOTH -- Smooth MULTISPEC parameters with the CURFIT package.
+# MS_SET_SMOOTH -- Initialize and define function for smoothing.
+# MS_FREE_SMOOTH -- Free allocated memory from smoothing.
+
+# This procedure is numerical and does not depend on the MULTISPEC
+# package.
+
+procedure ms_smooth (x, y)
+
+real	x[ARB]				# Array of x values
+real	y[ARB]				# Array of y values
+int	curve_type			# Curfit function
+int	order				# Order of function
+real	xmin				# Minimum x value
+real	xmax				# Maximum x value
+int	npoints				# Number of points in fits
+
+int	i, npts, ier
+real	xmn, xmx
+pointer	cv, w
+
+real	cveval()
+
+data	cv/NULL/, w/NULL/
+
+begin
+	# Check for a valid curfit pointer.
+	if (cv == NULL)
+	    call error (0, "param_smooth: Undefined smoothing function")
+	
+	# Zero and fit the data with uniform weights.
+	call cvzero (cv)
+	# call cvfit (cv, x, y, Memr[w], npts, WTS_UNIFORM, ier)
+
+	# Accumulate points and check for out of bounds points.
+	do i = 1, npts
+	    if ((x[i] >= xmn) && (x[i] <= xmx))
+		call cvaccum (cv, x[i], y[i], Memr[w+i-1], WTS_UNIFORM)
+	call cvsolve (cv, ier)
+
+	if (ier != OK)
+	    call error (0, "param_smooth: Error in function fit")
+
+	# Evaluate fit placing fit values back in y array.
+	# call cvvector (cv, x, y, npts)
+	do i = 1, npts
+	    if ((x[i] >= xmn) && (x[i] <= xmx))
+		y[i] = cveval (cv, x[i])
+
+	return
+
+entry ms_set_smooth (xmin, xmax, npoints)
+
+	# Set or reset curfit data structure and allocate memory for weights.
+	if (cv != NULL)
+	    call cvfree (cv)
+	if (w == NULL)
+	    call malloc (w, npoints, TY_REAL)
+
+	# Determine curve_type and order.
+	call clgcurfit ("function", "order", curve_type, order)
+
+	# Initialize curfit data structure and record number of points.
+	xmn = xmin
+	xmx = xmax
+	call cvinit (cv, curve_type, order, xmn, xmx)
+	npts = npoints
+
+	return
+
+entry ms_free_smooth ()
+
+	# Free allocated memory.
+	if (cv != NULL)
+	    call cvfree (cv)
+	if (w != NULL)
+	    call mfree (w, TY_REAL)
+
+end