Initial commit

1 files changed, 669 insertions, 0 deletions

diff --git a/noao/twodspec/apextract/aptrace.x b/noao/twodspec/apextract/aptrace.x
new file mode 100644
index 00000000..c38af01c
--- /dev/null
+++ b/noao/twodspec/apextract/aptrace.x
@@ -0,0 +1,669 @@
+include	<imhdr.h>
+include	<math/curfit.h>
+include	<pkg/center1d.h>
+include	<pkg/gtools.h>
+include	"apertures.h"
+
+define	MAXBUF	100000		# Column buffer size
+
+
+# AP_TRACE -- Trace features in a two dimensional image.
+#
+# Given an image pointer, the starting dispersion position, and a set
+# of apertures defining the centers of features, trace the feature
+# centers to other dispersion positions and fit a curve to the positions.
+# The user specifies the dispersion step size, the number of dispersion
+# lines to sum, and parameters for the feature centering function
+# fitting.
+
+procedure ap_trace (image, line, aps, naps, apedit)
+
+char	image[SZ_FNAME]		# Image name
+int	line			# Starting dispersion position
+pointer	aps[ARB]		# Apertures
+int	naps			# Number of apertures
+int	apedit			# Called from APEDIT?
+
+int	step			# Tracing step
+int	nsum			# Number of dispersion lines to sum
+int	nlost			# Number of steps lost before quitting
+real	cradius			# Centering radius
+real	cwidth			# Centering width
+real	cthreshold		# Detection threshold for centering
+
+int	i, na, dispaxis, apaxis
+real	center
+pointer	im, ic, ic1, sp, str
+data	ic1 /NULL/
+
+int	apgeti()
+real	apgetr()
+bool	clgetb(), ap_answer()
+pointer	ap_immap()
+
+errchk	ap_immap, ic_open, ap_ltrace, ap_ctrace, ap_default
+
+common	/apt_com/ ic
+
+begin
+	na = 0
+	do i = 1, naps
+	    if (AP_SELECT(aps[i]) == YES)
+		na = na + 1
+	if (naps > 0 && na == 0)
+	    return
+
+	# Query user.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	if (apedit == NO) {
+	    call sprintf (Memc[str], SZ_LINE, "Trace apertures for %s?")
+	        call pargstr (image)
+	    if (!ap_answer ("anstrace", Memc[str])) {
+	        call sfree (sp)
+	        return
+	    }
+
+	    call sprintf (Memc[str], SZ_LINE,
+	        "Fit traced positions for %s interactively?")
+	        call pargstr (image)
+	    if (ap_answer ("ansfittrace", Memc[str])) {
+	        call apgstr ("ansfittrace", Memc[str], SZ_LINE)
+	        call appstr ("ansfittrace1", Memc[str])
+	    } else
+	        call appstr ("ansfittrace1", "NO")
+
+	    if (clgetb ("verbose"))
+	        call printf ("Tracing apertures ...\n")
+	}
+
+	# Tracing parameters
+	step = apgeti ("t_step")
+	nsum = max (1, abs (apgeti ("t_nsum")))
+	nlost = apgeti ("t_nlost")
+	if (ic == NULL || ic1 == NULL) {
+	    call ic_open (ic)
+	    ic1 = ic
+	    call apgstr ("t_function", Memc[str], SZ_LINE)
+	    call ic_pstr (ic, "function", Memc[str])
+	    call ic_puti (ic, "order", apgeti ("t_order"))
+	    call apgstr ("t_sample", Memc[str], SZ_LINE)
+	    call ic_pstr (ic, "sample", Memc[str])
+	    call ic_puti (ic, "naverage", apgeti ("t_naverage"))
+	    call ic_puti (ic, "niterate", apgeti ("t_niterate"))
+	    call ic_putr (ic, "low", apgetr ("t_low_reject"))
+	    call ic_putr (ic, "high", apgetr ("t_high_reject"))
+	    call ic_putr (ic, "grow", apgetr ("t_grow"))
+	}
+
+	im = ap_immap (image, apaxis, dispaxis)
+
+	# If no apertures are defined default to the center of the image.
+	if (naps == 0) {
+	    naps = 1
+	    center = IM_LEN (im, apaxis) / 2.
+	    call ap_default (im, 1, 1, apaxis, center, real (line),
+		aps[naps])
+	    call sprintf (Memc[str], SZ_LINE,
+	       "TRACE - Default aperture defined centered on %s")
+	       call pargstr (image)
+	    call ap_log (Memc[str], YES, NO, YES)
+	}
+
+	# Centering parameters
+	cwidth = apgetr ("t_width")
+	cradius = apgetr ("radius")
+	cthreshold = apgetr ("threshold")
+
+	switch (dispaxis) {
+	case 1:
+	    call ap_ctrace (image, im, ic, line, step, nsum, nlost, cradius,
+		cwidth, cthreshold, aps, naps)
+	case 2:
+	    call ap_ltrace (image, im, ic, line, step, nsum, nlost, cradius,
+		cwidth, cthreshold, aps, naps)
+	}
+
+	# Log the tracing and write the traced apertures to the database.
+
+	call sprintf (Memc[str], SZ_LINE,
+	    "TRACE - %d apertures traced in %s.")
+       	    call pargi (na)
+       	    call pargstr (image)
+	if (apedit == NO)
+	    call ap_log (Memc[str], YES, YES, NO)
+	else
+	    call ap_log (Memc[str], YES, NO, NO)
+
+	call appstr ("ansdbwrite1", "yes")
+
+	call imunmap (im)
+	call sfree (sp)
+end
+
+
+procedure ap_trfree ()
+
+pointer	ic
+common	/apt_com/ ic
+
+begin
+	call ic_closer (ic)
+end
+
+
+# AP_CTRACE -- Trace feature positions for aperture axis 2.
+
+procedure ap_ctrace (image, im, ic, start, step, nsum, nlost, cradius, cwidth,
+    threshold, aps, naps)
+
+char	image[ARB]		# Image to be traced.
+pointer	im			# IMIO pointer
+pointer	ic			# ICFIT pointer
+int	start			# Starting column
+int	step			# Tracing step size
+int	nsum			# Number of lines or columns to sum
+int	nlost			# Number of steps lost before quiting
+real	cradius			# Centering radius
+real	cwidth			# Centering width
+real	threshold		# Detection threshold for centering
+pointer	aps[ARB]		# Apertures
+int	naps			# Number of apertures
+
+int	nlines, col, col1, col2, line1, line2
+int	i, j, n, nx, ny, ntrace, istart, lost, fd
+real	yc, yc1
+pointer	co, data, sp, str, x, y, wts, gp, gt
+
+real	center1d(), ap_cveval()
+bool	ap_answer()
+pointer comap(), gt_init()
+
+errchk	ap_cveval, xt_csum, xt_csumb, center1d, icg_fit, ic_fit
+errchk	ap_gopen, ap_popen
+
+begin
+	# Set up column access buffering.
+
+	co = comap (im, MAXBUF)
+
+	# Determine the number of lines to be traced and allocate memory.
+
+	nx = IM_LEN(im, 1)
+	ny = IM_LEN(im, 2)
+	if (IS_INDEFI (start))
+	    start = nx / 2
+	nlines = 5 * cwidth
+	istart = (start - 1) / step + 1
+	ntrace = istart + (nx - start) / step
+
+	# Allocate memory for the traced positions and the weights for fitting.
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (x, ntrace, TY_REAL)
+	call salloc (y, ntrace, TY_REAL)
+	call salloc (wts, ntrace, TY_REAL)
+	call aclrr (Memr[y], ntrace)
+	data = NULL
+
+	# Initialize the ICFIT limits and the GTOOLS parameters.
+	# Set initial interactive flag.
+
+	call ic_putr (ic, "xmin", 1.)
+	call ic_putr (ic, "xmax", real (nx))
+	call ic_pstr (ic, "xlabel", "Column")
+	call ic_pstr (ic, "ylabel", "Line")
+
+	gt = gt_init()
+	call gt_setr (gt, GTXMIN, 1. - step / 2)
+	call gt_setr (gt, GTXMAX, real (nx + step / 2))
+
+	# Trace each feature.
+
+	line1 = 0
+	line2 = 0
+	do j = 1, naps {
+	    if (AP_SELECT(aps[j]) == NO)
+		next
+
+	    # Trace from the starting column to the last column while the
+	    # position is not INDEF.
+
+	    lost = 0
+	    yc = AP_CEN(aps[j], 2) + ap_cveval (AP_CV(aps[j]), real (start))
+	    do i = istart, ntrace {
+		Memr[y+i-1] = INDEF
+		if (lost < nlost) {
+		    # Update the scrolling buffer if the feature center is less
+		    # than cwidth from the edge of the buffer.
+		    if (((yc-line1) < cwidth) || ((line2-yc) < cwidth)) {
+		        line1 = max (1, int (yc + .5 - nlines / 2))
+		        line2 = min (ny, line1 + nlines - 1)
+		        line1 = max (1, line2 - nlines + 1)
+		    }
+
+		    # Sum columns to form the 1D vector for centering.
+
+	            col = start + (i - istart) * step
+		    col1 = max (1, col - nsum / 2)
+		    col2 = min (nx, col1 + nsum - 1)
+		    col1 = max (1, col2 - nsum + 1)
+
+		    # If columns in the sum overlap then use buffering.
+
+		    if (step < nsum)
+		        call xt_csumb (co, col1, col2, line1, line2, data)
+		    else
+		        call xt_csum (co, col1, col2, line1, line2, data)
+
+		    # Center the feature for the new column using the previous
+		    # center as the starting point.  Convert to position
+		    # relative to the start of the data buffer for centering
+		    # and then convert back to position relative to the
+		    # edge of the image.
+
+		    yc1 = center1d (yc-line1+1, Memr[data], line2-line1+1,
+			cwidth, EMISSION, cradius, threshold)
+
+		    if (!IS_INDEF (yc1)) {
+			lost = 0
+		        yc = yc1 + line1 - 1
+			Memr[y+i-1] = yc
+			if (IS_INDEF (Memr[y+i-2])) {
+	    		    call sprintf (Memc[str], SZ_LINE,
+		   "TRACE - Trace of aperture %d in %s recovered at column %d.")
+			        call pargi (AP_ID(aps[j]))
+			        call pargstr (image)
+			        call pargi (col)
+			    call ap_log (Memc[str], YES, NO, YES)
+			}
+		    } else {
+			lost = lost + 1
+	    		call sprintf (Memc[str], SZ_LINE,
+		    "TRACE - Trace of aperture %d in %s lost at column %d.")
+			    call pargi (AP_ID(aps[j]))
+			    call pargstr (image)
+			    call pargi (col)
+			call ap_log (Memc[str], YES, NO, YES)
+		    }
+		}
+	    }
+
+	    # Trace from the starting column to the first column while the
+	    # position is not INDEF.
+
+	    lost = 0
+	    yc = AP_CEN(aps[j], 2) + ap_cveval (AP_CV(aps[j]), real (start))
+	    do i = istart - 1, 1, -1 {
+		Memr[y+i-1] = INDEF
+		if (lost < nlost) {
+		    # Update the scrolling buffer if the feature center is less
+		    # than cwidth from the edge of the buffer.
+
+		    if (((yc-line1) < cwidth) || ((line2-yc) < cwidth)) {
+		        line1 = max (1, int (yc + .5 - nlines / 2))
+		        line2 = min (ny, line1 + nlines - 1)
+		        line1 = max (1, line2 - nlines + 1)
+		    }
+
+		    # Sum columns to form the 1D vector for centering.
+
+	            col = start + (i - istart) * step
+		    col1 = max (1, col - nsum / 2)
+		    col2 = min (nx, col1 + nsum - 1)
+		    col1 = max (1, col2 - nsum + 1)
+
+		    # If columns in the sum overlap then use buffering.
+
+		    if (step < nsum)
+		        call xt_csumb (co, col1, col2, line1, line2, data)
+		    else
+		        call xt_csum (co, col1, col2, line1, line2, data)
+
+		    # Center the feature for the new column using the previous
+		    # center as the starting point.  Convert to position
+		    # relative to the start of the data buffer for centering
+		    # and then convert back to position relative to the
+		    # edge of the image.
+
+		    yc1 = center1d (yc-line1+1, Memr[data], line2-line1+1,
+			cwidth, EMISSION, cradius, threshold)
+
+		    if (!IS_INDEF (yc1)) {
+			lost = 0
+		        yc = yc1 + line1 - 1
+			Memr[y+i-1] = yc
+			if (IS_INDEF (Memr[y+i])) {
+	    		    call sprintf (Memc[str], SZ_LINE,
+		   "TRACE - Trace of aperture %d in %s recovered at column %d.")
+			        call pargi (AP_ID(aps[j]))
+			        call pargstr (image)
+			        call pargi ((i - 1) * step + 1)
+			    call ap_log (Memc[str], YES, NO, YES)
+			}
+		    } else {
+			lost = lost + 1
+	    		call sprintf (Memc[str], SZ_LINE,
+		    "TRACE - Trace of aperture %d in %s lost at column %d.")
+			    call pargi (AP_ID(aps[j]))
+			    call pargstr (image)
+			    call pargi ((i - 1) * step + 1)
+			call ap_log (Memc[str], YES, NO, YES)
+		    }
+		}
+	    }
+
+	    # Order the traced points and exclude INDEF positions.
+
+	    n = 0
+	    do i = 1, ntrace {
+		if (IS_INDEF (Memr[y+i-1]))
+		    next
+		n = n + 1
+		Memr[x+n-1] = start + (i - istart) * step
+		Memr[y+n-1] = Memr[y+i-1]
+		Memr[wts+n-1] = 1.
+	    }
+
+	    # If all positions are INDEF print a message and go on to the next
+	    # aperture.
+
+	    if (n < 2) {
+		call eprintf (
+		    "Not enough points traced for aperture %d of %s\n")
+		    call pargi (AP_ID(aps[j]))
+		    call pargstr (image)
+		next
+	    }
+		    
+	    # Fit a curve to the traced positions and graph the result.
+
+	    call sprintf (Memc[str], SZ_LINE, "Aperture %d of %s")
+	        call pargi (AP_ID(aps[j]))
+		call pargstr (image)
+	    call gt_sets (gt, GTTITLE, Memc[str])
+
+	    call sprintf (Memc[str], SZ_LINE,
+		"Fit curve to aperture %d of %s interactively")
+		call pargi (AP_ID(aps[j]))
+		call pargstr (image)
+	    if (ap_answer ("ansfittrace1", Memc[str])) {
+		call ap_gopen (gp)
+	        call icg_fit (ic, gp, "gcur", gt,
+		    AP_CV(aps[j]), Memr[x], Memr[y], Memr[wts], n)
+	    } else
+	        call ic_fit (ic, AP_CV(aps[j]), Memr[x], Memr[y], Memr[wts], n,
+		    YES, YES, YES, YES)
+
+	    call ap_popen (gp, fd, "trace")
+	    if (gp != NULL) {
+		call icg_graphr (ic, gp, gt, AP_CV(aps[j]),
+		    Memr[x], Memr[y], Memr[wts], n)
+		call ap_pclose (gp, fd)
+	    }
+
+	    call asubkr (Memr[y], AP_CEN(aps[j], 2), Memr[y], n)
+	    call ic_fit (ic, AP_CV(aps[j]), Memr[x], Memr[y], Memr[wts], n,
+		YES, YES, YES, YES)
+	}
+
+	# Free allocated memory.
+
+	call gt_free (gt)
+	call mfree (data, TY_REAL)
+	call counmap (co)
+	call sfree (sp)
+end
+
+
+# AP_LTRACE -- Trace feature positions for aperture axis 1.
+
+procedure ap_ltrace (image, im, ic, start, step, nsum, nlost, cradius, cwidth,
+    threshold, aps, naps)
+
+char	image[ARB]		# Image to be traced
+pointer	im			# IMIO pointer
+pointer	ic			# ICFIT pointer
+int	start			# Starting line
+int	step			# Tracing step size
+int	nsum			# Number of lines or columns to sum
+int	nlost			# Number of steps lost before quiting
+real	cradius			# Centering radius
+real	cwidth			# Centering width
+real	threshold		# Detection threshold for centering
+pointer	aps[ARB]		# Apertures
+int	naps			# Number of apertures
+
+real	xc1
+int	i, j, n, nx, ny, ntrace, istart, line, line1, line2, fd
+pointer	data, sp, str, x, y, wts, xc, lost, x1, x2, gp, gt
+
+real	center1d(), ap_cveval()
+bool	ap_answer()
+pointer	gt_init()
+
+errchk	ap_cveval, xt_lsum, xt_lsumb, center1d, icg_fit, ic_fit
+errchk	ap_gopen, ap_popen
+
+begin
+	# Determine the number of lines to be traced and allocate memory.
+
+	nx = IM_LEN(im, 1)
+	ny = IM_LEN(im, 2)
+	if (IS_INDEFI (start))
+	    start = ny / 2
+
+	istart = (start - 1) / step + 1
+	ntrace = istart + (ny - start) / step
+
+	# Allocate memory for the traced positions and the weights for fitting.
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (x, ntrace * naps, TY_REAL)
+	call salloc (y, ntrace, TY_REAL)
+	call salloc (wts, ntrace, TY_REAL)
+	call salloc (xc, naps, TY_REAL)
+	call salloc (lost, naps, TY_INT)
+	call aclrr ( Memr[x], ntrace * naps)
+	data = NULL
+
+	# Set the dispersion lines to be traced.
+
+	do i = 1, ntrace
+	    Memr[y+i-1] = start + (i - istart) * step
+
+	# Trace from the starting line to the last line.
+
+	x1 = x + istart - 1
+	do i = 1, naps {
+	    if (AP_SELECT(aps[i]) == NO)
+		next
+	    Memr[xc+i-1] = AP_CEN(aps[i], 1) +
+		ap_cveval (AP_CV(aps[i]), real (start))
+	    Memi[lost+i-1] =  0
+	}
+
+	do i = istart, ntrace {
+	    line = Memr[y+i-1]
+	    line1 = max (1, line - nsum / 2)
+	    line2 = min (ny, line1 + nsum - 1)
+	    line1 = max (1, line2 - nsum + 1)
+
+	    # If the sums overlap use buffering.
+
+	    if (step < nsum)
+	        call xt_lsumb (im, 1, nx, line1, line2, data)
+	    else
+	        call xt_lsum (im, 1, nx, line1, line2, data)
+
+	    do j = 1, naps {
+		if (AP_SELECT(aps[j]) == NO)
+		    next
+		x2 = x1 + (j - 1) * ntrace
+		Memr[x2] = INDEF
+		if (Memi[lost+j-1] < nlost) {
+		    xc1 = center1d (Memr[xc+j-1], Memr[data], nx,
+			cwidth, EMISSION, cradius, threshold)
+		    if (IS_INDEF(xc1)) {
+			Memi[lost+j-1] = Memi[lost+j-1] + 1
+			call sprintf (Memc[str], SZ_LINE,
+			"TRACE - Trace of aperture %d in %s lost at line %d.")
+			    call pargi (AP_ID(aps[j]))
+			    call pargstr (image)
+			    call pargi (line)
+			call ap_log (Memc[str], YES, NO, YES)
+		    } else {
+			Memi[lost+j-1] = 0
+			Memr[xc+j-1] = xc1
+			Memr[x2] = xc1
+			if (IS_INDEF (Memr[x2-1])) {
+			    call sprintf (Memc[str], SZ_LINE,
+		    "TRACE - Trace of aperture %d in %s recovered at line %d.")
+				call pargi (AP_ID(aps[j]))
+				call pargstr (image)
+				call pargi (line)
+			    call ap_log (Memc[str], YES, NO, YES)
+			}
+		    }
+		}
+	    }
+	    x1 = x1 + 1
+	}
+
+	# Trace from the starting line to the first line.
+
+	x1 = x + istart - 2
+	do i = 1, naps {
+	    if (AP_SELECT(aps[i]) == NO)
+		next
+	    Memr[xc+i-1] = AP_CEN(aps[i], 1) +
+		ap_cveval (AP_CV(aps[i]), real (start))
+	    Memi[lost+i-1] = 0
+	}
+
+	do i = istart - 1, 1, -1 {
+	    line = Memr[y+i-1]
+	    line1 = max (1, line - nsum / 2)
+	    line2 = min (ny, line1 + nsum - 1)
+	    line1 = max (1, line2 - nsum + 1)
+
+	    # If the sums overlap use buffering.
+
+	    if (step < nsum)
+	        call xt_lsumb (im, 1, nx, line1, line2, data)
+	    else
+	        call xt_lsum (im, 1, nx, line1, line2, data)
+
+	    do j = 1, naps {
+		if (AP_SELECT(aps[j]) == NO)
+		    next
+		x2 = x1 + (j - 1) * ntrace
+		Memr[x2] = INDEF
+		if (Memi[lost+j-1] < nlost) {
+		    xc1 = center1d (Memr[xc+j-1], Memr[data], nx,
+			cwidth, EMISSION, cradius, threshold)
+		    if (IS_INDEF(xc1)) {
+			Memi[lost+j-1] = Memi[lost+j-1] + 1
+			call sprintf (Memc[str], SZ_LINE,
+			    "TRACE - Trace of aperture %d in %s lost at line %d.")
+			    call pargi (AP_ID(aps[j]))
+			    call pargstr (image)
+			    call pargi (line)
+			call ap_log (Memc[str], YES, NO, YES)
+		    } else {
+			Memi[lost+j-1] = 0
+			Memr[xc+j-1] = xc1
+			Memr[x2] = xc1
+			if (IS_INDEF (Memr[x2+1])) {
+			    call sprintf (Memc[str], SZ_LINE,
+		    "TRACE - Trace of aperture %d in %s recovered at line %d.")
+				call pargi (AP_ID(aps[j]))
+				call pargstr (image)
+				call pargi (line)
+			    call ap_log (Memc[str], YES, NO, YES)
+			}
+		    }
+		}
+	    }
+	    x1 = x1 - 1
+	}
+
+	# Initialize the the GTOOLS parameters.
+	call ic_putr (ic, "xmin", 1.)
+	call ic_putr (ic, "xmax", real (ny))
+	call ic_pstr (ic, "xlabel", "Line")
+	call ic_pstr (ic, "ylabel", "Column")
+
+	gt = gt_init()
+	call gt_setr (gt, GTXMIN, 1. - step / 2)
+	call gt_setr (gt, GTXMAX, real (ny + step / 2))
+
+	do j = 1, naps {
+	    if (AP_SELECT(aps[j]) == NO)
+		next
+
+	    # Order the traced points and exclude INDEF positions.
+
+	    x1 = x + (j - 1) * ntrace
+	    n = 0
+
+	    do i = 1, ntrace {
+		if (IS_INDEF (Memr[x1+i-1]))
+		    next
+		n = n + 1
+		Memr[x1+n-1] = Memr[x1+i-1]
+		Memr[y+n-1] = start + (i - istart) * step
+		Memr[wts+n-1] = 1.
+	    }
+
+	    # If all positions are INDEF print a message and go on to the next
+	    # aperture.
+
+	    if (n < 2) {
+		call eprintf (
+		    "Not enough points traced for aperture %d of %s\n")
+		    call pargi (AP_ID(aps[j]))
+		    call pargstr (image)
+		next
+	    }
+
+	    # Fit a curve to the traced positions and graph the result.
+
+	    call sprintf (Memc[str], SZ_LINE, "Aperture %d of %s")
+	        call pargi (AP_ID(aps[j]))
+		call pargstr (image)
+	    call gt_sets (gt, GTTITLE, Memc[str])
+
+	    call sprintf (Memc[str], SZ_LINE,
+		"Fit curve to aperture %d of %s interactively")
+		call pargi (AP_ID(aps[j]))
+		call pargstr (image)
+	    if (ap_answer ("ansfittrace1", Memc[str])) {
+		call ap_gopen (gp)
+	        call icg_fit (ic, gp, "gcur", gt,
+		    AP_CV(aps[j]), Memr[y], Memr[x1], Memr[wts], n)
+	    } else
+	        call ic_fit (ic, AP_CV(aps[j]), Memr[y], Memr[x1], Memr[wts], n,
+		    YES, YES, YES, YES)
+
+	    call ap_popen (gp, fd, "trace")
+	    if (gp != NULL) {
+		call icg_graphr (ic, gp, gt, AP_CV(aps[j]),
+		    Memr[y], Memr[x1], Memr[wts], n)
+		call ap_pclose (gp, fd)
+	    }
+
+	    # Subtract the aperture center and refit offset curve.
+	    call asubkr (Memr[x1], AP_CEN(aps[j], 1), Memr[x1], n)
+	    call ic_fit (ic, AP_CV(aps[j]), Memr[y], Memr[x1], Memr[wts], n,
+		YES, YES, YES, YES)
+	}
+
+	# Free allocated memory.
+
+	call gt_free (gt)
+	call mfree (data, TY_REAL)
+	call sfree (sp)
+end