Repatch (from linux) of OSX IRAF

1 files changed, 172 insertions, 0 deletions

diff --git a/noao/digiphot/daophot/psf/dppsubrast.x b/noao/digiphot/daophot/psf/dppsubrast.x
new file mode 100644
index 00000000..1fafb408
--- /dev/null
+++ b/noao/digiphot/daophot/psf/dppsubrast.x
@@ -0,0 +1,172 @@
+include <mach.h>
+include <imhdr.h>
+include "../lib/daophotdef.h"
+include "../lib/psfdef.h"
+
+# DP_PSUBRAST -- Fetch the prospective PSF star data and check it for bad
+# pixel values.
+
+pointer procedure dp_psubrast (dao, im, lowbad, highbad, x1, x2, y1, y2,
+	saturated)
+
+pointer	dao		# pointer to the daophot structure
+pointer	im		# pointer to the input image
+real	lowbad, highbad	# minimum and maximum good data values
+int	x1, x2		# output x limits of the extracted subraster
+int	y1, y2		# output y limits of the extracted subraster
+int	saturated	# is the star saturated ?
+
+pointer	psf, buf
+real	psfrad, fitrad
+int	dp_chksr()
+pointer	dp_subrast()
+
+begin
+	# Initialize.
+	psf = DP_PSF(dao)
+	buf = NULL
+	psfrad = DP_PSFRAD(dao)
+	fitrad = DP_FITRAD(dao)
+
+	# Get the data.
+ 	buf = dp_subrast (im, DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), psfrad,
+	    x1, x2, y1, y2)
+	if (buf == NULL)
+	    return (NULL)
+
+	# Check for bad pixels in subraster, compute the min and max.
+	if (dp_chksr (DP_CUR_PSFX(psf), DP_CUR_PSFY(psf), Memr[buf],
+	    x2 - x1 + 1, y2 - y1 + 1, x1, y1, psfrad, fitrad,
+	    lowbad, highbad, saturated, DP_CUR_PSFMIN(psf),
+	        DP_CUR_PSFMAX(psf), DP_CUR_PSFGMAX(psf)) == ERR) {
+	    call mfree (buf, TY_REAL)
+	    return (NULL)
+	}
+	
+	return (buf)
+end
+
+
+# DL_LSUBRAST -- Give a valid PSF star compute the limits of the data to
+# be extracted around it.
+
+int procedure dp_lsubrast (im, xcen, ycen, radius, x1, x2, y1, y2)
+
+pointer	im			# input image descriptor
+real	xcen, ycen		# center of subraster
+real	radius			# radius of the box
+int	x1, y1, x2, y2		# boundaries of subraster
+
+begin
+	# Calculate start position of extraction box.
+	x1 = int (xcen - radius) - 2
+	x2 = int (xcen + radius) + 3
+	y1 = int (ycen - radius) - 2
+	y2 = int (ycen + radius) + 3
+	if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1)
+	    return (ERR)
+
+	x1 = max (1, x1)
+	x2 = min (IM_LEN(im,1), x2)
+	y1 = max (1, y1)
+	y2 = min (IM_LEN(im,2), y2)
+	return (OK)
+end
+
+
+# DP_SUBRAST -- Given a valid PSF star extract the data around it.
+
+pointer procedure dp_subrast (im, xcen, ycen, radius, x1, x2, y1, y2)
+
+pointer	im			# input image descriptor
+real	xcen, ycen		# center of subraster
+real	radius			# radius of the box
+int	x1, y1, x2, y2		# boundaries of subraster
+
+int	j, ncols
+pointer	buf, ptr, imbuf
+pointer imgs2r()
+
+begin
+	# Calculate start position of extraction box.
+	x1 = int (xcen - radius) - 2
+	x2 = int (xcen + radius) + 3
+	y1 = int (ycen - radius) - 2
+	y2 = int (ycen + radius) + 3
+	if (x1 > IM_LEN(im,1) || x2 < 1 || y1 > IM_LEN(im,2) || y2 < 1)
+	    return (NULL)
+
+	x1 = max (1, x1)
+	x2 = min (IM_LEN(im,1), x2)
+	y1 = max (1, y1)
+	y2 = min (IM_LEN(im,2), y2)
+	call malloc (buf, (x2 - x1 + 1) * (y2 - y1 + 1), TY_REAL)
+
+	ptr = buf
+	ncols = x2 - x1 + 1
+	do j = y1, y2 {
+	    imbuf = imgs2r (im, x1, x2, j, j)
+	    call amovr (Memr[imbuf], Memr[ptr], ncols) 
+	    ptr = ptr + ncols
+	}
+
+	return (buf)
+end
+
+
+# DP_CHKSR -- Check the input subraster for bad pixels.
+
+int procedure dp_chksr (x, y, sr, xdim, ydim, x1, y1, psfrad, fitrad, lowbad,
+        highbad, saturated, dmin, dmax, gmax)
+
+real	x, y			# position of the star
+real	sr[xdim,ydim]		# the data subraster
+int	xdim, ydim		# the dimensions of the subraster
+int	x1, y1			# the lower left hand coordinates of the array
+real	psfrad			# the psf radius
+real	fitrad			# the fitting radius
+real	lowbad, highbad		# the good data limits
+int	saturated		# is the star saturated
+real	dmin, dmax		# output data limits
+real	gmax			# maximum good data limit
+
+int	i,j
+real	pradsq, fradsq, dy2, r2
+
+begin
+	pradsq = psfrad * psfrad
+	fradsq = fitrad * fitrad
+	dmin = MAX_REAL
+	dmax = -MAX_REAL
+	gmax = -MAX_REAL
+	saturated = NO
+
+	# Loop through the pixels looking for bad values.
+	do j = 1, ydim {
+	    dy2 = (y - (y1 + j -1)) ** 2 
+	    if (dy2 > pradsq)
+		next
+	    do i = 1, xdim {
+		r2 = (x - (x1 + i - 1)) ** 2 + dy2
+	        if (r2 > pradsq)
+		    next
+		if (sr[i,j] < dmin)
+		    dmin = sr[i,j]
+		if (sr[i,j] > dmax)
+		    dmax = sr[i,j]
+		if (sr[i,j] < lowbad || sr[i,j] > highbad) {
+		    if (r2 <= fradsq) {
+		        if (sr[i,j] < lowbad)
+		    	    return (ERR)
+			else if (saturated == NO)
+			    saturated = YES
+		    }
+		} else {
+		    if (sr[i,j] > gmax)
+		        gmax = sr[i,j]
+		}
+	    }
+	}
+
+	return (OK)
+end