1 files changed, 158 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/t_fitfunc.x b/noao/twodspec/multispec/t_fitfunc.x
new file mode 100644
index 00000000..9f6209ad
--- /dev/null
+++ b/noao/twodspec/multispec/t_fitfunc.x
@@ -0,0 +1,158 @@
+include <math/curfit.h>
+include	"ms.h"
+
+# T_FIT_FUNCTION -- Fit a function to selected spectra parameters.
+#
+# A function is fit to the parameter values determined at the sample
+# lines for selected spectra.  The function coefficients are stored in
+# the database and the fitted values replace the original values at
+# the sample lines.  The type of function, the parameter to be fitted,
+# the sample lines used in the fit, and the spectra to be fitted
+# are all selected by the user.
+
+procedure t_fit_function()
+
+char	image[SZ_FNAME]			# Image affected
+char	parameter[SZ_LINE]		# Parameter to be fit
+int	function			# Type of fitting function
+int	order				# Order of the fitting function
+int	spectra[3, MAX_RANGES]		# Spectra to be fitted
+pointer	samples				# Sample lines to be fitted.
+
+int	i, param_id, nsamples
+pointer	ms, sp
+
+int	ms_db_id(), clgranges(), get_sample_lines()
+pointer	msmap()
+
+begin
+	# Access database and determine parameter to be fit and the
+	# fitting function and order.
+
+	call clgstr ("image", image, SZ_FNAME)
+	ms = msmap (image, READ_WRITE, 0)
+	call clgstr ("parameter", parameter, SZ_LINE)
+	param_id = ms_db_id (ms, parameter)
+	call clgcurfit ("function", "order", function, order)
+
+	# Get the image lines to be used in the fit and convert to sample
+	# lines.  Get the spectra to be fit.
+
+	i = clgranges ("lines", 1, MS_LEN(ms, 2), spectra, MAX_RANGES)
+	call smark (sp)
+	call salloc (samples, MS_NSAMPLES(ms), TY_INT)
+	nsamples = get_sample_lines (ms, spectra, Memi[samples])
+	i = clgranges ("spectra", 1, MS_NSPECTRA(ms), spectra,
+	    MAX_RANGES)
+
+	# Fit the parameters for each spectrum, store the fits in the database,
+	# and substitute the fitted values for the parameter values at all
+	# the sample lines.
+
+	call fit_function (ms, Memi[samples], nsamples, spectra, param_id,
+	    function, order)
+
+	# Finish up.
+	call msphdr (ms)
+	call msunmap (ms)
+	call sfree (sp)
+end
+
+
+
+# FIT_FUNCTION -- Fit a function to the parameter data.
+#
+# If the fit coefficients for the specified parameter are not in
+# the database then the database entry is created.
+
+procedure fit_function (ms, lines, nlines, spectra, param_id, function, order)
+
+pointer	ms				# MULTISPEC data structure
+int	lines[nlines]			# Sample lines to be used
+int	nlines				# Number of sample lines
+int	spectra[ARB]			# Spectra to be fitted
+int	param_id			# Parameter being fit
+int	function			# Function to be fit
+int	order				# Order of the function
+
+char	comment[SZ_LINE]
+int	i, spectrum, fit_id, ier
+real	x, wt
+
+int	ms_fit_id(), get_next_number()
+real	cveval()
+bool	dbaccess()
+
+begin
+	# Determine the MULTISPEC fit id from the parameter id.
+	fit_id = ms_fit_id (param_id)
+	if (fit_id == ERR)
+	    call error (MS_ERROR, "Unknown fit identifier")
+
+	# Enter the fit records in the database if necessary.
+	if (!dbaccess (MS_DB(ms), NAME(ms, fit_id)))
+	    call dbenter (MS_DB(ms), NAME(ms, fit_id),
+		(7 + MS_NSAMPLES(ms)) * SZ_REAL, MS_NSPECTRA(ms))
+
+	# Allocate memory for the curfit data structures pointers.
+	if (MS_DATA(ms, fit_id) == NULL)
+	    call malloc (MS_DATA(ms, fit_id), MS_NSPECTRA(ms), TY_INT)
+
+	# Initialize the curfit data structures.
+	# If the order is INDEF then use maximum order assuming no INDEF points.
+	spectrum = 0
+	while (get_next_number (spectra, spectrum) != EOF) {
+	    if (IS_INDEFI (order)) {
+	        switch (function) {
+	        case LEGENDRE, CHEBYSHEV:
+		    order =  nlines
+	        case SPLINE3:
+		    order = nlines - 3
+	        }
+	    }
+	    call cvinit (CV(ms, fit_id, spectrum), function, order, 1.,
+		real (MS_LEN(ms, 2)))
+	}
+
+	# Accumulate the parameter values.
+	do i = 1, nlines {
+	    x = LINE(ms, lines[i])
+	    call msgparam (ms, param_id, lines[i])
+
+	    spectrum = 0
+	    while (get_next_number (spectra, spectrum) != EOF)
+		call cvaccum (CV(ms, fit_id, spectrum), x,
+		    PARAMETER(ms, param_id, spectrum), wt, WTS_UNIFORM)
+	}
+
+	# Compute and write the fit coeffients to the database.
+
+	spectrum = 0
+	while (get_next_number (spectra, spectrum) != EOF) {
+	    call cvsolve (CV(ms, fit_id, spectrum), ier)
+	    if (ier == NO_DEG_FREEDOM)
+		call error (MS_ERROR, "Error fitting parameters")
+	    call mspfit (ms, fit_id, spectrum)
+	}
+
+	# For each sample line and each selected spectrum replace the
+	# selected parameter value with the fit evaluation.
+
+	do i = 1, MS_NSAMPLES(ms) {
+	    x = LINE(ms, i)
+	    call msgparam (ms, param_id, i)
+
+	    spectrum = 0
+	    while (get_next_number (spectra, spectrum) != EOF)
+		PARAMETER(ms, param_id, spectrum) =
+		    cveval (CV(ms, fit_id, spectrum), x)
+
+	    call mspparm (ms, param_id, i)
+	}
+
+	# Add a comment to the database comments.
+
+	call sprintf (comment, SZ_LINE, "Fit a function to parameter %s.")
+	    call pargstr (NAME(ms, param_id))
+	call history (ms, comment)
+end