Initial commit

1 files changed, 181 insertions, 0 deletions

diff --git a/noao/onedspec/ecidentify/t_ecreid.x b/noao/onedspec/ecidentify/t_ecreid.x
new file mode 100644
index 00000000..5e9769ff
--- /dev/null
+++ b/noao/onedspec/ecidentify/t_ecreid.x
@@ -0,0 +1,181 @@
+include	<error.h>
+include <smw.h>
+include	"ecidentify.h"
+
+# T_ECREIDENTIFY -- Reidentify echelle features starting from reference.
+# If no initial shift is given then the procedure ec_shift computes a
+# shift between the reference features and the features in the image.
+# The purpose of the shift is to get the feature positions from the
+# reference image close enough to those of the image being identified
+# that the centering algorithm will determine the exact positions of the
+# features.  The recentered features are then fit with either a shift
+# of a full echelle function and written to database.
+
+procedure t_ecreidentify ()
+
+int	images			# List of images
+pointer	ref			# Reference image
+double	shift			# Initial shift
+
+int	i, j, fd, nfeatures1, nfeatures2
+double	shift1, pix, fit, pix_shift, fit_shift, z_shift
+pointer	sp, log, ec
+
+int	imtopenp(), ec_getim(), clpopnu(), clgfil(), open(), btoi()
+double	ec_fitpt(), ec_fittopix(), ec_shift(), ec_center(), ec_rms()
+double	clgetd()
+bool	clgetb()
+real	clgetr()
+errchk	ec_dbread(), ec_gdata(), ec_fitdata()
+
+begin
+	call smark (sp)
+	call salloc (ref, SZ_FNAME, TY_CHAR)
+	call salloc (log, SZ_FNAME, TY_CHAR)
+
+	# Allocate the basic data structure.
+	call ec_init (ec)
+
+	# Initialize fitting
+	call ecf_seti ("niterate", 0)
+	call ecf_setd ("low", 3.D0)
+	call ecf_setd ("high", 3.D0)
+
+	# Get task parameters.
+	images = imtopenp ("images")
+	call clgstr ("reference", Memc[ref], SZ_FNAME)
+	shift = clgetd ("shift")
+	call clgstr ("database", Memc[EC_DATABASE(ec)], SZ_FNAME)
+	EC_CRADIUS(ec) = clgetr ("cradius")
+	EC_THRESHOLD(ec) = clgetr ("threshold")
+	EC_LOGFILES(ec) = clpopnu ("logfiles")
+	EC_REFIT(ec) = btoi (clgetb ("refit"))
+
+	# Write logfile header.
+	while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) {
+	    iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    call sysid (Memc[log], SZ_LINE)
+	    call fprintf (fd, "\nECREIDENTIFY: %s\n")
+		call pargstr (Memc[log])
+	    call fprintf (fd,
+		"  Reference image = %s, Refit = %b\n")
+		call pargstr (Memc[ref])
+		call pargb (EC_REFIT(ec) == YES)
+	    call fprintf (fd, "%20s  %7s %7s %9s  %10s  %7s  %7s\n")
+		call pargstr ("Image")
+		call pargstr ("Found")
+		call pargstr ("Fit")
+		call pargstr ("Pix Shift")
+		call pargstr ("User Shift")
+		call pargstr ("Z Shift")
+		call pargstr ("RMS")
+	    call close (fd)
+	}
+
+	# Reidentify features in each spectrum.
+	while (ec_getim (images, Memc[EC_IMAGE(ec)], SZ_FNAME) != EOF) {
+	    call ec_gdata (ec)
+	    call ec_dbread (ec, Memc[ref], NO)
+	    call ec_fitdata (ec)
+	    call ec_fitfeatures (ec)
+
+	    if (IS_INDEFD (shift)) {
+		EC_FWIDTH(ec) = FWIDTH(ec,1)
+		EC_FTYPE(ec) = abs (FTYPE(ec,1))
+		EC_MINSEP(ec) = 1.
+		EC_MAXFEATURES(ec) = 20
+	        shift1 = ec_shift (ec)
+	    } else
+	        shift1 = shift
+
+	    # Recenter features.
+	    pix_shift = 0.
+	    fit_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) + shift1/ORDER(ec,i))
+	        pix = ec_center (ec, pix, FWIDTH(ec,i), FTYPE(ec,i))
+	        if (IS_INDEFD (pix))
+		    next
+	        fit = ec_fitpt (ec, APN(ec,i), pix)
+
+	        pix_shift = pix_shift + pix - PIX(ec,i)
+		fit_shift = fit_shift + (fit - FIT(ec,i)) * ORDER(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) = abs (FTYPE(ec,i))
+	    }
+	    EC_NFEATURES(ec) = j
+
+	    # If refitting the coordinate function is requested and there
+	    # is more than one feature and there is a previously defined
+	    # coordinate function then refit.  Otherwise compute a coordinate
+	    # shift.
+
+	    if ((EC_REFIT(ec)==YES)&&(EC_NFEATURES(ec)>1)&&(EC_ECF(ec)!=NULL)) {
+	        iferr (call ec_dofit (ec, NO, YES)) {
+		    call erract (EA_WARN)
+		    next
+		}
+	    } else
+	        call ec_doshift (ec, NO)
+	    if (EC_NEWECF(ec) == YES)
+	        call ec_fitfeatures (ec)
+
+	    nfeatures2 = 0
+	    do i = 1, EC_NFEATURES(ec)
+		if (FTYPE(ec,i) > 0)
+		    nfeatures2 = nfeatures2 + 1
+
+	    # Write a database entry for the reidentified image.
+	    if (EC_NFEATURES(ec) > 0)
+	        call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], NO)
+
+	    # Record log information if a log file descriptor is given.
+	    call clprew (EC_LOGFILES(ec))
+	    while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) {
+	        iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) {
+		    call erract (EA_WARN)
+		    next
+	        }
+	        call fprintf (fd,
+		    "%20s  %3d/%-3d %3d/%-3d %9.3g  %10.3g  %7.3g  %7.3g\n")
+		    call pargstr (Memc[EC_IMAGE(ec)])
+		    call pargi (EC_NFEATURES(ec))
+		    call pargi (nfeatures1)
+		    call pargi (nfeatures2)
+		    call pargi (EC_NFEATURES(ec))
+		    call pargd (pix_shift / max (1, EC_NFEATURES(ec)))
+		    call pargd (fit_shift / max (1, EC_NFEATURES(ec)))
+		    call pargd (z_shift / max (1, EC_NFEATURES(ec)))
+		    call pargd (ec_rms(ec))
+	        call close (fd)
+	    }
+
+	    call smw_close (MW(EC_SH(ec)))
+	    do i = 1, EC_NLINES(ec)
+		MW(SH(ec,i)) = NULL
+	}
+
+	call dgsfree (EC_ECF(ec))
+	call clpcls (EC_LOGFILES(ec))
+	call ec_free (ec)
+	call imtclose (images)
+	call sfree (sp)
+end