Initial commit

1 files changed, 737 insertions, 0 deletions

diff --git a/noao/twodspec/apextract/apfit.x b/noao/twodspec/apextract/apfit.x
new file mode 100644
index 00000000..67bf149d
--- /dev/null
+++ b/noao/twodspec/apextract/apfit.x
@@ -0,0 +1,737 @@
+include	<imhdr.h>
+include	<imset.h>
+include	<pkg/gtools.h>
+include	"apertures.h"
+
+
+# AP_FITSPEC -- Fit a spectrum by a smoothing function.
+
+procedure ap_fitspec (ap, in, spec, ny)
+
+pointer	ap			# Aperture (used for labels)
+pointer	in			# Input image (used for labels)
+real	spec[ny]		# spectrum
+int	ny			# Number of points in spectra
+
+int	i, fd, apaxis, clgeti()
+real	clgetr()
+pointer	sp, str, x, wts, cv, gp, gt, ic, ic1, gt_init()
+bool	ap_answer()
+data	ic1 /NULL/
+errchk	icg_fit, ic_fit
+
+common	/apn_com/ ic, gt
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (x, ny, TY_REAL)
+	call salloc (wts, ny, TY_REAL)
+
+	do i = 1, ny {
+	    Memr[x+i-1] = i
+	    Memr[wts+i-1] = 1
+	}
+
+	if (ic == NULL || ic1 == NULL) {
+	    call ic_open (ic)
+	    ic1 = ic
+	    call clgstr ("function", Memc[str], SZ_LINE)
+	    call ic_pstr (ic, "function", Memc[str])
+	    call ic_puti (ic, "order", clgeti ("order"))
+	    call clgstr ("sample", Memc[str], SZ_LINE)
+	    call ic_pstr (ic, "sample", Memc[str])
+	    call ic_puti (ic, "naverage", clgeti ("naverage"))
+	    call ic_puti (ic, "niterate", clgeti ("niterate"))
+	    call ic_putr (ic, "low", clgetr ("low_reject"))
+	    call ic_putr (ic, "high", clgetr ("high_reject"))
+	    call ic_putr (ic, "grow", clgetr ("grow"))
+	    call ic_pstr (ic, "ylabel", "")
+
+	    gt = gt_init()
+	}
+
+	call ic_putr (ic, "xmin", 1.)
+	call ic_putr (ic, "xmax", real (ny))
+	apaxis = AP_AXIS(ap)
+	switch (apaxis) {
+	case 1:
+	    call ic_pstr (ic, "xlabel", "Line")
+	case 2:
+	    call ic_pstr (ic, "xlabel", "Column")
+	}
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Fit spectrum by a smoothing function.
+	call sprintf (Memc[str], SZ_LINE,
+	    "%s: %s - Aperture %s")
+	    call pargstr (IM_HDRFILE(in))
+	    call pargstr (IM_TITLE(in))
+	    call pargi (AP_ID(ap))
+	call gt_sets (gt, GTTITLE, Memc[str])
+
+	# Query the user to fit the spectrum interactively.
+	call sprintf (Memc[str], SZ_LINE,
+	    "Fit spectrum for aperture %d for %s interactively?")
+	    call pargi (AP_ID(ap))
+	    call pargstr (IM_HDRFILE(in))
+	if (ap_answer ("ansfitspec1", Memc[str])) {
+	    call ap_gopen (gp)
+	    call icg_fit (ic, gp, "gcur", gt, cv, Memr[x], spec,
+		Memr[wts], ny)
+	    call amovkr (1., Memr[wts], ny)
+	} else
+	    call ic_fit (ic, cv, Memr[x], spec, Memr[wts], ny,
+		YES, YES, YES, YES)
+
+	# Make a graph to the plot log.
+	call ap_popen (gp, fd, "fitspec")
+	if (gp != NULL) {
+	    call icg_graphr (ic, gp, gt, cv, Memr[x], spec, Memr[wts], ny)
+	    call ap_pclose (gp, fd)
+	}
+
+	call cvvector (cv, Memr[x], spec, ny)
+	call cvfree (cv)
+end
+
+
+procedure ap_fitfree ()
+
+pointer	ic, gt
+common	/apn_com/ ic, gt
+
+begin
+	call ic_closer (ic)
+	call gt_free (gt)
+end
+
+
+# AP_LNORM -- Normalize the input line apertures by the norm spectra.
+
+procedure ap_lnorm (ap, out, gain, dbuf, nc, nl, c1, l1, spec, ny, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+int	ny			# Size of profile array
+int	ys			# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+bool	clgetb()		# Center normalize?
+real	threshold, clgetr()	# Division threshold
+
+int	i, ncols, nlines, ix1, ix2, iy, nsum
+real	cen, low, high, s, x1, x2, sum, ap_cveval(), asumr()
+pointer	cv, datain, dataout, imps2r(), impl2r()
+
+begin
+	threshold = clgetr ("threshold")
+
+	cen = AP_CEN(ap,1)
+	low = AP_CEN(ap,1) + AP_LOW(ap,1)
+	high = AP_CEN(ap,1) + AP_HIGH(ap,1)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	# Normalize by the aperture width and apply threshold.
+	call adivkr (spec, high - low, spec, ny)
+	if (clgetb ("cennorm")) {
+	    sum = 0.
+	    nsum = 0
+	    do i = 1, nlines {
+		iy = i - ys + 1
+		if (iy < 1 || iy > ny)
+		    next
+		s = cen + ap_cveval (cv, real (i))
+		ix1 = max (1, int (s))
+		ix2 = min (ncols, int (s + 1))
+		if (ix1 > ix2)
+		    next
+		datain = dbuf + (i - l1) * nc + ix1 - c1
+		if (ix1 == ix2)
+		    sum = sum + Memr[datain]
+		else
+		    sum = sum + (ix2-s)*Memr[datain] + (s-ix1)*Memr[datain+1]
+		nsum = nsum + 1
+	    }
+	    if (nsum > 0) {
+	        sum = (asumr (spec, ny) / ny) / (sum / nsum / gain)
+		call adivkr (spec, sum, spec, ny)
+	    }
+	}
+	if (!IS_INDEF (threshold))
+	    call arltr (spec, ny, threshold, threshold)
+
+	do i = 1, nlines {
+	    if (init == YES) {
+	        dataout = impl2r (out, i)
+	        call amovkr (1., Memr[dataout], ncols)
+	    }
+
+	    iy = i - ys + 1
+	    if (iy < 1 || iy > ny)
+		next
+	    s = ap_cveval (cv, real (i))
+	    x1 = max (0.5, low + s) 
+	    x2 = min (ncols + 0.49, high + s) 
+	    if (x1 > x2)
+		next
+
+	    ix1 = nint (x1)
+	    ix2 = nint (x2)
+
+	    datain = dbuf + (i - l1) * nc + ix1 - c1
+	    if (init == YES)
+		dataout = dataout + ix1 - 1
+	    else
+	        dataout = imps2r (out, ix1, ix2, i, i)
+	    call adivkr (Memr[datain], spec[iy] * gain, Memr[dataout],
+		ix2-ix1+1)
+	}
+
+	call imaddr (out, "CCDMEAN", 1.)
+end
+
+
+# AP_CNORM -- Normalize the input column apertures by the norm spectra.
+
+procedure ap_cnorm (ap, out, gain, dbuf, nc, nl, c1, l1, spec, ny, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+int	ny			# Size of profile array
+int	ys			# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+bool	clgetb()		# Center normalize?
+real	threshold, clgetr()	# Division threshold
+
+int	ncols, nlines, ix, iy, ix1, ix2, iy1, iy2, nsum
+real	cen, low, high, s, sum, ap_cveval(), asumr()
+pointer	sp, y1, y2, cv, datain, dataout, buf, imps2r(), impl2r()
+
+begin
+	threshold = clgetr ("threshold")
+
+	call smark (sp)
+	call salloc (y1, 2 * ny, TY_INT)
+	y1 = y1 - ys
+	y2 = y1 + ny
+
+	cen = AP_CEN(ap,2)
+	low = AP_CEN(ap,2) + AP_LOW(ap,2)
+	high = AP_CEN(ap,2) + AP_HIGH(ap,2)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	# Normalize by the aperture width and apply threshold.
+	call adivkr (spec, high - low, spec, ny)
+	if (clgetb ("cennorm")) {
+	    sum = 0.
+	    nsum = 0
+	    do ix = ys, ys+ny-1 {
+		s = cen + ap_cveval (cv, real (ix))
+		iy1 = max (1, int (s))
+		iy2 = min (nlines, int (s + 1))
+		if (iy1 > iy2)
+		    next
+		datain = dbuf + (ix - l1) * nc + iy1 - c1
+		if (iy1 == iy2)
+		    sum = sum + Memr[datain]
+		else
+		    sum = sum + (iy2-s)*Memr[datain] + (s-iy1)*Memr[datain+1]
+		nsum = nsum + 1
+	    }
+	    if (nsum > 0) {
+	        sum = (asumr (spec, ny) / ny) / (sum / nsum / gain)
+		call adivkr (spec, sum, spec, ny)
+	    }
+	}
+	if (!IS_INDEF (threshold))
+	    call arltr (spec, ny, threshold, threshold)
+
+	do ix = ys, ys+ny-1 {
+	    s = ap_cveval (cv, real (ix))
+	    Memi[y1+ix] = nint (low + s) 
+	    Memi[y2+ix] = nint (high + s) 
+	}
+	call alimi (Memi[y1+ys], 2 * ny, iy1, iy2)
+
+	do iy = 1, nlines {
+	    if (init == YES) {
+	        buf = impl2r (out, iy)
+	        call amovkr (1., Memr[buf], ncols)
+	    }
+
+	    if (iy < iy1 || iy > iy2)
+		next
+
+	    for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) {
+		if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1])
+		    next
+		for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1)
+		    if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2])
+			break
+		ix2 = ix2 - 1
+
+	        datain = dbuf + (ix1 - l1) * nc + iy - c1
+	        if (init == YES)
+		    dataout = buf + ix1 - 1
+	        else
+	            dataout = imps2r (out, ix1, ix2, iy, iy)
+		do ix = ix1, ix2 {
+		    Memr[dataout] = Memr[datain] / spec[ix-ys+1] / gain
+		    datain = datain +  nc
+		    dataout = dataout + 1
+		}
+		ix1 = ix2
+	    }
+	}
+
+	call imaddr (out, "CCDMEAN", 1.)
+
+	call sfree (sp)
+end
+
+
+# AP_LFLAT -- Flatten the input line apertures by the norm spectra.
+
+procedure ap_lflat (ap, out, dbuf, nc, nl, c1, l1, spec, sbuf, profile, nx, ny,
+	xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+pointer	sbuf			# Sky buffer
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+real	threshold, clgetr()	# Division threshold
+
+int	i, ncols, nlines, ix, iy, ix1, ix2
+real	low, high, s, x1, x2, ap_cveval()
+pointer	cv, datain, dataout, sky, imps2r(), impl2r()
+
+begin
+	threshold = clgetr ("threshold")
+	if (IS_INDEF(threshold))
+	    threshold = 0.
+	threshold = max (0., threshold)
+
+	low = AP_CEN(ap,1) + AP_LOW(ap,1)
+	high = AP_CEN(ap,1) + AP_HIGH(ap,1)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do i = 1, nlines {
+	    if (init == YES) {
+	        dataout = impl2r (out, i)
+	        call amovkr (1., Memr[dataout], ncols)
+	    }
+
+	    iy = i - ys + 1
+	    if (iy < 1 || iy > ny)
+		next
+	    s = ap_cveval (cv, real (i))
+	    x1 = max (0.5, low + s) 
+	    x2 = min (ncols + 0.49, high + s) 
+	    if (x1 > x2)
+		next
+
+	    ix1 = nint (x1)
+	    ix2 = nint (x2)
+
+	    datain = dbuf + (i - l1) * nc + ix1 - c1
+	    if (init == YES)
+		dataout = dataout + ix1 - 1
+	    else
+	        dataout = imps2r (out, ix1, ix2, i, i)
+	    if (sbuf != NULL)
+		sky = sbuf + (iy - 1) * nx - xs[iy]
+	    do ix = ix1, ix2 {
+		s = spec[iy] * profile[iy, ix-xs[iy]+1]
+		if (sbuf != NULL)
+		    s = s + Memr[sky+ix]
+		if (s > threshold)
+		    Memr[dataout] = Memr[datain] / s
+		else
+		    Memr[dataout] = 1.
+		datain = datain + 1
+		dataout = dataout + 1
+	    }
+	}
+
+	call imaddr (out, "CCDMEAN", 1.)
+end
+
+
+# AP_CFLAT -- Flatten the input column apertures by the norm spectra.
+
+procedure ap_cflat (ap, out, dbuf, nc, nl, c1, l1, spec, sbuf, profile, nx, ny,
+    xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+pointer	sbuf			# Sky buffer
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+real	threshold, clgetr()	# Division threshold
+
+int	ncols, nlines, ix, iy, ix1, ix2, iy1, iy2
+real	low, high, s, ap_cveval()
+pointer	sp, y1, y2, cv, datain, dataout, sky, buf, imps2r(), impl2r()
+
+begin
+	threshold = clgetr ("threshold")
+	if (IS_INDEF(threshold))
+	    threshold = 0.
+	threshold = max (0., threshold)
+
+	call smark (sp)
+	call salloc (y1, 2 * ny, TY_INT)
+	y1 = y1 - ys
+	y2 = y1 + ny
+
+	low = AP_CEN(ap,2) + AP_LOW(ap,2)
+	high = AP_CEN(ap,2) + AP_HIGH(ap,2)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do ix = ys, ys+ny-1 {
+	    s = ap_cveval (cv, real (ix))
+	    Memi[y1+ix] = nint (low + s) 
+	    Memi[y2+ix] = nint (high + s) 
+	}
+	call alimi (Memi[y1+ys], 2 * ny, iy1, iy2)
+
+	do iy = 1, nlines {
+	    if (init == YES) {
+	        buf = impl2r (out, iy)
+	        call amovkr (1., Memr[buf], ncols)
+	    }
+
+	    if (iy < iy1 || iy > iy2)
+		next
+
+	    for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) {
+		if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1])
+		    next
+		for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1)
+		    if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2])
+			break
+		ix2 = ix2 - 1
+
+	        datain = dbuf + (ix1 - l1) * nc + iy - c1
+	        if (init == YES)
+		    dataout = buf + ix1 - 1
+	        else
+	            dataout = imps2r (out, ix1, ix2, iy, iy)
+	        if (sbuf != NULL)
+		    sky = sbuf - ys * nx + iy - xs[iy]
+	        do ix = ix1, ix2 {
+		    s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1]
+		    if (sbuf != NULL)
+			s = s + Memr[sky+ix*nx]
+		    if (s > threshold)
+		        Memr[dataout] = Memr[datain] / s
+		    else
+		        Memr[dataout] = 1.
+		    datain = datain + nc
+		    dataout = dataout + 1
+	        }
+		ix1 = ix2
+	    }
+	}
+
+	call imaddr (out, "CCDMEAN", 1.)
+
+	call sfree (sp)
+end
+
+
+# AP_LDIFF -- Model residuals.
+
+procedure ap_ldiff (ap, out, gain, dbuf, nc, nl, c1, l1, spec, profile, nx, ny,
+	xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+int	i, ncols, nlines, ix, iy, ix1, ix2
+real	low, high, s, x1, x2, ap_cveval()
+pointer	cv, datain, dataout, imps2r(), impl2r()
+
+begin
+	low = AP_CEN(ap,1) + AP_LOW(ap,1)
+	high = AP_CEN(ap,1) + AP_HIGH(ap,1)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do i = 1, nlines {
+	    if (init == YES) {
+	        dataout = impl2r (out, i)
+	        call aclrr (Memr[dataout], ncols)
+	    }
+
+	    iy = i - ys + 1
+	    if (iy < 1 || iy > ny)
+		next
+	    s = ap_cveval (cv, real (i))
+	    x1 = max (0.5, low + s) 
+	    x2 = min (ncols + 0.49, high + s) 
+	    if (x1 > x2)
+		next
+
+	    ix1 = nint (x1)
+	    ix2 = nint (x2)
+
+	    datain = dbuf + (i - l1) * nc + ix1 - c1
+	    if (init == YES)
+		dataout = dataout + ix1 - 1
+	    else
+	        dataout = imps2r (out, ix1, ix2, i, i)
+	    do ix = ix1, ix2 {
+		s = spec[iy] * profile[iy, ix-xs[iy]+1]
+		Memr[dataout] = (Memr[datain] - s) / gain
+		datain = datain + 1
+		dataout = dataout + 1
+	    }
+	}
+end
+
+
+# AP_CDIFF -- Model residuals
+
+procedure ap_cdiff (ap, out, gain, dbuf, nc, nl, c1, l1, spec, profile, nx, ny,
+    xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+pointer	dbuf			# Data buffer
+int	nc, nl			# Size of data buffer
+int	c1, l1			# Origin of data buffer
+real	spec[ny]		# Normalization spectrum
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+int	ncols, nlines, ix, iy, ix1, ix2, iy1, iy2
+real	low, high, s, ap_cveval()
+pointer	sp, y1, y2, cv, datain, dataout, buf, imps2r(), impl2r()
+
+begin
+	call smark (sp)
+	call salloc (y1, 2 * ny, TY_INT)
+	y1 = y1 - ys
+	y2 = y1 + ny
+
+	low = AP_CEN(ap,2) + AP_LOW(ap,2)
+	high = AP_CEN(ap,2) + AP_HIGH(ap,2)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do ix = ys, ys+ny-1 {
+	    s = ap_cveval (cv, real (ix))
+	    Memi[y1+ix] = nint (low + s) 
+	    Memi[y2+ix] = nint (high + s) 
+	}
+	call alimi (Memi[y1+ys], 2 * ny, iy1, iy2)
+
+	do iy = 1, nlines {
+	    if (init == YES) {
+	        buf = impl2r (out, iy)
+	        call aclrr (Memr[buf], ncols)
+	    }
+
+	    if (iy < iy1 || iy > iy2)
+		next
+
+	    for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) {
+		if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1])
+		    next
+		for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1)
+		    if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2])
+			break
+		ix2 = ix2 - 1
+
+	        datain = dbuf + (ix1 - l1) * nc + iy - c1
+	        if (init == YES)
+		    dataout = buf + ix1 - 1
+	        else
+	            dataout = imps2r (out, ix1, ix2, iy, iy)
+	        do ix = ix1, ix2 {
+		    s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1]
+		    Memr[dataout] = (Memr[datain] - s) / gain
+		    datain = datain + nc
+		    dataout = dataout + 1
+	        }
+		ix1 = ix2
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_LFIT -- Model fit
+
+procedure ap_lfit (ap, out, gain, spec, profile, nx, ny, xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+real	spec[ny]		# Normalization spectrum
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+int	i, ncols, nlines, ix, iy, ix1, ix2
+real	low, high, s, x1, x2, ap_cveval()
+pointer	cv, dataout, imps2r(), impl2r()
+
+begin
+	low = AP_CEN(ap,1) + AP_LOW(ap,1)
+	high = AP_CEN(ap,1) + AP_HIGH(ap,1)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do i = 1, nlines {
+	    if (init == YES) {
+	        dataout = impl2r (out, i)
+	        call aclrr (Memr[dataout], ncols)
+	    }
+
+	    iy = i - ys + 1
+	    if (iy < 1 || iy > ny)
+		next
+	    s = ap_cveval (cv, real (i))
+	    x1 = max (0.5, low + s) 
+	    x2 = min (ncols + 0.49, high + s) 
+	    if (x1 > x2)
+		next
+
+	    ix1 = nint (x1)
+	    ix2 = nint (x2)
+
+	    if (init == YES)
+		dataout = dataout + ix1 - 1
+	    else
+	        dataout = imps2r (out, ix1, ix2, i, i)
+	    do ix = ix1, ix2 {
+		s = spec[iy] * profile[iy, ix-xs[iy]+1]
+		Memr[dataout] = s / gain
+		dataout = dataout + 1
+	    }
+	}
+end
+
+
+# AP_CFIT -- Model fit
+
+procedure ap_cfit (ap, out, gain, spec, profile, nx, ny, xs, ys, init)
+
+pointer	ap			# Aperture structure
+pointer	out			# Output IMIO pointer
+real	gain			# Gain
+real	spec[ny]		# Normalization spectrum
+real	profile[ny,nx]		# Profile
+int	nx, ny			# Size of profile array
+int	xs[ny], ys		# Start of spectrum in image
+int	init			# Fill between apertures with 1?
+
+int	ncols, nlines, ix, iy, ix1, ix2, iy1, iy2
+real	low, high, s, ap_cveval()
+pointer	sp, y1, y2, cv, dataout, buf, imps2r(), impl2r()
+
+begin
+	call smark (sp)
+	call salloc (y1, 2 * ny, TY_INT)
+	y1 = y1 - ys
+	y2 = y1 + ny
+
+	low = AP_CEN(ap,2) + AP_LOW(ap,2)
+	high = AP_CEN(ap,2) + AP_HIGH(ap,2)
+	cv = AP_CV(ap)
+	ncols = IM_LEN(out, 1)
+	nlines = IM_LEN(out, 2)
+
+	do ix = ys, ys+ny-1 {
+	    s = ap_cveval (cv, real (ix))
+	    Memi[y1+ix] = nint (low + s) 
+	    Memi[y2+ix] = nint (high + s) 
+	}
+	call alimi (Memi[y1+ys], 2 * ny, iy1, iy2)
+
+	do iy = 1, nlines {
+	    if (init == YES) {
+	        buf = impl2r (out, iy)
+	        call aclrr (Memr[buf], ncols)
+	    }
+
+	    if (iy < iy1 || iy > iy2)
+		next
+
+	    for (ix1=ys; ix1<=ys+ny-1; ix1=ix1+1) {
+		if (iy < Memi[y1+ix1] || iy > Memi[y2+ix1])
+		    next
+		for (ix2=ix1+1; ix2<=ys+ny-1; ix2=ix2+1)
+		    if (iy < Memi[y1+ix2] || iy > Memi[y2+ix2])
+			break
+		ix2 = ix2 - 1
+
+	        if (init == YES)
+		    dataout = buf + ix1 - 1
+	        else
+	            dataout = imps2r (out, ix1, ix2, iy, iy)
+	        do ix = ix1, ix2 {
+		    s = spec[ix-ys+1] * profile[ix-ys+1, iy-xs[ix-ys+1]+1]
+		    Memr[dataout] = s / gain
+		    dataout = dataout + 1
+	        }
+		ix1 = ix2
+	    }
+	}
+
+	call sfree (sp)
+end