Repatch (from linux) of OSX IRAF

1 files changed, 304 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/center/apcplot.x b/noao/digiphot/apphot/center/apcplot.x
new file mode 100644
index 00000000..6e2845c3
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcplot.x
@@ -0,0 +1,304 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+# AP_CPLOT -- Procedure to compute radial profile plots for the centering
+# routine.
+
+procedure ap_cplot (ap, sid, gd, makeplot)
+
+pointer	ap		# the pointer to the apphot structure
+int	sid		# the output sequence number of the star
+pointer	gd		# the graphics stream
+int	makeplot	# the make a plot ?
+
+int	nx, ny
+pointer	ctr, sp, str, r, gt
+real	xcenter, ycenter, xc, yc, rmin, rmax, imin, imax
+real	u1, u2, v1, v2, x1, x2, y1, y2
+int	apstati()
+pointer	ap_gtinit()
+real	apstatr()
+
+begin
+	# Check for enabled graphics stream.
+	if (gd == NULL || makeplot == NO)
+	    return
+
+	# Check for defined center.
+	xcenter = apstatr (ap, XCENTER)
+	ycenter = apstatr (ap, YCENTER)
+	if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter))
+	    return
+
+	# Check for a centering algorithm.
+	if (apstati (ap, CENTERFUNCTION) == AP_NONE)
+	    return
+
+	# Get the pixel buffer parameters.
+	ctr = AP_PCENTER(ap)
+	nx = AP_CNX(ctr)
+	ny = AP_CNY(ctr)
+	if (AP_CTRPIX(ctr) == NULL || nx <= 0 || ny <= 0)
+	    return
+	xc = AP_CXC(ctr) + (xcenter - apstatr (ap, CXCUR))
+	yc = AP_CYC(ctr) + (ycenter - apstatr (ap, CYCUR))
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (r, nx * ny, TY_REAL)
+
+	# Compute the radii and the plot limits.
+	call ap_ijtor2 (Memr[r], nx, ny, xc, yc)
+	call alimr (Memr[r], nx * ny, rmin, rmax)
+	call alimr (Memr[AP_CTRPIX(ctr)], nx * ny, imin, imax)
+
+	# Reactivate the work station.
+	call greactivate (gd, 0)
+
+	# Save old viewport and window coordinates.
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Initialize the plot.
+	call apstats (ap, IMROOT, Memc[str], SZ_LINE)
+	call sprintf (Memc[str], SZ_LINE, "%s  Star %d")
+	    call pargstr (Memc[str])
+	    call pargi (sid)
+	gt = ap_gtinit (Memc[str], apstatr (ap,OXINIT), apstatr (ap, OYINIT))
+
+	# Make the plot.
+	call gclear (gd)
+	call ap_cpset (gd, gt, ap, rmin, rmax, imin, imax)
+	if (apstati (ap, POSITIVE) == YES)
+	    call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny,
+	        rmin, rmax, imin + EPSILONR, imax, "plus")
+	else
+	    call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny,
+	        rmin, rmax, imin, imax - EPSILONR, "plus")
+	call ap_cpreset (gd, gt, ap, rmin, rmax, imin, imax)
+	call ap_cpannotate (gd, ap)
+
+	# Restore the viewport and window coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	# Free the space.
+	call ap_gtfree (gt)
+	call gdeactivate (gd, 0)
+	call sfree (sp)
+end
+
+
+# AP_CPSET -- Procedure to set up the parameters for the center radial profile
+# plot.
+
+procedure ap_cpset (gd, gt, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# the graphics stream
+pointer	gt		# the gtools pointer
+pointer	ap		# the apphot pointer
+real	xmin, xmax	# the minimum and maximum radial distance
+real	ymin, ymax	# the minimum and maximum of y axis (ymin not used)
+
+int	fd
+pointer	sp, str, title
+real	scale, aspect, datalimit, skysigma, threshold, vx1, vx2, vy1, vy2
+real	cthreshold
+int	stropen(), apstati()
+real	apstatr(), gstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	# Encode the parameter string.
+	fd = stropen (Memc[str], SZ_LINE, WRITE_ONLY)
+	call sysid (Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+	call fprintf (fd,
+	    "Center: xc=%0.2f yc=%0.2f xerr=%0.2f yerr=%0.2f\n")
+	    call pargr (apstatr (ap, OXCENTER))
+	    call pargr (apstatr (ap, OYCENTER))
+	    call pargr (apstatr (ap, XERR))
+	    call pargr (apstatr (ap, YERR))
+	call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+	call strclose (fd)
+
+	# Store some default plotting parameters.
+	scale = apstatr (ap, SCALE)
+	aspect = gstatr (gd, G_ASPECT)
+	call gsetr (gd, G_ASPECT, 0.75)
+
+	# Set the labels and window.
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gseti (gd, G_XDRAWAXES, 2)
+	    call gswind (gd, xmin / scale, xmax / scale, datalimit + threshold,
+		datalimit + threshold + ymax)
+	    call glabax (gd, Memc[str], "", "Intensity")
+	    call gseti (gd, G_YDRAWAXES, 0)
+	    call gseti (gd, G_XDRAWAXES, 1)
+	    call ggview (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, xmin, xmax, datalimit + threshold,
+		datalimit + threshold + ymax)
+	    call glabax (gd, "",
+	        "Radial Distance (lower-pixels, upper-scale units)", "")
+	} else {
+	    call gseti (gd, G_XDRAWAXES, 2)
+	    call gswind (gd, xmin / scale, xmax / scale, datalimit -
+	        threshold - ymax, datalimit - threshold)
+	    call glabax (gd, Memc[str], "", "Intensity")
+	    call gseti (gd, G_YDRAWAXES, 0)
+	    call gseti (gd, G_XDRAWAXES, 1)
+	    call ggview (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd,  xmin,  xmax, datalimit - threshold - ymax,
+	        datalimit - threshold)
+	    call glabax (gd, "",
+	        "Radial Distance (lower-pixels, upper-scale units)", "")
+	}
+
+	# Set the window up for plotting the data.
+	call gseti (gd, G_YDRAWAXES, 3)
+	call gseti (gd, G_XDRAWAXES, 3)
+	call gsetr (gd, G_ASPECT, aspect)
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_setr (gt, GTXMIN, xmin)
+	call gt_setr (gt, GTXMAX, xmax)
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gt_setr (gt, GTYMIN, ymin)
+	    call gt_setr (gt, GTYMAX, ymax)
+	} else {
+	    call gt_setr (gt, GTYMIN, ymin)
+	    call gt_setr (gt, GTYMAX, ymax)
+	}
+	call gt_swind (gd, gt)
+
+	call sfree (sp)
+end
+
+
+# AP_CPANNOTATE -- Procedure to annotate the radial plot in center.
+
+procedure ap_cpannotate (gd, ap)
+
+pointer	gd		# the graphics stream
+pointer	ap		# the apphot structure
+
+pointer	sp, str
+real	fwhmpsf, capert, datalimit, threshold, skysigma, cthreshold
+real	xmin, xmax, ymin, ymax
+int	apstati()
+real	apstatr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	call ggwind (gd, xmin, xmax, ymin, ymax)
+
+	fwhmpsf = 0.5 * apstatr (ap, FWHMPSF) * apstatr (ap, SCALE)
+	capert = 2.0 * fwhmpsf * apstatr (ap, CAPERT)
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	if (apstati (ap, POSITIVE) == YES)
+	    threshold = datalimit + threshold
+	else
+	    threshold = datalimit - threshold
+
+	# Plot the full width half maximum of the radial profile.
+	if (fwhmpsf >= xmin && fwhmpsf <= xmax) {
+	    call gamove (gd, fwhmpsf, ymin)
+	    call gadraw (gd, fwhmpsf, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "hwhm = %0.2f")
+	        call pargr (fwhmpsf)
+	    call gtext (gd, fwhmpsf, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the centering aperture.
+	if (capert >= xmin && capert <= xmax) {
+	    call gamove (gd, capert, ymin)
+	    call gadraw (gd, capert, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "cbox half-width = %0.2f")
+	        call pargr (capert)
+	    call gtext (gd, capert, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the threshold level for centering.
+        call sprintf (Memc[str], SZ_LINE, "threshold = %g")
+	    call pargr (threshold)
+	call gtext (gd, xmin, ymin, Memc[str], "q=h")
+
+	# Mark the sky sigma if defined.
+	if (! IS_INDEFR(skysigma) && skysigma >= ymin && skysigma <= ymax) {
+	    call gmark (gd, (xmin + xmax) / 2.0, (ymin + ymax) / 2.0, 
+		GM_VEBAR, -0.25, -skysigma)
+	    call sprintf (Memc[str], SZ_LINE, "sigma = %g")
+		call pargr (skysigma)
+	    call gtext (gd, (xmin + xmax) / 2.0, (ymin + ymax + skysigma) / 2.0,
+		Memc[str], "q=h;h=c")
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_CPRSET -- Procedure to reset the plot window after the data points
+# have been plotted.
+
+procedure ap_cpreset (gd, gt, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# the graphics stream
+pointer	gt		# the gtools pointer
+pointer	ap		# the apphot pointer
+real	xmin, xmax	# the minimum and maximum radial distance
+real	ymin, ymax	# the minimum and maximum of the y axis (ymin not used)
+
+real	datalimit, skysigma, threshold, cthreshold
+int	apstati()
+real	apstatr()
+
+begin
+	# Set the data window.
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	call gt_setr (gt, GTXMIN, xmin)
+	call gt_setr (gt, GTXMAX, xmax)
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gt_setr (gt, GTYMIN, datalimit + threshold)
+	    call gt_setr (gt, GTYMAX, ymax + datalimit + threshold)
+	} else {
+	    call gt_setr (gt, GTYMIN, datalimit - ymax - threshold)
+	    call gt_setr (gt, GTYMAX, datalimit - threshold)
+	}
+	call gt_swind (gd, gt)
+end