Initial commit

1 files changed, 486 insertions, 0 deletions

diff --git a/noao/imred/ccdred/src/cosmic/crexamine.x b/noao/imred/ccdred/src/cosmic/crexamine.x
new file mode 100644
index 00000000..d84961bc
--- /dev/null
+++ b/noao/imred/ccdred/src/cosmic/crexamine.x
@@ -0,0 +1,486 @@
+include	<error.h>
+include	<syserr.h>
+include <imhdr.h>
+include	<gset.h>
+include	<mach.h>
+include	<pkg/gtools.h>
+include	"crlist.h"
+
+# CR_EXAMINE  -- Examine cosmic ray candidates interactively.
+# CR_GRAPH    -- Make a graph
+# CR_NEAREST  -- Find the nearest cosmic ray to the cursor.
+# CR_DELETE   -- Set replace flag for cosmic ray candidate nearest cursor.
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+# CR_UPDATE   -- Change replacement flags, thresholds, and graphs.
+# CR_PLOT     -- Make log plot
+
+define	HELP	"noao$lib/scr/cosmicrays.key"
+define	PROMPT	"cosmic ray options"
+
+# CR_EXAMINE -- Examine cosmic ray candidates interactively.
+
+procedure cr_examine (cr, gp, gt, im, fluxratio, first)
+
+pointer	cr			# Cosmic ray list
+pointer	gp			# GIO pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+int	first			# Initial key
+
+char	cmd[SZ_LINE]
+int	i, newgraph, wcs, key, nc, nl, c1, c2, l1, l2, show
+real	wx, wy
+pointer	data
+
+int	clgcur()
+pointer	imgs2r()
+
+begin
+	# Set up the graphics.
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+
+	# Set image limits
+	nc = IM_LEN(im, 1)
+	nl = IM_LEN(im, 2)
+
+	# Enter cursor loop.
+	key = first
+	repeat {
+	    switch (key) {
+	    case '?': # Print help text.
+		call gpagefile (gp, HELP, PROMPT)
+	    case ':': # Colon commands.
+		switch (cmd[1]) {
+		case '/':
+		    call gt_colon (cmd, gp, gt, newgraph)
+		default:
+		    call printf ("\007")
+		}
+	    case 'a': # Toggle show all
+		if (show == 0)
+		    show = 1
+		else
+		    show = 0
+		newgraph = YES
+	    case 'd': # Delete candidate
+		call cr_delete (gp, wx, wy, cr, i, show)
+	    case 'q': # Quit
+		break
+	    case 'r': # Redraw the graph.
+		newgraph = YES
+	    case 's': # Make surface plots
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		c1 = max (1, int (Memr[CR_COL(cr)+i-1]) - 5)
+		c2 = min (nc, int (Memr[CR_COL(cr)+i-1]) + 5)
+		l1 = max (1, int (Memr[CR_LINE(cr)+i-1]) - 5)
+		l2 = min (nl, int (Memr[CR_LINE(cr)+i-1]) + 5)
+		data = imgs2r (im, c1, c2, l1, l2)
+		call gclear (gp)
+		call gsview (gp, 0.03, 0.48, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -33., 25.)
+		call gsview (gp, 0.53, 0.98, 0.53, 0.98)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -123., 25.)
+		call gsview (gp, 0.03, 0.48, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 57., 25.)
+		call gsview (gp, 0.53, 0.98, 0.03, 0.48)
+		call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 147., 25.)
+		call fprintf (STDERR, "[Type any key to continue]")
+		i = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+		newgraph = YES
+	    case 't': # Set threshold
+		call cr_update (gp, wy, cr, fluxratio, show)
+		call clputr ("fluxratio", fluxratio)
+	    case 'u': # Undelete candidate
+		call cr_undelete (gp, wx, wy, cr, i, show)
+	    case 'w':# Window the graph.
+		call gt_window (gt, gp, "cursor", newgraph)
+	    case ' ': # Print info
+		call cr_nearest (gp, wx, wy, cr, i, show)
+		call printf ("%d %d\n")
+		    call pargr (Memr[CR_COL(cr)+i-1])
+		    call pargr (Memr[CR_LINE(cr)+i-1])
+	    case 'z': # NOP
+		newgraph = NO
+	    default: # Ring bell for unrecognized commands.
+		call printf ("\007")
+	    }
+
+	    # Update the graph if needed.
+	    if (newgraph == YES) {
+		call cr_graph (gp, gt, cr, fluxratio, show)
+	        newgraph = NO
+	    }
+	} until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+end
+
+
+# CR_GRAPH -- Make a graph
+
+procedure cr_graph (gp, gt, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+pointer	gt		# GTOOLS pointers
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x1, x2, y1, y2
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	call gclear (gp)
+	call gt_ascale (gp, gt, Memr[x+1], Memr[y+1], ncr)
+	call gt_swind (gp, gt)
+	call gt_labax (gp, gt)
+
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_PLUS, 2., 2.)
+	    else
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_CROSS, 2., 2.)
+	    if (Memr[w+i] != 0.)
+		call gmark (gp, Memr[x+i], Memr[y+i], GM_BOX, 2., 2.)
+	}
+
+	call ggwind (gp, x1, x2, y1, y2)
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	call sfree (sp)
+end
+
+
+# CR_NEAREST -- Find the nearest cosmic ray to the cursor.
+
+procedure cr_nearest (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+	if (index != NULL)
+	    nearest = Memi[index+nearest]
+
+	# Move the cursor to the selected point.
+	call gscur (gp, x2, y2)
+
+	call sfree (sp)
+end
+
+
+# CR_DELETE -- Set replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_delete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the deleted point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSYES
+	    Memi[CR_WT(cr)+nearest-1] = -1
+	    call gscur (gp, x2, y2)
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_undelete (gp, wx, wy, cr, nearest, show)
+
+pointer	gp		# GIO pointer
+real	wx, wy		# Cursor position
+pointer	cr		# Cosmic ray list
+int	nearest		# Index of nearest point (returned)
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	sp, x, y, w, flag, index
+
+begin
+	call smark (sp)
+
+	call cr_show (show, cr, x, y, w, flag, index, ncr)
+	if (ncr == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Search for nearest point in NDC.
+	nearest = 0
+	r2min = MAX_REAL
+	call gctran (gp, wx, wy, wx, wy, 1, 0)
+	do i = 1, ncr {
+	    if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    call gctran (gp, x1, y1, x0, y0, 1, 0)
+	    r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	    if (r2 < r2min) {
+		r2min = r2
+		x2 = x1
+		y2 = y1
+		nearest = i
+	    }
+	}
+
+	# Move the cursor to the selected point and mark the delete point.
+	if (nearest > 0) {
+	    if (index != NULL)
+		nearest = Memi[index+nearest]
+	    Memi[CR_FLAG(cr)+nearest-1] = ALWAYSNO
+	    Memi[CR_WT(cr)+nearest-1] = 1
+	    call gscur (gp, x2, y2)
+
+	    call gseti (gp, G_PMLTYPE, 0)
+	    y2 = Memr[CR_RATIO(cr)+nearest-1]
+	    call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+	    call gseti (gp, G_PMLTYPE, 1)
+	    call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+	}
+
+	call sfree (sp)
+end
+
+
+# CR_UPDATE -- Change replacement flags, thresholds, and graphs.
+
+procedure cr_update (gp, wy, cr, fluxratio, show)
+
+pointer	gp		# GIO pointer
+real	wy		# Y cursor position
+pointer	cr		# Cosmic ray list
+real	fluxratio	# Flux ratio threshold
+int	show		# Show (0=all, 1=train)
+
+int	i, ncr, flag
+real	x1, x2, y1, y2
+pointer	x, y, f
+
+begin
+	call gseti (gp, G_PLTYPE, 0)
+	call ggwind (gp, x1, x2, y1, y2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+	fluxratio = wy
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, fluxratio, x2, fluxratio)
+
+	if (show == 1)
+	    return
+
+	ncr = CR_NCR(cr)
+	x = CR_FLUX(cr) - 1
+	y = CR_RATIO(cr) - 1
+	f = CR_FLAG(cr) - 1
+
+	do i = 1, ncr {
+	    flag = Memi[f+i]
+	    if ((flag == ALWAYSYES) || (flag == ALWAYSNO))
+		next
+	    x1 = Memr[x+i]
+	    y1 = Memr[y+i]
+	    if (flag == NO) {
+	        if (y1 < fluxratio) {
+		    Memi[f+i] = YES
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		}
+	    } else {
+	        if (y1 >= fluxratio) {
+		    Memi[f+i] = NO
+		    call gseti (gp, G_PMLTYPE, 0)
+		    call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+		    call gseti (gp, G_PMLTYPE, 1)
+		    call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+		}
+	    }
+	}
+end
+
+
+# CR_PLOT -- Make log plot
+
+procedure cr_plot (cr, im, fluxratio)
+
+pointer	cr			# Cosmic ray list
+pointer	im			# Image pointer
+real	fluxratio		# Flux ratio threshold
+
+int	fd, open(), errcode()
+pointer	sp, fname, gp, gt, gopen(), gt_init()
+errchk	gopen
+
+begin
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Open the plotfile.
+	call clgstr ("plotfile", Memc[fname], SZ_FNAME)
+	iferr (fd = open (Memc[fname], APPEND, BINARY_FILE)) {
+	    if (errcode() != SYS_FNOFNAME)
+		call erract (EA_WARN)
+	    return
+	}
+
+	# Set up the graphics.
+	gp = gopen ("stdplot", NEW_FILE, fd)
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTXTRAN, "log")
+	call gt_setr (gt, GTXMIN, 10.)
+	call gt_setr (gt, GTYMIN, 0.)
+	call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+	call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+	call gt_sets (gt, GTXLABEL, "Flux")
+	call gt_sets (gt, GTYLABEL, "Flux Ratio")
+
+	call cr_graph (gp, gt, cr, fluxratio, 'r')
+
+	call gt_free (gt)
+	call gclose (gp)
+	call close (fd)
+	call sfree (sp)
+end
+
+
+# CR_SHOW -- Select data to show.
+# This returns pointers to the data.  Note the pointers are salloc from
+# the last smark which is done by the calling program.
+
+procedure cr_show (show, cr, x, y, w, flag, index, ncr)
+
+int	show		#I Data to show (0=all, 1=train)
+pointer	cr		#I CR data
+pointer	x		#O Fluxes
+pointer	y		#O Ratios
+pointer	w		#O Weights
+pointer	flag		#O Flags
+pointer	index		#O Index into CR data (if not null)
+int	ncr		#O Number of selected data points
+
+int	i
+
+begin
+	switch (show) {
+	case 0:
+	    ncr = CR_NCR(cr)
+	    x = CR_FLUX(cr) - 1
+	    y = CR_RATIO(cr) - 1
+	    w = CR_WT(cr) - 1
+	    flag = CR_FLAG(cr) - 1
+	    index = NULL
+	case 1:
+	    ncr = CR_NCR(cr)
+	    call salloc (x, ncr, TY_REAL)
+	    call salloc (y, ncr, TY_REAL)
+	    call salloc (w, ncr, TY_REAL)
+	    call salloc (flag, ncr, TY_INT)
+	    call salloc (index, ncr, TY_INT)
+
+	    ncr = 0
+	    x = x - 1
+	    y = y - 1
+	    w = w - 1
+	    flag = flag - 1
+	    index = index - 1
+
+	    do i = 1, CR_NCR(cr) {
+		if (Memr[CR_WT(cr)+i-1] == 0.)
+		    next
+		ncr = ncr + 1
+		Memr[x+ncr] = Memr[CR_FLUX(cr)+i-1]
+		Memr[y+ncr] = Memr[CR_RATIO(cr)+i-1]
+		Memr[w+ncr] = Memr[CR_WT(cr)+i-1]
+		Memi[flag+ncr] = Memi[CR_FLAG(cr)+i-1]
+		Memi[index+ncr] = i
+	    }
+	}
+end