Initial commit

1 files changed, 535 insertions, 0 deletions

diff --git a/noao/onedspec/ecidentify/ecidentify.x b/noao/onedspec/ecidentify/ecidentify.x
new file mode 100644
index 00000000..827568d1
--- /dev/null
+++ b/noao/onedspec/ecidentify/ecidentify.x
@@ -0,0 +1,535 @@
+include	<error.h>
+include	<imhdr.h>
+include	<gset.h>
+include	<smw.h>
+include	"ecidentify.h"
+
+define	HELP		"noao$onedspec/ecidentify/ecidentify.key"
+define	PROMPT		"ecidentify options"
+
+define	PAN		1	# Pan graph
+define	ZOOM		2	# Zoom graph
+
+# EC_IDENTIFY -- Identify echelle features in an image.
+# This is the basic interactive loop.
+
+procedure ec_identify (ec)
+
+pointer	ec			# EC pointer
+
+real	wx, wy
+int	wcs, key
+char	cmd[SZ_LINE]
+
+char	newimage[SZ_FNAME]
+int	i, j, last, all, prfeature, nfeatures1, npeaks
+bool	answer
+double	pix, fit, user, shift, pix_shift, z_shift
+pointer	peaks
+
+bool	clgetb()
+int	clgcur(), scan(), nscan(), find_peaks(), ec_next(), ec_previous()
+int	ec_line()
+double	ec_center(), ec_fittopix(), ec_fitpt(), ec_shift(), ec_rms()
+double	ecf_getd()
+errchk	ec_gdata(), ec_graph(), ec_dbread(), xt_mk1d(), ec_line()
+
+define	newim_		10
+define	newkey_		20
+define	beep_		99
+
+begin
+newim_	# Start here for each new image.
+
+	# Get the image data.  Return if there is an error.
+	iferr (call ec_gdata (ec)) {
+	    call erract (EA_WARN)
+	    return
+	}
+
+	# Look for a database entry for the image.
+	iferr {
+	    call ec_dbread (ec, Memc[EC_IMAGE(ec)], NO)
+	    EC_NEWDBENTRY(ec) = NO
+	} then
+	    if ((EC_NFEATURES(ec) > 0) || (EC_ECF(ec) != NULL))
+	        EC_NEWDBENTRY(ec) = YES
+
+	# Set the coordinate array and the feature data.
+	iferr (call ec_fitdata (ec))
+	    call erract (EA_WARN)
+	call ec_fitfeatures (ec)
+
+	# Begin with the first image line.
+	EC_LINE(ec) = 1
+	EC_AP(ec) = APS(ec,EC_LINE(ec))
+	EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec))
+	call ec_gline (ec, EC_LINE(ec))
+
+	# Initialize.
+	EC_GTYPE(ec) = PAN
+	EC_REFIT(ec) = NO
+	EC_NEWFEATURES(ec) = NO
+	EC_NEWECF(ec) = NO
+	EC_CURRENT(ec) = 0
+	i = ec_next (ec, EC_CURRENT(ec))
+	last = EC_CURRENT(ec)
+	all = 0
+	newimage[1] = EOS
+	key = 'r'
+
+	repeat {
+	    prfeature = YES
+	    if (all != 0)
+		all = mod (all + 1, 3)
+
+	    switch (key) {
+	    case '?':	# Page help
+		call gpagefile (EC_GP(ec), HELP, PROMPT)
+	    case ':':	# Execute colon commands
+		if (cmd[1] == '/')
+		    call gt_colon (cmd, EC_GP(ec), EC_GT(ec), EC_NEWGRAPH(ec))
+		else
+		    call ec_colon (ec, cmd, newimage, prfeature)
+	    case ' ':	# Go to the current feature
+	    case '.':	# Go to the nearest feature
+		if (EC_NFEATURES(ec) == 0)
+		    goto beep_
+		call ec_nearest (ec, double (wx))
+	    case '-':	# Go to the previous feature
+		if (ec_previous (ec, EC_CURRENT(ec)) == EOF)
+		    goto beep_
+	    case '+', 'n':	# Go to the next feature
+		if (ec_next (ec, EC_CURRENT(ec)) == EOF)
+		    goto beep_
+	    case 'a':	# Set the all flag for the next key
+		all = 1
+	    case 'c':	# Center features on data
+		if (all != 0) {
+		    call eprintf ("Recentering features ...\n")
+		    for (i = 1; i <= EC_NFEATURES(ec); i = i + 1) {
+			call ec_gline (ec, LINE(ec,i))
+		        call gseti (EC_GP(ec), G_PLTYPE, 0)
+		        call ec_mark (ec, i)
+		        call gseti (EC_GP(ec), G_PLTYPE, 1)
+			FWIDTH(ec,i) = EC_FWIDTH(ec)
+		        PIX(ec,i) = ec_center (ec, PIX(ec,i), FWIDTH(ec,i),
+			    FTYPE(ec,i))
+		        if (!IS_INDEFD (PIX(ec,i))) {
+			    FIT(ec,i) = ec_fitpt (ec, APN(ec,i), PIX(ec,i))
+		            call ec_mark (ec, i)
+			} else {
+			    call ec_delete (ec, i)
+			    i = i - 1
+		        }
+		    }
+		    call ec_gline (ec, EC_LINE(ec))
+		    EC_NEWFEATURES(ec) = YES
+		} else {
+		    if (EC_NFEATURES(ec) == 0)
+		       goto beep_
+
+		    call ec_nearest (ec, double (wx))
+		    pix = PIX(ec,EC_CURRENT(ec))
+		    pix = ec_center (ec, pix, EC_FWIDTH(ec),
+			FTYPE(ec,EC_CURRENT(ec)))
+		    if (!IS_INDEFD (pix)) {
+		        call gseti (EC_GP(ec), G_PLTYPE, 0)
+		        call ec_mark (ec, EC_CURRENT(ec))
+			PIX(ec,EC_CURRENT(ec)) = pix
+		        FWIDTH(ec,EC_CURRENT(ec)) = EC_FWIDTH(ec)
+			FIT(ec,EC_CURRENT(ec)) =
+			    ec_fitpt (ec, APN(ec,EC_CURRENT(ec)), pix)
+		        call gseti (EC_GP(ec), G_PLTYPE, 1)
+			call ec_mark (ec, EC_CURRENT(ec))
+		        EC_NEWFEATURES(ec) = YES
+		    } else {
+			call eprintf ("Centering failed\n")
+			prfeature = NO
+		    }
+		}
+	    case 'd':	# Delete features
+		if (all != 0) {
+		    EC_NFEATURES(ec) = 0
+		    EC_CURRENT(ec) = 0
+		    EC_NEWFEATURES(ec) = YES
+		    EC_NEWGRAPH(ec) = YES
+		} else {
+		    if (EC_NFEATURES(ec) == 0)
+			goto beep_
+
+		    call ec_nearest (ec, double (wx))
+		    call gseti (EC_GP(ec), G_PLTYPE, 0)
+		    call ec_mark (ec, EC_CURRENT(ec))
+		    call gseti (EC_GP(ec), G_PLTYPE, 1)
+		    call ec_delete (ec, EC_CURRENT(ec))
+		    call ec_nearest (ec, double (wx))
+		    last = 0
+		}
+	    case 'f':	# Fit dispersion function
+		iferr (call ec_dofit (ec, YES, NO)) {
+		    call erract (EA_WARN)
+		    prfeature = NO
+		    goto beep_
+		}
+	    case 'g':	# Fit shift
+		call ec_doshift (ec, YES)
+		prfeature = NO
+	    case 'i':	# Initialize
+	        call dgsfree (EC_ECF(ec))
+		call ecf_setd ("shift", 0.D0)
+	        EC_NEWECF(ec) = YES
+		EC_NFEATURES(ec) = 0
+		EC_CURRENT(ec) = 0
+		EC_NEWFEATURES(ec) = YES
+		EC_NEWGRAPH(ec) = YES
+	    case 'j':	# Go to the previous order
+		EC_LINE(ec) =
+		    mod (EC_LINE(ec)+EC_NLINES(ec)-2, EC_NLINES(ec)) + 1
+		EC_AP(ec) = APS(ec,EC_LINE(ec))
+		EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec))
+		call ec_gline (ec, EC_LINE(ec))
+		EC_NEWGRAPH(ec) = YES
+		EC_CURRENT(ec) = 0
+		i = ec_next (ec, EC_CURRENT(ec))
+	    case 'k':	# Go to the next order
+		EC_LINE(ec) = mod (EC_LINE(ec), EC_NLINES(ec)) + 1
+		EC_AP(ec) = APS(ec,EC_LINE(ec))
+		EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec))
+		call ec_gline (ec, EC_LINE(ec))
+		EC_NEWGRAPH(ec) = YES
+		EC_CURRENT(ec) = 0
+		i = ec_next (ec, EC_CURRENT(ec))
+	    case 'l':	# Find features using a line list
+		if (EC_ECF(ec) == NULL) {
+		    call eprintf ("Doing initial fit ...\n")
+		    iferr (call ec_dofit (ec, NO, NO)) {
+		        call erract (EA_WARN)
+		        prfeature = NO
+			goto beep_
+		    }
+		    if (EC_NEWECF(ec) == YES) {
+		        iferr (call ec_fitdata (ec)) {
+			    call erract (EA_WARN)
+			    prfeature = NO
+		        }
+		        call ec_fitfeatures (ec)
+		        EC_NEWECF(ec) = NO
+		    }
+		}
+
+		call eprintf ("Searching coordinate list ...\n")
+		call ec_linelist (ec)
+		EC_CURRENT(ec) = 0
+		i = ec_next (ec, EC_CURRENT(ec))
+		if (EC_NEWFEATURES(ec) == YES)
+		    EC_NEWGRAPH(ec) = YES
+	    case 'm':	# Mark a new feature
+		fit = wx
+		pix = ec_fittopix (ec, fit)
+		pix = ec_center (ec, pix, EC_FWIDTH(ec), EC_FTYPE(ec))
+		if (IS_INDEFD (pix))
+		    goto beep_
+		fit = ec_fitpt (ec, EC_AP(ec), pix)
+		user = fit
+		call ec_newfeature (ec, EC_AP(ec), pix, fit, user,
+		    EC_FWIDTH(ec), EC_FTYPE(ec))
+		USER(ec,EC_CURRENT(ec)) = INDEFD
+		call ec_match (ec, FIT(ec,EC_CURRENT(ec)),
+		    USER(ec,EC_CURRENT(ec)))
+		call ec_mark (ec, EC_CURRENT(ec))
+		call printf ("%3d %10.2f %10.8g (%10.8g): ")
+		    call pargi (APN(ec,EC_CURRENT(ec)))
+		    call pargd (PIX(ec,EC_CURRENT(ec)))
+		    call pargd (FIT(ec,EC_CURRENT(ec)))
+		    call pargd (USER(ec,EC_CURRENT(ec)))
+		call flush (STDOUT)
+		if (scan() != EOF) {
+		    call gargd (user)
+		    if (nscan() == 1) {
+			USER(ec,EC_CURRENT(ec)) = user
+			call ec_match (ec, user, USER(ec,EC_CURRENT(ec)))
+		    }
+		}
+	    case 'o':	# Go to a specified order
+		call printf ("Aperture (%d): ")
+		    call pargi (EC_AP(ec))
+		    call flush (STDOUT)
+		if (scan() != EOF) {
+		    call gargi (j)
+		    if (nscan() == 1) {
+			if (j != EC_AP(ec)) {
+			    iferr {
+				i = ec_line (ec, j)
+				call ec_gline (ec, i)
+				EC_LINE(ec) = i
+				EC_AP(ec) = j
+				EC_ORDER(ec) = ORDERS(ec,i)
+				EC_NEWGRAPH(ec) = YES
+				EC_CURRENT(ec) = 0
+				i = ec_next (ec, EC_CURRENT(ec))
+			    } then
+				goto beep_
+			}
+		    }
+		}
+	    case 'p':	# Go to pan graph mode
+		if (EC_GTYPE(ec) == PAN)
+		    goto beep_
+
+		EC_GTYPE(ec) = PAN
+		EC_NEWGRAPH(ec) = YES
+	    case 'q':	# Quit
+		break
+	    case 'r':	# Redraw the current graph
+		EC_NEWGRAPH(ec) = YES
+	    case 's', 'x':	# Shift or cross correlate features
+		# Get coordinate shift.
+		switch (key) {
+		case 's':
+		    call printf ("User coordinate (%10.8g): ")
+		        call pargr (wx)
+		        call flush (STDOUT)
+		    if (scan() != EOF) {
+		        call gargd (user)
+		        if (nscan() == 1)
+			    shift = (wx - user) * EC_ORDER(ec)
+		    } else
+		        shift = 0.
+		case 'x':
+		    if (EC_NFEATURES(ec) > 5) {
+			call eprintf (
+			    "Cross correlating features with peaks ...\n")
+			shift = ec_shift (ec)
+		    } else
+			goto beep_
+		}
+
+		EC_NEWFEATURES(ec) = YES
+		EC_NEWECF(ec) = YES
+		EC_NEWGRAPH(ec) = YES
+		prfeature = NO
+
+		if (EC_NFEATURES(ec) < 1) {
+		    call printf ("User coordinate shift=%5f")
+			call pargd (shift / EC_ORDER(ec))
+		    call ecf_setd ("shift", ecf_getd ("shift") - shift)
+		    goto newkey_
+		}
+
+		# Recenter features.
+		call eprintf ("Recentering features ...\n")
+		pix_shift = 0.
+		z_shift = 0.
+		nfeatures1 = EC_NFEATURES(ec)
+
+		j = 0.
+		do i = 1, EC_NFEATURES(ec) {
+		    call ec_gline (ec, LINE(ec,i))
+		    pix = ec_fittopix (ec, FIT(ec,i) + shift/ORDER(ec,i))
+		    pix = ec_center (ec, pix, FWIDTH(ec,i), FTYPE(ec,i))
+		    if (IS_INDEFD (pix)) {
+			if (EC_CURRENT(ec) == i)
+			     EC_CURRENT(ec) = 0
+			next
+		    }
+		    fit = ec_fitpt (ec, APN(ec,i), pix)
+
+		    pix_shift = pix_shift + pix - PIX(ec,i)
+		    if (FIT(ec,i) != 0.)
+		        z_shift = z_shift + (fit - FIT(ec,i)) / FIT(ec,i)
+
+		    j = j + 1
+		    APN(ec,j) = APN(ec,i)
+		    LINE(ec,j) = LINE(ec,i)
+		    ORDER(ec,j) = ORDER(ec,i)
+		    PIX(ec,j) = pix
+		    FIT(ec,j) = FIT(ec,i)
+		    USER(ec,j) = USER(ec,i)
+		    FWIDTH(ec,j) = FWIDTH(ec,i)
+		    FTYPE(ec,j) = FTYPE(ec,i)
+		    if (EC_CURRENT(ec) == i)
+			EC_CURRENT(ec) = j
+		}
+		call ec_gline (ec, EC_LINE(ec))
+		EC_NFEATURES(ec) = j
+		if (EC_CURRENT(ec) == 0)
+		    i = ec_next (ec, EC_CURRENT(ec))
+
+		if (EC_NFEATURES(ec) < 1) {
+		    call printf ("User coordinate shift=%5f")
+			call pargd (shift / EC_ORDER(ec))
+		    call printf (", No features found during recentering")
+		    call ecf_setd ("shift", ecf_getd ("shift") - shift)
+		    goto newkey_
+		}
+
+		# Adjust shift.
+		pix = ecf_getd ("shift")
+		call ec_doshift (ec, NO)
+		call ec_fitfeatures (ec)
+
+		# Print results.
+		call printf ("Recentered=%d/%d")
+		    call pargi (EC_NFEATURES(ec))
+		    call pargi (nfeatures1)
+		call printf (
+		    ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g")
+		    call pargd (pix_shift / EC_NFEATURES(ec))
+		    call pargd ((pix - ecf_getd ("shift")) / EC_ORDER(ec))
+		    call pargd (z_shift / EC_NFEATURES(ec))
+		    call pargd (ec_rms(ec))
+	    case 't':	# Move current feature
+		if (EC_CURRENT(ec) == 0)
+		    goto beep_
+
+		call gseti (EC_GP(ec), G_PLTYPE, 0)
+		call ec_mark (ec, EC_CURRENT(ec))
+		pix = ec_fittopix (ec, double (wx))
+		PIX(ec,EC_CURRENT(ec)) = pix
+		FIT(ec,EC_CURRENT(ec)) =
+		    ec_fitpt (ec, APN(ec,EC_CURRENT(ec)), pix)
+		call gseti (EC_GP(ec), G_PLTYPE, 1)
+		call ec_mark (ec, EC_CURRENT(ec))
+		EC_NEWFEATURES(ec) = YES
+	    case 'u':	# Set uesr coordinate value
+		if (EC_NFEATURES(ec) == 0)
+		    goto beep_
+
+		call printf ("%3d %10.2f %10.8g (%10.8g): ")
+		    call pargi (APN(ec,EC_CURRENT(ec)))
+		    call pargd (PIX(ec,EC_CURRENT(ec)))
+		    call pargd (FIT(ec,EC_CURRENT(ec)))
+		    call pargd (USER(ec,EC_CURRENT(ec)))
+		    call flush (STDOUT)
+		if (scan() != EOF) {
+		    call gargd (user)
+		    if (nscan() == 1) {
+			USER(ec,EC_CURRENT(ec)) = user
+			EC_NEWFEATURES(ec) = YES
+		    }
+		}
+	    case 'w':	# Window graph
+		call gt_window (EC_GT(ec), EC_GP(ec), "cursor", EC_NEWGRAPH(ec))
+	    case 'y':	# Find peaks in order
+		call malloc (peaks, EC_NPTS(ec), TY_REAL)
+		npeaks = find_peaks (IMDATA(ec,1), Memr[peaks],
+		    EC_NPTS(ec), 0., int (EC_MINSEP(ec)), 0, EC_MAXFEATURES(ec),
+		    0., false)
+		for (j = 1; j <= EC_NFEATURES(ec); j = j + 1) {
+		    for (i = 1; i <= npeaks; i = i + 1) {
+			if (!IS_INDEF(pix)) {
+			    pix = Memr[peaks + i - 1]
+			    if (abs (pix - PIX(ec,j)) < EC_MINSEP(ec))
+			        Memr[peaks + i - 1] = INDEF
+			}
+		    }
+		}
+		for (i = 1; i <= npeaks; i = i + 1) {
+		    pix = Memr[peaks+i-1]
+		    pix = ec_center (ec, pix, EC_FWIDTH(ec), EC_FTYPE(ec))
+		    if (IS_INDEFD (pix))
+			next
+		    fit = ec_fitpt (ec, EC_AP(ec), pix)
+		    user = INDEFD
+		    call ec_match (ec, fit, user)
+		    call ec_newfeature (ec, EC_AP(ec), pix, fit, user,
+			EC_FWIDTH(ec), EC_FTYPE(ec))
+		    call ec_mark (ec, EC_CURRENT(ec))
+		}
+		call mfree (peaks, TY_REAL)
+	    case 'z':	# Go to zoom mode
+		if (EC_CURRENT(ec) == 0)
+		    goto beep_
+
+		if (EC_GTYPE(ec) == PAN)
+		    EC_NEWGRAPH(ec) = YES
+		EC_GTYPE(ec) = ZOOM
+		call ec_nearest (ec, double (wx))
+	    case 'I':	# Interrupt
+		call fatal (0, "Interrupt")
+	    default:	# Beep
+beep_		 call printf ("\007\n")
+	    }
+
+newkey_
+	    # Set database update flag if there has been a change.
+	    if ((EC_NEWFEATURES(ec) == YES) || (EC_NEWECF(ec) == YES))
+		EC_NEWDBENTRY(ec) = YES
+
+	    # Exit loop and then start new image.
+	    if (newimage[1] != EOS)
+		break
+
+	    # Refit the dispersion function if needed.
+	    if (EC_REFIT(ec) == YES) {
+		iferr (call ec_dofit (ec, NO, NO)) {
+		    call erract (EA_WARN)
+		    prfeature = NO
+		}
+		EC_REFIT(ec) = NO
+	    }
+
+	    # Recompute the coordinate information.
+	    if (EC_NEWECF(ec) == YES) {
+		iferr (call ec_fitdata (ec)) {
+		    call erract (EA_WARN)
+		    prfeature = NO
+		}
+		call ec_fitfeatures (ec)
+		EC_NEWECF(ec) = NO
+	    }
+
+	    # Redraw new feature in zoom mode.
+	    if ((EC_GTYPE(ec) == ZOOM) && (last != EC_CURRENT(ec)))
+		EC_NEWGRAPH(ec) = YES
+
+	    # Redraw graph.
+	    if (EC_NEWGRAPH(ec) == YES) {
+		call ec_graph (ec, EC_GTYPE(ec))
+	        EC_NEWGRAPH(ec) = NO
+	    }
+
+	    # Set cursor and print current feature on status (unless canceled).
+	    if (EC_CURRENT(ec) > 0) {
+		call gscur (EC_GP(ec), real (FIT(ec,EC_CURRENT(ec))), wy)
+		if (prfeature == YES) {
+	            call printf ("%d %10.2f %10.8g %10.8g\n")
+		        call pargi (APN(ec,EC_CURRENT(ec)))
+		        call pargd (PIX(ec,EC_CURRENT(ec)))
+		        call pargd (FIT(ec,EC_CURRENT(ec)))
+		        call pargd (USER(ec,EC_CURRENT(ec)))
+		}
+	    }
+
+	    last = EC_CURRENT(ec)
+	} until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+
+	# Warn user that feature data is newer than database entry.
+	if (EC_NEWDBENTRY(ec) == YES) {
+	    answer = true
+	    if (!clgetb ("autowrite")) {
+		call printf ("Write feature data to the database (yes)? ")
+		call flush (STDOUT)
+		if (scan() != EOF)
+		    call gargb (answer)
+	    }
+	    if (answer)
+	        call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], NO)
+	}
+
+	call flush (STDOUT)
+
+	# Free image data and MWCS
+	call mfree (EC_PIXDATA(ec), TY_DOUBLE)
+	call mfree (EC_FITDATA(ec), TY_DOUBLE)
+	call smw_close (MW(EC_SH(ec)))
+	do i = 1, EC_NLINES(ec)
+	    MW(SH(ec,i)) = NULL
+
+	# If a new image was specified by a colon command don't return.
+	if (newimage[1] != EOS) {
+	    call strcpy (newimage, Memc[EC_IMAGE(ec)], SZ_FNAME)
+	    goto newim_
+	}
+end