Initial commit

1 files changed, 273 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/phot/appplot.x b/noao/digiphot/apphot/phot/appplot.x
new file mode 100644
index 00000000..5172533d
--- /dev/null
+++ b/noao/digiphot/apphot/phot/appplot.x
@@ -0,0 +1,273 @@
+include <pkg/gtools.h>
+include <gset.h>
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/center.h"
+include "../lib/fitsky.h"
+include "../lib/phot.h"
+
+# AP_PPLOT -- Procedure to compute radial profile plots for the centering
+# routine.
+
+procedure ap_pplot (ap, im, sid, gd, makeplot)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# pointer to the iraf image
+int	sid		# id number of the star
+pointer	gd		# graphics stream
+int	makeplot	# make a plot ?
+
+int	apert, nx, ny
+pointer	buf, sp, str, r, gt
+real	xcenter, ycenter, xc, yc, rmin, rmax, imin, imax
+real	u1, u2, v1, v2, x1, x2, y1, y2
+int	ap_ctrpix()
+pointer	ap_gtinit()
+real	apstatr()
+
+begin
+	# Initialize
+	if (gd == NULL || makeplot == NO)
+	    return
+
+	# Check for defined center and get the pixels.
+	xcenter = apstatr (ap, XCENTER)
+	ycenter = apstatr (ap, YCENTER)
+	if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter))
+	    return
+
+	# Fetch the pixels.
+	apert = 2 * int (apstatr (ap, SCALE) * (apstatr (ap, ANNULUS) +
+	    apstatr (ap, DANNULUS))) + 1
+	buf = ap_ctrpix (im, xcenter, ycenter, apert, xc, yc, nx, ny)
+	if (buf == NULL)
+	    return
+
+	# Reactivate the work station.
+	call greactivate (gd, 0)
+
+	# 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[buf], nx * ny, imin, imax)
+
+	# Store the 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))
+	call gclear (gd)
+
+	# Label and annotate the plot.
+	call ap_ppset (gd, gt, ap, rmin, rmax, imin, imax)
+	call ap_ppannotate (gd, ap, rmin, rmax, imin, imax)
+
+	# Plot the coordinates.
+	call ap_plotrad (gd, gt, Memr[r], Memr[buf], nx * ny, "plus")
+
+	# Store the viewport and window coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	call ap_gtfree (gt)
+	call gdeactivate (gd, 0)
+	call sfree (sp)
+end
+
+
+# AP_PPSET -- Procedure to set up the parameters for the phot radial profile
+# plot.
+
+procedure ap_ppset (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 axes
+
+int	fd, naperts
+pointer	sp, str, title, temp
+real	aspect, scale, vx1, vx2, vy1, vy2
+int	stropen(), apstati()
+real	apstatr(), gstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, 5 * SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	naperts = apstati (ap, NAPERTS)
+	call salloc (temp, naperts, TY_REAL)
+
+	fd = stropen (Memc[str], 5 * SZ_LINE, WRITE_ONLY)
+	call sysid (Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+
+	# Encode the center parameter string.
+	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))
+
+	# Encode the sky fitting parameter string
+	call fprintf (fd,
+	    "Sky: value=%0.2f sigma=%0.2f skew=%0.2f nsky=%d nrej=%d\n")
+	    call pargr (apstatr (ap, SKY_MODE))
+	    call pargr (apstatr (ap, SKY_SIGMA))
+	    call pargr (apstatr (ap, SKY_SKEW))
+	    call pargi (apstati (ap, NSKY))
+	    call pargi (apstati (ap, NSKY_REJECT))
+
+	# Encode the apertures and magnitudes.
+	call ap_arrayr (ap, APERTS, Memr[temp])
+	call amulkr (Memr[temp], apstatr (ap, SCALE), Memr[temp], naperts)
+	call fprintf (fd, "Photometry: maxapert=")
+	call fprintf (fd, "%0.2f  mag=")
+	    call pargr (Memr[temp+naperts-1])
+	call ap_arrayr (ap, MAGS, Memr[temp])
+	call fprintf (fd, "%0.3f  merr=")
+	    call pargr (Memr[temp+naperts-1])
+	call ap_arrayr (ap, MAGERRS, Memr[temp])
+	call fprintf (fd, "%0.3f\n")
+	    call pargr (Memr[temp+naperts-1])
+
+	# Encode the title.
+	call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+
+	call strclose (fd)
+
+	aspect = gstatr (gd, G_ASPECT)
+	call gsetr (gd, G_ASPECT, 0.70)
+	scale = apstatr (ap, SCALE)
+
+	# Draw three axes.
+	call gseti (gd, G_XDRAWAXES, 2)
+	call gswind (gd, xmin / scale, xmax / scale, ymin, ymax)
+	call glabax (gd, Memc[str], "", "Intensity")
+
+	# Draw the bottom x axis.
+	call gseti (gd, G_YDRAWAXES, 0)
+	call gseti (gd, G_XDRAWAXES, 1)
+	call ggview (gd, vx1, vx2, vy1, vy2)
+	call gsview (gd, vx1, vx2, vy1, vy2)
+	call gswind (gd, xmin, xmax, ymin, ymax)
+	call glabax (gd,
+	    "","Radial Distance (lower-pixels, upper-scale units)", "")
+
+	# Restore the default draw axis parameters.
+	call gseti (gd, G_YDRAWAXES, 3)
+	call gseti (gd, G_XDRAWAXES, 3)
+	call gsetr (gd, G_ASPECT, aspect)
+
+	# Set the mark type.
+	call gt_sets (gt, GTTYPE, "mark")
+
+	call sfree (sp)
+end
+
+
+# AP_PPANNOTATE -- Procedure to annotate the radial plot in phot.
+
+procedure ap_ppannotate (gd, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# graphics stream
+pointer	ap		# apphot structure
+real	xmin, xmax	# minimum and maximum of the x axis
+real	ymin, ymax	# minimum and maximum of the y axis
+
+int	i, naperts
+pointer	sp, str, temp
+real	annulus, dannulus, skyval, skysigma
+int	apstati()
+real	apstatr ()
+
+begin
+	naperts = apstati (ap, NAPERTS)
+	call smark (sp)
+	call salloc (str, SZ_LINE + 1, TY_CHAR)
+	call salloc (temp, naperts, TY_REAL)
+
+	# Define some temporary variables
+	annulus = apstatr (ap, SCALE) * apstatr (ap, ANNULUS)
+	dannulus = annulus + apstatr (ap, SCALE) * apstatr (ap, DANNULUS)
+
+	# Mark the inner sky annulus.
+	if (annulus >= xmin && annulus <= xmax) {
+	    call gamove (gd, annulus, ymin)
+	    call gadraw (gd, annulus, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "inner sky radius = %0.2f")
+		call pargr (annulus)
+	    call gtext (gd, annulus, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+     
+	# Mark the outer sky annulus.
+	if (dannulus >= xmin && dannulus <= xmax) {
+	    call gamove (gd, dannulus, ymin)
+	    call gadraw (gd, dannulus, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "outer sky radius = %0.2f")
+	        call pargr (dannulus)
+	    call gtext (gd, dannulus, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the sky value.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	skyval = apstatr (ap, SKY_MODE)
+	if (skyval >= ymin && skyval <= ymax) {
+	    call gamove (gd, xmin, skyval)
+	    call gadraw (gd, xmax, skyval)
+	}
+
+	# Mark the upper sky sigma.
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	if (! IS_INDEFR(apstatr (ap, SKY_SIGMA)))
+	    skysigma = skyval + apstatr (ap, SHIREJECT) * apstatr (ap,
+	        SKY_SIGMA)
+	else
+	    skysigma = INDEFR
+	if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) {
+	    call gamove (gd, xmin, skysigma)
+	    call gadraw (gd, xmax, skysigma)
+	}
+
+	# Mark the lower sky sigma
+	if (! IS_INDEFR(apstatr (ap, SKY_SIGMA)))
+	    skysigma = skyval - apstatr (ap, SLOREJECT) * apstatr (ap,
+	        SKY_SIGMA)
+	else
+	    skysigma = INDEFR
+	if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) {
+	    call gamove (gd, xmin, skysigma)
+	    call gadraw (gd, xmax, skysigma)
+	}
+
+	# Mark the appertures.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	call ap_arrayr (ap, APERTS, Memr[temp])
+	call amulkr (Memr[temp], apstatr (ap, SCALE), Memr[temp], naperts)
+	do i = 1, naperts {
+	    call gamove (gd, Memr[temp+i-1], ymin)
+	    call gadraw (gd, Memr[temp+i-1], ymax)
+	    call sprintf (Memc[str], SZ_LINE, "apert[%d] = %0.2f")
+		call pargi (i)
+		call pargr (Memr[temp+i-1])
+	    call gtext (gd, Memr[temp+i-1], ymax, Memc[str],
+	        "q=h;u=180;v=t;p=r")
+	}
+
+	call sfree (sp)
+end