Initial commit

1 files changed, 254 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/fitsky/apskybuf.x b/noao/digiphot/apphot/fitsky/apskybuf.x
new file mode 100644
index 00000000..cacd4737
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apskybuf.x
@@ -0,0 +1,254 @@
+include <imhdr.h>
+include <math.h>
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APSKYBUF -- Procedure to fetch the sky pixels given the pointer to the
+# IRAF image, the coordinates of the center and the size of the apphot
+# sky annulus.
+
+int procedure apskybuf (ap, im, wx, wy)
+
+pointer	ap		# pointer to apphot structure
+pointer	im		# pointer to the IRAF image
+real	wx, wy		# center coordinates
+
+int	lenbuf
+pointer	sky
+real	annulus, dannulus, datamin, datamax
+int	ap_skypix(), ap_bskypix()
+
+begin
+	# Check for 0 radius annulus.
+	sky = AP_PSKY(ap)
+	annulus = AP_ANNULUS(sky) * AP_SCALE(ap)
+	dannulus = AP_DANNULUS(sky) * AP_SCALE(ap)
+	if (dannulus <= 0.0)
+	    return (AP_NOSKYAREA)
+
+	# Allocate space for sky pixels.
+	lenbuf = PI * (2.0 * annulus + dannulus + 1.0) * (dannulus + 0.5)
+
+	if (lenbuf != AP_LENSKYBUF(sky)) {
+	    if (AP_SKYPIX(sky) != NULL)
+		call mfree (AP_SKYPIX(sky), TY_REAL)
+	    call malloc (AP_SKYPIX(sky), lenbuf, TY_REAL)
+	    if (AP_COORDS(sky) != NULL)
+		call mfree (AP_COORDS(sky), TY_INT)
+	    call malloc (AP_COORDS(sky), lenbuf, TY_INT)
+	    if (AP_INDEX(sky) != NULL)
+		call mfree (AP_INDEX(sky), TY_INT)
+	    call malloc (AP_INDEX(sky), lenbuf, TY_INT)
+	    if (AP_SWGT(sky) != NULL)
+		call mfree (AP_SWGT(sky), TY_REAL)
+	    call malloc (AP_SWGT(sky), lenbuf, TY_REAL)
+	    AP_LENSKYBUF(sky) = lenbuf
+	}
+
+	# Fetch the sky pixels.
+	if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap))) {
+	    AP_NSKYPIX(sky) = ap_skypix (im, wx, wy, annulus, (annulus +
+	        dannulus), Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+		AP_SXC(sky), AP_SYC(sky), AP_SNX(sky), AP_SNY(sky))
+	    AP_NBADSKYPIX(sky) = 0
+	} else {
+	    if (IS_INDEFR(AP_DATAMIN(ap)))
+		datamin = -MAX_REAL
+	    else
+		datamin = AP_DATAMIN(ap)
+	    if (IS_INDEFR(AP_DATAMAX(ap)))
+		datamax = MAX_REAL
+	    else
+		datamax = AP_DATAMAX(ap)
+	    AP_NSKYPIX(sky) = ap_bskypix (im, wx, wy, annulus, (annulus +
+	        dannulus), datamin, datamax, Memr[AP_SKYPIX(sky)],
+		Memi[AP_COORDS(sky)], AP_SXC(sky), AP_SYC(sky), AP_SNX(sky),
+		AP_SNY(sky), AP_NBADSKYPIX(sky))
+	}
+
+	if (AP_NSKYPIX(sky) <= 0) {
+	    if (AP_NBADSKYPIX(sky) <= 0)
+	        return (AP_SKY_OUTOFBOUNDS)
+	    else
+	        return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
+
+
+# AP_SKYPIX -- Procedure to fetch the sky pixels from the image
+
+int procedure ap_skypix (im, wx, wy, rin, rout, skypix, coords, xc, yc,
+        nx, ny)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of sky annulus
+real	rin, rout	# inner and outer radius of sky annulus
+real	skypix[ARB]	# skypixels
+int	coords[ARB]	# sky subraster coordinates [i + nx * (j - 1)]
+real	xc, yc		# center of sky subraster
+int	nx, ny		# max dimensions of sky subraster (output)
+
+int	i, j, ncols, nlines, c1, c2, l1, l2, nskypix
+pointer	buf
+real	xc1, xc2, xl1, xl2, rin2, rout2, rj2, r2
+pointer	imgs2r()
+
+#pointer	tbuf
+
+begin
+	if (rout <= rin)
+	    return (0)
+
+	# Test for out of bounds sky regions.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - rout
+	xc2 = wx + rout
+	xl1 = wy - rout
+	xl2 = wy + rout
+	if (xc2 < 1.0 || xc1 > real (ncols) || xl2 < 1.0 || xl1 > real (nlines))
+	    return (0)
+
+	# Compute the column and line limits.
+	c1 = max (1.0, min (real (ncols), wx - rout)) + 0.5
+	c2 = min (real (ncols), max (1.0, wx + rout)) + 0.5
+	l1 = max (1.0, min (real (nlines), wy - rout)) + 0.5
+	l2 = min (real (nlines), max (1.0, wy + rout)) + 0.5
+	nx = c2 - c1 + 1
+	ny = l2 - l1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+
+	# Fetch the sky pixels.
+	rin2 = rin ** 2
+	rout2 = rout ** 2
+	nskypix = 0
+
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    rj2 = (wy - j) ** 2
+	    do i = c1, c2 {
+	        r2 = (wx - i) ** 2 + rj2
+		if (r2 > rin2 && r2 <= rout2) {
+		    skypix[nskypix+1] = Memr[buf+i-c1]
+		    coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		    nskypix = nskypix + 1
+		}
+	    }
+	}
+
+	#buf = imgs2r (im, c1, c2, l1, l2)
+	#tbuf = buf
+	#do j = l1,  l2 {
+	    #rj2 = (wy - j) ** 2
+	    #do i = c1, c2 {
+	        #r2 = (wx - i) ** 2 + rj2
+		#if (r2 > rin2 && r2 <= rout2) {
+		    #skypix[nskypix+1] = Memr[tbuf+i-c1]
+		    #coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		    #nskypix = nskypix + 1
+		#}
+	    #}
+	    #tbuf = tbuf + nx
+	#}
+
+	return (nskypix)
+end
+
+
+# AP_BSKYPIX -- Procedure to fetch the sky pixels from the image
+
+int procedure ap_bskypix (im, wx, wy, rin, rout, datamin, datamax,
+	skypix, coords, xc, yc, nx, ny, nbad)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of sky annulus
+real	rin, rout	# inner and outer radius of sky annulus
+real	datamin		# minimum good value
+real	datamax		# maximum good value
+real	skypix[ARB]	# skypixels
+int	coords[ARB]	# sky subraster coordinates [i + nx * (j - 1)]
+real	xc, yc		# center of sky subraster
+int	nx, ny		# max dimensions of sky subraster (output)
+int	nbad		# number of bad pixels
+
+int	i, j, ncols, nlines, c1, c2, l1, l2, nskypix
+pointer	buf
+real	xc1, xc2, xl1, xl2, rin2, rout2, rj2, r2, pixval
+pointer	imgs2r()
+
+#pointer	tbuf
+
+begin
+	if (rout <= rin)
+	    return (0)
+
+	# Test for out of bounds sky regions.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - rout
+	xc2 = wx + rout
+	xl1 = wy - rout
+	xl2 = wy + rout
+	if (xc2 < 1.0 || xc1 > real (ncols) || xl2 < 1.0 || xl1 > real (nlines))
+	    return (0)
+
+	# Compute the column and line limits.
+	c1 = max (1.0, min (real (ncols), wx - rout)) + 0.5
+	c2 = min (real (ncols), max (1.0, wx + rout)) + 0.5
+	l1 = max (1.0, min (real (nlines), wy - rout)) + 0.5
+	l2 = min (real (nlines), max (1.0, wy + rout)) + 0.5
+	nx = c2 - c1 + 1
+	ny = l2 - l1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+
+	rin2 = rin ** 2
+	rout2 = rout ** 2
+	nskypix = 0
+	nbad = 0
+
+	# Fetch the sky pixels.
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    rj2 = (wy - j) ** 2
+	    do i = c1, c2 {
+	        r2 = (wx - i) ** 2 + rj2
+		if (r2 > rin2 && r2 <= rout2) {
+		    pixval = Memr[buf+i-c1] 
+		    if (pixval < datamin || pixval > datamax)
+		        nbad = nbad + 1
+		    else {
+		        skypix[nskypix+1] = pixval
+		        coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		        nskypix = nskypix + 1
+		    }
+		}
+	    }
+	}
+
+	#buf = imgs2r (im, c1, c2, l1, l2)
+	#tbuf = buf
+	#do j = l1, l2 {
+	    #rj2 = (wy - j) ** 2
+	    #do i = c1, c2 {
+	        #r2 = (wx - i) ** 2 + rj2
+		#if (r2 > rin2 && r2 <= rout2) {
+		    #pixval = Memr[tbuf+i-c1] 
+		    #if (pixval < datamin || pixval > datamax)
+		        #nbad = nbad + 1
+		    #else {
+		        #skypix[nskypix+1] = pixval
+		        #coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		        #nskypix = nskypix + 1
+		    #}
+		#}
+	    #}
+	    #tbuf = tbuf + nx
+	#}
+
+	return (nskypix)
+end