1 files changed, 368 insertions, 0 deletions
diff --git a/noao/twodspec/apextract/apskyeval.x b/noao/twodspec/apextract/apskyeval.x
new file mode 100644
index 00000000..05f47f14
--- /dev/null
+++ b/noao/twodspec/apextract/apskyeval.x
@@ -0,0 +1,368 @@
+include	<math/iminterp.h>
+include <mach.h>
+include	"apertures.h"
+
+# Background fitting types
+define	BACKGROUND	"|none|average|median|minimum|fit|"
+define	B_NONE		1
+define	B_AVERAGE	2
+define	B_MEDIAN	3
+define	B_MINIMUM	4
+define	B_FIT		5
+
+define	NSAMPLE		20	# Maximum number of background sample regions
+
+
+# AP_SKYEVAL -- Evaluate sky within aperture.
+#
+# The sky pixels specified by the background sample string are used to
+# determine a sky function at each line which is then evaluated for each
+# pixel in the aperture as given by the SBUF array with starting offsets
+# given by XS.  The fitting consists of either a straight average or a
+# function fit using ICFIT.  The sky regions are specified relative to the
+# aperture center.  To avoid systematics due to shifting of the aperture
+# relative to the integer pixel positions the sky regions are linearly
+# interpolated.  The average uses the integral of the interpolation
+# function within the sample region endpoints.  The fit samples the
+# interpolation on a pixel grid with the aperture exactly centered on
+# a pixel.  A crude sky variance is computed for each line based solely
+# on the variance model and the square root of the number of "pixels"
+# used for the fit.  This variance is used to boost the variance of
+# the sky subtracted spectrum during variance weighting.  Because sky
+# noise may be significant in short slits a box car smoothing may be
+# used giving a lower variance per pixel but bad errors near sky lines.
+# An unweighted aperture sum of the sky is returned in case the user
+# wants to save the subtracted 1D sky spectrum.
+
+procedure ap_skyeval (im, ap, dbuf, nc, nl, c1, l1, sbuf, svar, sky, nx, ny,
+	xs, ys, nsubaps, rdnoise)
+
+pointer	im		# IMIO pointer
+pointer	ap		# Aperture structure
+pointer	dbuf		# Data buffer
+int	nc, nl		# Size of data buffer
+int	c1, l1		# Origin of data buffer
+real	sbuf[nx,ny]	# Sky values
+real	svar[ny]	# Sky variances
+real	sky[ny,nsubaps]	# Extracted sky (out)
+int	nx, ny		# Size of profile array
+int	xs[ny], ys	# Origin of profile array
+int	nsubaps		# Number of subapertures
+real	rdnoise		# Readout noise in RMS data numbers.
+
+int	bkg		# Background type
+int	skybox		# Sky box car smoothing
+
+int	i, j, ix1, ix2, nsample, nsky, nfit, ix, iy
+real	center, xmin, xmax, a, b, c, s, avg
+pointer	ic, cv, cv1, asi, sp, str, data, as, bs, x, y, w
+
+int	apgwrd(), apgeti(), ctor()
+real	ic_getr(), ap_cveval(), asieval(), asigrl(), amedr()
+errchk	salloc, ic_fit
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get CL parameters and set shift and fitting function pointers.
+	bkg = apgwrd ("background", Memc[str], SZ_LINE, BACKGROUND)
+	skybox = apgeti ("skybox")
+
+	cv = AP_CV(ap)
+	ic = AP_IC(ap)
+
+	# Set center and maximum limits relative to data buffer.
+	# The limits are required to overlap the aperture and include
+	# an extra point at each end for interpolation.  Shifts
+	# and boundary limits will be enforced later.
+
+	i = AP_AXIS(ap)
+	center = AP_CEN(ap,i)
+	xmin = center + min (AP_LOW(ap,i), ic_getr (ic, "xmin"))
+	xmax = center + max (AP_HIGH(ap,i), ic_getr (ic, "xmax"))
+	ix1 = nint (xmin) - 1
+	ix2 = nint (xmax) + 1
+	nfit = ix2 - ix1 + 1
+
+	# Allocate memory and parse sample string.
+	# The colons in the sample string must be changed to avoid
+	# sexigesimal interpretation.
+
+	call salloc (as, NSAMPLE, TY_REAL)
+	call salloc (bs, NSAMPLE, TY_REAL)
+
+	call ic_gstr (ic, "sample", Memc[str], SZ_LINE)
+	for (i=str; Memc[i]!=EOS; i=i+1)
+	    if (Memc[i] == ':')
+		Memc[i] = '$'
+
+	nsample = 0
+	for (i=1; Memc[str+i-1]!=EOS; i=i+1) {
+	    if (ctor (Memc[str], i, a) > 0) {
+		i = i - 1
+		if (Memc[str+i] == '$') {
+		    i = i + 2
+		    if (ctor (Memc[str], i, b) > 0) {
+			i = i - 1
+			Memr[as+nsample] = center + min (a, b)
+			Memr[bs+nsample] = center + max (a, b)
+			nsample = nsample + 1
+			if (nsample == NSAMPLE)
+			    break
+		    }
+	        }
+	    }
+	}
+
+	if (nsample == 0) {
+	    Memr[as] = xmin
+	    Memr[bs] = xmax
+	    nsample = 1
+	}
+
+	if (bkg == B_MEDIAN)
+	    call salloc (y, nfit, TY_REAL)
+	else if (bkg == B_FIT) {
+	    call salloc (x, nfit, TY_REAL)
+	    call salloc (y, nfit, TY_REAL)
+	    call salloc (w, nfit, TY_REAL)
+	}
+
+	# Initialize the image interpolator.
+	call asiinit (asi, II_LINEAR)
+
+	# Determine sky at each dispersion point.
+	call aclrr (svar, ny)
+	do iy = 1, ny {
+
+	    # Fit image interpolation function including extra points
+	    # and apply image boundary limits.
+
+	    i = iy + ys - 1
+	    s = ap_cveval (cv, real (i))
+	    ix1 = max (c1, nint (xmin + s) - 1)
+	    ix2 = min (c1+nc-1, nint (xmax + s) + 1)
+	    nfit = ix2 - ix1 + 1
+	    if (nfit < 3) {
+		call aclrr (sbuf[1,iy], nx)
+		svar[iy] = 0.
+		next
+	    }
+	    data = dbuf + (i - l1) * nc + ix1 - c1
+	    if (bkg == B_AVERAGE || bkg == B_FIT) {
+		iferr (call asifit (asi, Memr[data], nfit)) {
+		    call aclrr (sbuf[1,iy], nx)
+		    svar[iy] = 0.
+		    next
+		}
+	    }
+
+	    # Determine background
+	    switch (bkg) {
+	    case B_AVERAGE:
+		# The background is computed by integrating the interpolator 
+	        avg = 0.
+		nsky = 0
+	        c = 0.
+		for (i=0; i < nsample; i=i+1) {
+		    a = max (real (ix1), Memr[as+i] + s) - ix1 + 1
+		    b = min (real (ix2), Memr[bs+i] + s) - ix1 + 1
+		    if (b - a > 0.) {
+			avg = avg + asigrl (asi, a, b)
+			c = c + b - a
+			nsky = nsky + nint (b) - nint(a) + 1
+		    }
+	        }
+	        if (c > 0.)
+	            avg = avg / c
+		call amovkr (avg, sbuf[1,iy], nx) 
+		if (nsky > 1)
+		    svar[iy] = max (0., (rdnoise + avg) / (nsky - 1))
+	    case B_MEDIAN:
+		# The background is computed by the median pixel
+	        avg = 0.
+		nsky = 0
+		for (i=0; i < nsample; i=i+1) {
+		    a = max (real (ix1), Memr[as+i] + s) - ix1 + 1
+		    b = min (real (ix2), Memr[bs+i] + s) - ix1 + 1
+		    do j = nint (a), nint (b) {
+			Memr[y+nsky] = Memr[data+j-1]
+			nsky = nsky + 1
+		    }
+	        }
+		if (nsky > 0)
+		    avg = amedr (Memr[y], nsky)
+		call amovkr (avg, sbuf[1,iy], nx) 
+		if (nsky > 1)
+		    svar[iy] = max (0., (rdnoise + avg) / (nsky - 1))
+	    case B_MINIMUM:
+		# The background is computed by the minimum pixel
+	        avg = MAX_REAL
+		nsky = 0
+		for (i=0; i < nsample; i=i+1) {
+		    a = max (real (ix1), Memr[as+i] + s) - ix1 + 1
+		    b = min (real (ix2), Memr[bs+i] + s) - ix1 + 1
+		    do j = nint (a), nint (b) {
+			avg = min (avg, Memr[data+j-1])
+			nsky = nsky + 1
+		    }
+	        }
+		if (nsky == 0)
+		    avg = 0
+		call amovkr (avg, sbuf[1,iy], nx) 
+		if (nsky > 1)
+		    svar[iy] = max (0., (rdnoise + avg) / (nsky - 1))
+	    case B_FIT:
+		# The fitting is done in a coordinate system relative to
+		# aperture center.
+
+		c = center + s
+	        a = ix1 + c - int (c)
+	        do i = 1, nfit-1 {
+	            Memr[x+i-1] = nint (1000. * (a - c)) / 1000.
+		    Memr[y+i-1] = asieval (asi, a-ix1+1)
+		    Memr[w+i-1] = 1.
+		    a = a + 1.
+	        }
+
+	        iferr {
+	            call ic_fit (ic, cv1, Memr[x], Memr[y], Memr[w], nfit-1,
+	                YES, YES, YES, YES)
+
+		    avg = 0.
+	            do i = 1, nx {
+		        a = xs[iy] + i - 1
+			b = ap_cveval (cv1, a - c)
+			avg = avg + b
+			sbuf[i,iy] = b
+	            }
+		    avg = avg / nx
+	        } then {
+		    avg = 0.
+		    call aclrr (sbuf[1,iy], nx)
+		}
+
+		nsky = 0.
+		for (i=0; i < nsample; i=i+1) {
+		    a = max (real (ix1), Memr[as+i] + s) - ix1 + 1
+		    b = min (real (ix2), Memr[bs+i] + s) - ix1 + 1
+		    nsky = nsky + nint (b) - nint (a) + 1
+		}
+		if (nsky > 1)
+		    svar[iy] = max (0., (rdnoise + avg) / (nsky - 1))
+	    }
+	}
+
+	# Do box car smoothing if desired.
+	if (skybox > 1) {
+	    ix2 = skybox ** 2
+	    avg = 0.
+	    a = 0.
+	    iy = 1
+	    for (i=1; i<=skybox; i=i+1) {
+		avg = avg + sbuf[1,i]
+		a = a + svar[i]
+	    }
+	    for (; i<=ny; i=i+1) {
+		b = sbuf[1,iy]
+		c = svar[iy]
+		sbuf[1,iy] = avg / skybox
+		svar[iy] = a / ix2
+		avg = avg + sbuf[1,i] - b
+		a = a + svar[i] - c
+		iy = iy + 1
+	    }
+	    sbuf[1,iy] = avg / skybox
+	    svar[iy] = a / ix2
+	    i = ny - skybox + 1
+	    for (iy=ny; iy > ny-skybox/2; iy=iy-1)
+		svar[iy] = svar[i]
+	    for (; i > 1; i=i-1) {
+		svar[iy] = svar[i]
+		iy = iy - 1
+	    }
+	    for (; iy > 1; iy=iy-1)
+		svar[iy] = svar[1]
+
+	    switch (bkg) {
+	    case B_AVERAGE, B_MEDIAN, B_MINIMUM:
+		i = ny - skybox + 1
+	        for (iy=ny; iy > ny-skybox/2; iy=iy-1)
+		    call amovkr (sbuf[1,i], sbuf[1,iy], nx)
+	        for (; i > 1; i=i-1) {
+		    call amovkr (sbuf[1,i], sbuf[1,iy], nx)
+		    iy = iy - 1
+	        }
+	        for (; iy > 1; iy=iy-1)
+		    call amovkr (sbuf[1,1], sbuf[1,iy], nx)
+	    case B_FIT:
+	        i = ny - skybox + 1
+	        for (iy=ny; iy > ny-skybox/2; iy=iy-1)
+		    sbuf[1,iy] = sbuf[1,i]
+	        for (; i > 1; i=i-1) {
+		    sbuf[1,iy] = sbuf[1,i]
+		    iy = iy - 1
+	        }
+	        for (; iy > 1; iy=iy-1)
+		    sbuf[1,iy] = sbuf[1,1]
+		do ix1 = 2, nx {
+		    avg = 0.
+		    iy = 1
+		    for (i=1; i<=skybox; i=i+1)
+			avg = avg + sbuf[ix1,i]
+		    for (; i<=ny; i=i+1) {
+			b = sbuf[ix1,iy]
+			sbuf[ix1,iy] = avg / skybox
+			avg = avg + sbuf[ix1,i] - b
+			iy = iy + 1
+		    }
+		    sbuf[ix1,iy] = avg / skybox
+		    i = ny - skybox + 1
+		    for (iy=ny; iy > ny-skybox/2; iy=iy-1)
+			sbuf[ix1,iy] = sbuf[ix1,i]
+		    for (; i > 1; i=i-1) {
+			sbuf[ix1,iy] = sbuf[ix1,i]
+			iy = iy - 1
+		    }
+		    for (; iy > 1; iy=iy-1)
+			sbuf[ix1,iy] = sbuf[ix1,1]
+		}
+	    }
+	}
+
+	# Compute the unweighted aperture sky spectrum.
+	i = AP_AXIS(ap)
+	a = AP_CEN(ap,i) + AP_LOW(ap,i)
+	b = AP_CEN(ap,i) + AP_HIGH(ap,i)
+	c = (b - a) / nsubaps
+
+	do iy = 1, ny {
+	    data = dbuf + (iy + ys - 1 - l1) * nc + xs[iy] - c1 - 1
+	    s = ap_cveval (cv, real (iy + ys - 1)) - c1 + 1
+	    do i = 1, nsubaps {
+		xmin = max (0.5, a + (i - 1) * c + s) + c1 - xs[iy]
+		xmax = min (nc + 0.49, a + i * c + s) + c1 - xs[iy]
+		if (xmin >= xmax) {
+		    sky[iy,i] = 0.
+		    next
+		}
+		ix1 = nint (xmin)
+		ix2 = nint (xmax)
+
+		if (ix1 == ix2)
+		    sky[iy,i] = (xmax - xmin) * sbuf[ix1,iy]
+		else {
+		    sky[iy,i] = (ix1 - xmin + 0.5) * sbuf[ix1,iy]
+		    sky[iy,i] = sky[iy,i] + (xmax - ix2 + 0.5) * sbuf[ix2,iy]
+		}
+		do ix = ix1+1, ix2-1
+		    sky[iy,i] = sky[iy,i] + sbuf[ix,iy]
+	    }
+	}
+
+	if (bkg == B_FIT)
+	    call cvfree (cv1)
+	call asifree (asi)
+	call sfree (sp)
+end