Initial commit

1 files changed, 1092 insertions, 0 deletions

diff --git a/noao/rv/rvidlines/t_reidentify.x b/noao/rv/rvidlines/t_reidentify.x
new file mode 100644
index 00000000..5a724325
--- /dev/null
+++ b/noao/rv/rvidlines/t_reidentify.x
@@ -0,0 +1,1092 @@
+include	<error.h>
+include	<fset.h>
+include	<gset.h>
+include	<pkg/gtools.h>
+include	<smw.h>
+include	"identify.h"
+
+define	ICFITHELP	"noao$lib/scr/idicgfit.key"
+define	VLIGHT	2.997925e5	# Speed of light, Km/sec
+
+
+# T_REIDENTIFY -- Reidentify features starting from reference features.
+
+procedure t_reidentify ()
+
+begin
+	call reiden (IDENTIFY)
+end
+
+
+# T_RVREIDLINES -- Reidentify lines for radial velocities
+
+procedure t_rvreidlines ()
+
+begin
+	call reiden (RVIDLINES)
+end
+
+
+# REIDEN -- Reidentify features starting from reference features.
+# A reference spectrum is specified and the same features are identified
+# in other images.  Some lines may be lost due to bad centering.  Additional
+# lines may be excluded from a new fit to the dispersion function.  Instead
+# of refitting the dispersion function the user may elect to determine only
+# a shift in the reference dispersion function.  Additional features may
+# be added given a coordinate list.
+#
+# In 2D images a starting line or column is selected.  A number of lines
+# or columns may be averaged before identifying features.  If a positive step
+# size is given then additional lines or columns may be reidentified in
+# the reference image.  This may be done either by tracing or by reidentifying
+# starting from the same reference features.  Reidentification between images
+# is done by taking the same line or column from the reference image.
+# 
+# Multispec format images are matched by aperture number and the spectra
+# need not be in the same order in each image.
+
+procedure reiden (taskid)
+
+int	taskid			#I Task ID
+
+pointer	reference		# Reference image
+int	list			# List of images
+char	ans[3]			# Interactive?
+
+int	i, fd, nlogfd
+pointer	sp, logfile, str, id, logfd, pd
+
+int	clscan(), clgeti(), clpopnu(), clgfil(), clgwrd(), btoi()
+int	nscan(), open(), nowhite(), imtopenp(), imtgetim()
+bool	clgetb()
+real	clgetr()
+pointer	gopen(), gt_init()
+
+begin
+	call smark (sp)
+	call salloc (reference, SZ_FNAME, TY_CHAR)
+	call salloc (logfile, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Allocate the basic data structures.
+	call id_init (id, taskid)
+	call ic_open (ID_IC(id))
+
+	# Get task parameters.
+	call clgstr ("reference", Memc[reference], SZ_FNAME)
+	list = imtopenp ("images")
+	i = nowhite (Memc[reference], Memc[reference], SZ_FNAME)
+
+	if (ID_TASK(id) == IDENTIFY)
+	    ID_REFIT(id) = btoi (clgetb ("refit"))
+	else
+	    ID_REFIT(id) = 3
+
+	if (clscan ("nsum") != EOF) {
+	    call gargi (ID_NSUM(id,1))
+	    call gargi (ID_NSUM(id,2))
+	    if (nscan() == 0)
+		call error (1, "Error in 'nsum' parameter")
+	    if (nscan() == 1)
+		ID_NSUM(id,2) = ID_NSUM(id,1)
+	    ID_NSUM(id,1) = max (1, ID_NSUM(id,1))
+	    ID_NSUM(id,2) = max (1, ID_NSUM(id,2))
+	}
+	ID_MAXFEATURES(id) = clgeti ("maxfeatures")
+	ID_MINSEP(id) = clgetr ("minsep")
+	ID_MATCH(id) = clgetr ("match")
+	if (ID_TASK(id) == IDENTIFY) {
+	    ID_ZWIDTH(id) = clgetr ("identify.zwidth")
+	    ID_FTYPE(id) = clgwrd ("identify.ftype", Memc[str], SZ_LINE, FTYPES)
+	    ID_FWIDTH(id) = clgetr ("identify.fwidth")
+	} else {
+	    ID_ZWIDTH(id) = clgetr ("rvidlines.zwidth")
+	    ID_FTYPE(id) = clgwrd ("rvidlines.ftype", Memc[str], SZ_LINE,FTYPES)
+	    ID_FWIDTH(id) = clgetr ("rvidlines.fwidth")
+	}
+	ID_CRADIUS(id) = clgetr ("cradius")
+	ID_THRESHOLD(id) = clgetr ("threshold")
+	call clgstr ("database", Memc[ID_DATABASE(id)], SZ_FNAME)
+	call clgstr ("coordlist", Memc[ID_COORDLIST(id)], SZ_FNAME)
+	ID_LABELS(id) = 1
+
+	call id_mapll (id)
+	ID_LOGFILES(id) = clpopnu ("logfiles")
+
+	switch (clgwrd ("interactive", ans, SZ_FNAME, "|no|yes|NO|YES|")) {
+	case 1, 3:
+	    call strcpy ("NO", ans, 3)
+	    ID_GP(id) = NULL
+	case 2, 4:
+	    # Open graphics
+	    call clgstr ("graphics", Memc[logfile], SZ_FNAME)
+	    ID_GP(id) = gopen (Memc[logfile], NEW_FILE+AW_DEFER, STDGRAPH)
+	    if (ID_TASK(id) == IDENTIFY) {
+		call ic_pstr (ID_IC(id), "help", ICFITHELP)
+		call ic_pstr (ID_IC(id), "xlabel", "Feature positions")
+		call ic_pstr (ID_IC(id), "xunits", "pixels")
+		call ic_pstr (ID_IC(id), "ylabel", "")
+		call ic_pkey (ID_IC(id), 1, 'y', 'x')
+		call ic_pkey (ID_IC(id), 2, 'y', 'v')
+		call ic_pkey (ID_IC(id), 3, 'y', 'r')
+		call ic_pkey (ID_IC(id), 4, 'y', 'd')
+		call ic_pkey (ID_IC(id), 5, 'y', 'n')
+		call ic_puti (ID_IC(id), "key", 3)
+	    }
+	}
+
+	# Open log and plot files.
+	nlogfd = 0
+	if (clgetb ("verbose")) {
+	    nlogfd = 1
+	    call malloc (logfd, nlogfd, TY_INT)
+	    Memi[logfd] = STDOUT
+	}
+	while (clgfil (ID_LOGFILES(id), Memc[logfile], SZ_FNAME) != EOF) {
+	    fd = open (Memc[logfile], APPEND, TEXT_FILE)
+	    call fseti (fd, F_FLUSHNL, YES)
+	    nlogfd = nlogfd + 1
+	    if (nlogfd == 1)
+		call malloc (logfd, nlogfd, TY_INT)
+	    else
+		call realloc (logfd, nlogfd, TY_INT)
+	    Memi[logfd+nlogfd-1] = fd
+	}
+	call ri_loghdr (id, Memc[reference], Memi[logfd], nlogfd, 1)
+
+	pd = NULL
+	if (ID_TASK(id) == IDENTIFY) {
+	    call clgstr ("plotfile", Memc[logfile], SZ_FNAME)
+	    if (nowhite (Memc[logfile], Memc[logfile], SZ_FNAME) > 0) {
+		fd = open (Memc[logfile], APPEND, BINARY_FILE)
+		pd = gopen ("stdvdm", NEW_FILE, fd)
+	    }
+	}
+
+	ID_GT(id) = gt_init()
+	call gt_sets (ID_GT(id), GTTYPE, "line")
+
+	# Get and trace the reference solutions.
+	call ri_reference (id, Memc[reference], ans, Memi[logfd], nlogfd, pd)
+
+	# Expand the image template and reidentify features.
+	while (imtgetim (list, Memc[ID_IMAGE(id)], SZ_FNAME) != EOF)
+	    call ri_image (id, Memc[reference], Memc[ID_IMAGE(id)], ans,
+		Memi[logfd], nlogfd, pd)
+
+	# Finish up.
+	if (nlogfd > 0) {
+	    do i = 1, nlogfd
+	        call close (Memi[logfd+i-1])
+	    call mfree (logfd, TY_INT)
+	}
+
+	if (ID_GP(id) != NULL)
+	    call gclose (ID_GP(id))
+	if (pd != NULL) {
+	    call gclose (pd)
+	    call close (fd)
+	}
+	call clpcls (ID_LOGFILES(id))
+	call imtclose (list)
+	call id_free (id)
+	call smw_daxis (NULL, NULL, 0, 0, 0)
+	call sfree (sp)
+end
+
+
+# RI_REFERENCE -- Set reference features.  Trace if needed.
+
+procedure ri_reference (id, reference, ans, logfd, nlogfd, pd)
+
+pointer	id			# ID pointer
+char	reference[ARB]		# Reference image
+char	ans[3]			# Interactive?
+int	logfd[ARB]		# Logfiles
+int	nlogfd			# Number of logfiles
+pointer	pd			# Plot file pointer
+
+int	step[2]
+double	shift[2]
+int	nreid
+bool	override
+bool	trace
+
+int	i, start[2], line[2], loghdr
+double	fit_shift[2]
+pointer	ic, ic1
+bool	clgetb()
+int	clscan(), clgeti(), nscan(), id_gid(), id_dbcheck()
+errchk	id_dbread
+
+begin
+	# Open the image and return if there is an error.
+	call strcpy (reference, Memc[ID_IMAGE(id)], SZ_FNAME)
+	iferr (call id_map (id)) {
+	    call erract (EA_WARN)
+	    return
+	}
+
+	# Set parameters
+	start[1] = ID_LINE(id,1)
+	start[2] = ID_LINE(id,2)
+	if (clscan ("step") != EOF) {
+	    call gargi (step[1])
+	    call gargi (step[2])
+	    if (nscan() == 0)
+		call error (1, "Error in 'step' parameter")
+	    if (nscan() == 1)
+		step[2] = step[1]
+	}
+	if (clscan ("shift") != EOF) {
+	    call gargd (shift[1])
+	    call gargd (shift[2])
+	    if (nscan() == 0)
+		call error (1, "Error in 'shift' parameter")
+	    if (nscan() == 1)
+		shift[2] = shift[1]
+	}
+	nreid = max (1, ID_NFEATURES(id) - clgeti ("nlost"))
+	override = clgetb ("override")
+	trace = clgetb ("trace")
+
+	# Get and save the reference database entry.
+	call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO, NO)
+	call id_saveid (id, ID_LINE(id,1))
+
+	# Get and save other entries.
+	if (!override) {
+	    for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) {
+		ID_LINE(id,2) = line[2]
+		ID_AP(id,2) = line[2]
+		for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+		    if (line[1]==start[1] && line[2]==start[2])
+			next
+		    ID_LINE(id,1) = line[1]
+		    ID_AP(id,1) = line[1]
+		    if (ID_APS(id) != NULL)
+			ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		    ifnoerr (
+			call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+			    NO, NO)) {
+			call id_saveid (id, ID_LINE(id,1))
+		    }
+		}
+		for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+		    line[1]=line[1]+step[1]) {
+		    ID_LINE(id,1) = line[1]
+		    ID_AP(id,1) = line[1]
+		    if (ID_APS(id) != NULL)
+			ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		    ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)],
+			ID_AP(id,1), NO, NO)) {
+			call id_saveid (id, ID_LINE(id,1))
+		    }
+		}
+	    }
+	    for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2);
+		line[2]=line[2]+step[2]) {
+		ID_LINE(id,2) = line[2]
+		ID_AP(id,2) = line[2]
+		for (line[1]=start[1]-step[1]; line[1]>0;
+		    line[1]=line[1]-step[1]) {
+		    ID_LINE(id,1) = line[1]
+		    ID_AP(id,1) = line[1]
+		    if (ID_APS(id) != NULL)
+			ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		    ifnoerr (
+			call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+			    NO, NO)) {
+			call id_saveid (id, ID_LINE(id,1))
+		    }
+		}
+		for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+		    line[1]=line[1]+step[1]) {
+		    ID_LINE(id,1) = line[1]
+		    ID_AP(id,1) = line[1]
+		    if (ID_APS(id) != NULL)
+			ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		    ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)],
+			ID_AP(id,1), NO, NO)) {
+			call id_saveid (id, ID_LINE(id,1))
+		    }
+		}
+	    }
+	}
+
+	# Reidentify.
+	loghdr = 2
+	ic = ID_IC(id)
+	if (ans[1] == 'N')
+	    ic1 = ic
+	else {
+	    call ic_open (ic1)
+	    call ic_copy (ic, ic1)
+	}
+
+	fit_shift[2] = shift[2]
+	for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) {
+	    ID_LINE(id,2) = line[2]
+	    ID_AP(id,2) = line[2]
+	    ID_IC(id) = ic
+
+	    if (IS_INDEFD(shift[2]))
+		fit_shift[2] = INDEFD
+	    else {
+		if (!trace)
+		    fit_shift[2] = fit_shift[2] - shift[2]
+		else
+		    fit_shift[2] = -shift[2]
+	    }
+
+	    fit_shift[1] = fit_shift[2]
+	    for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+		if (line[1]==start[1] && line[2]==start[2])
+		    next
+		ID_LINE(id,1) = line[1]
+		ID_AP(id,1) = line[1]
+		ID_IC(id) = ic
+		if (ID_APS(id) != NULL)
+		    ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+
+		if (!trace) {
+		    i = id_gid (id, start)
+		    ID_LINE(id,1) = line[1]
+		    ID_LINE(id,2) = line[2]
+		}
+
+		if (IS_INDEFD(shift[1]))
+		    fit_shift[1] = INDEFD
+		else {
+		    if (!trace)
+			fit_shift[1] = fit_shift[1] - shift[1]
+		    else
+			fit_shift[1] = -shift[1]
+		}
+
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) < nreid && trace) {
+		    call ri_loghdr (id, reference, logfd, nlogfd, 3)
+		    break
+		}
+
+		if (ID_NFEATURES(id) > 0) {
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		    call id_saveid (id, line)
+		}
+	    }
+
+	    ID_IC(id) = ic
+	    i = id_gid (id, start)
+	    fit_shift[1] = fit_shift[2]
+	    for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+		line[1]=line[1]+step[1]) {
+		ID_LINE(id,1) = line[1]
+		ID_AP(id,1) = line[1]
+		ID_IC(id) = ic
+		if (ID_APS(id) != NULL)
+		    ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+
+		if (!trace) {
+		    i = id_gid (id, start)
+		    ID_LINE(id,1) = line[1]
+		    ID_LINE(id,2) = line[2]
+		}
+
+		if (IS_INDEFD(shift[1]))
+		    fit_shift[1] = INDEFD
+		else {
+		    if (!trace)
+			fit_shift[1] = fit_shift[1] + shift[1]
+		    else
+			fit_shift[1] = shift[1]
+		}
+
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) < nreid && trace) {
+		    call ri_loghdr (id, reference, logfd, nlogfd, 3)
+		    break
+		}
+
+		if (ID_NFEATURES(id) > 0) {
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		    call id_saveid (id, line)
+		}
+	    }
+	}
+
+
+	fit_shift[2] = 0.
+	for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2);
+	    line[2]=line[2]+step[2]) {
+	    ID_LINE(id,2) = line[2]
+	    ID_AP(id,2) = line[2]
+	    ID_IC(id) = ic
+
+	    if (IS_INDEFD(shift[2]))
+		fit_shift[2] = INDEFD
+	    else {
+		if (!trace)
+		    fit_shift[2] = fit_shift[2] + shift[2]
+		else
+		    fit_shift[2] = shift[2]
+	    }
+
+	    fit_shift[1] = fit_shift[2]
+	    for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+		ID_LINE(id,1) = line[1]
+		ID_AP(id,1) = line[1]
+		ID_IC(id) = ic
+		if (ID_APS(id) != NULL)
+		    ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+
+		if (!trace) {
+		    i = id_gid (id, start)
+		    ID_LINE(id,1) = line[1]
+		    ID_LINE(id,2) = line[2]
+		}
+
+		if (IS_INDEFD(shift[1]))
+		    fit_shift[1] = INDEFD
+		else {
+		    if (!trace)
+			fit_shift[1] = fit_shift[1] - shift[1]
+		    else
+			fit_shift[1] = -shift[1]
+		}
+
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) < nreid && trace) {
+		    call ri_loghdr (id, reference, logfd, nlogfd, 3)
+		    break
+		}
+
+		if (ID_NFEATURES(id) > 0) {
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		    call id_saveid (id, line)
+		}
+	    }
+
+	    ID_IC(id) = ic
+	    i = id_gid (id, start)
+	    fit_shift[1] = fit_shift[2]
+	    for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+		line[1]=line[1]+step[1]) {
+		ID_LINE(id,1) = line[1]
+		ID_AP(id,1) = line[1]
+		ID_IC(id) = ic
+		if (ID_APS(id) != NULL)
+		    ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+
+		if (!trace) {
+		    i = id_gid (id, start)
+		    ID_LINE(id,1) = line[1]
+		    ID_LINE(id,2) = line[2]
+		}
+
+		if (IS_INDEFD(shift[1]))
+		    fit_shift[1] = INDEFD
+		else {
+		    if (!trace)
+			fit_shift[1] = fit_shift[1] + shift[1]
+		    else
+			fit_shift[1] = shift[1]
+		}
+
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) < nreid && trace) {
+		    call ri_loghdr (id, reference, logfd, nlogfd, 3)
+		    break
+		}
+
+		if (ID_NFEATURES(id) > 0) {
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		    call id_saveid (id, line)
+		}
+	    }
+	}
+
+	ID_IC(id) = ic
+	if (ic != ic1)
+	    call ic_closed (ic1)
+
+	call smw_close (MW(ID_SH(id)))
+	call imunmap (IM(ID_SH(id)))
+	call shdr_close (ID_SH(id))
+end
+
+
+# RI_IMAGE -- Reidentify an image.
+
+procedure ri_image (id, reference, image, ans, logfd, nlogfd, pd)
+
+pointer	id			# ID pointer
+char	reference[ARB]		# Reference image
+char	image[ARB]		# Image to be reidentified
+char	ans[3]			# Interactive?
+int	logfd[ARB]		# Logfiles
+int	nlogfd			# Number of logfiles
+pointer	pd			# Plot file pointer
+
+bool	newaps			# Add new apertures not in reference?
+bool	override		# Override previous identifications?
+bool	verbose			# Verbose output?
+
+int	i, j, ap, loghdr, id_getid(), id_dbcheck()
+double	shift, fit_shift, clgetd()
+pointer	ic, ic1
+bool	clgetb()
+
+begin
+	# Open the image and return if there is an error.
+	call strcpy (image, Memc[ID_IMAGE(id)], SZ_FNAME)
+	iferr (call id_map (id)) {
+	    call erract (EA_WARN)
+	    return
+	}
+	call dtunmap (ID_DT(id))
+
+	newaps = clgetb ("newaps")
+	override = clgetb ("override")
+	verbose = clgetb ("verbose")
+
+	ic = ID_IC(id)
+	if (ans[1] == 'N')
+	    ic1 = ic
+	else
+	    call ic_open (ic1)
+
+	loghdr = 2
+	shift = clgetd ("shift")
+
+	# For MULTISPEC search the reference list of each aperture.  If
+	# a reference of the same aperture is not found and the newaps
+	# flag is set use the initial reference and then add the
+	# reidentification to the reference list.
+	# For NDSPEC apply each reference to the image.
+
+	if (SMW_FORMAT(MW(ID_SH(id))) == SMW_ES ||
+	    SMW_FORMAT(MW(ID_SH(id))) == SMW_MS) {
+	    for (i=1; i<=ID_MAXLINE(id,1); i=i+1) {
+		ap = Memi[ID_APS(id)+i-1]
+		for (j=ID_NID(id); id_getid (id,j)!=EOF; j=j-1)
+		    if (ap == ID_AP(id,1))
+			break
+
+		if (j == 0 && !newaps) {
+		    if (verbose) {
+			call printf (
+			    "%s: Reference for aperture %d not found\n")
+			    call pargstr (image)
+			    call pargi (ap)
+		    }
+		    next
+		}
+
+		ID_AP(id,1) = ap
+		ID_LINE(id,1) = i
+
+		if (i == 1 && ic != ic1)
+		    call ic_copy (ic, ic1)
+
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+		
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+
+		fit_shift = shift
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) > 0) {
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		    if (j == 0 && newaps) {
+			call id_sid (id, ID_NID(id)+1)
+			if (verbose) {
+			    call printf (
+				"%s: New reference for aperture %d\n")
+				call pargstr (image)
+				call pargi (ap)
+			} 
+		    }
+		}
+		ID_IC(id) = ic
+	    }
+
+	} else {
+	    for (i=1; id_getid (id,i)!=EOF; i=i+1) {
+		if (i == 1 && ic != ic1)
+		    call ic_copy (ic, ic1)
+
+		if (!override)
+		    if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+			next
+		
+		ID_IC(id) = ic1
+		call id_gdata (id)
+		iferr (call id_fitdata (id))
+		    ;
+
+		call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+		loghdr = 0
+
+		fit_shift = shift
+		call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+		if (ID_NFEATURES(id) > 0)
+		    call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+		ID_IC(id) = ic
+	    }
+	}
+
+	ID_IC(id) = ic
+	if (ic != ic1)
+	    call ic_closed (ic1)
+	call smw_close (MW(ID_SH(id)))
+	call imunmap (IM(ID_SH(id)))
+	call shdr_close (ID_SH(id))
+end
+
+
+# RI_REIDENTIFY -- Reidentify features using a reference image database entry.
+
+procedure ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+pointer	id			# ID pointer
+double	fit_shift		# Shift in fit coords (input and output)
+char	ans[3]			# Interactive?
+int	logfd[ARB]		# Logfiles
+int	nlogfd			# Number of logfiles
+pointer	pd			# Plot file pointer
+
+int	i, j, nfeatures1, nfeatures2, nfit, iden, mono, clgwrd()
+double	shift, pix_shift, z_shift, v, vrms
+double	id_fitpt(), fit_to_pix(), id_shift(), id_center(), id_rms(), id_zval()
+pointer	sp, pix, fit
+bool	clgetb()
+
+begin
+	nfeatures1 = ID_NFEATURES(id)
+	call smark (sp)
+	call salloc (pix, nfeatures1, TY_DOUBLE)
+	call salloc (fit, nfeatures1, TY_DOUBLE)
+	call amovd (PIX(id,1), Memd[pix], nfeatures1)
+	call amovd (FIT(id,1), Memd[fit], nfeatures1)
+
+	# If no initial shift is given then the procedure id_shift
+	# computes a shift between the reference features and the
+	# features in the image.  The purpose of the shift is to get the
+	# reference feature positions close enough to those of the image
+	# being identified that the centering algorithm will determine
+	# the exact positions of the features.  An initial shift of zero
+	# is used if the two images are very nearly aligned as in the
+	# case of tracing features in a two dimensional image or for a
+	# set of images taken with the same observing setup.
+
+	if (IS_INDEFD(fit_shift)) {
+	    ID_FWIDTH(id) = FWIDTH(id,1)
+	    ID_FTYPE(id) = FTYPE(id,1)
+	    shift = id_shift (id)
+	} else
+	    shift = fit_shift
+
+	# For each reference feature a shift is added to bring the pixel
+	# position near that for the image being identified and then the
+	# centering algorithm is used.  If the centering algorithm fails
+	# the feature is discarded.  A mean shift is computed for the
+	# features which have been reidentified.
+
+	do i = 1, ID_NFEATURES(id) {
+	    PIX(id,i) = fit_to_pix (id, FIT(id,i) + shift)
+	    PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i), FTYPE(id,i),
+		NO)
+	    if (!IS_INDEFD(PIX(id,i)))
+	        FIT(id,i) = id_fitpt (id, PIX(id,i))
+	}
+	for (i=1; i<ID_NFEATURES(id); i=i+1) {
+	    if (IS_INDEFD(PIX(id,i)))
+		next
+	    for (j=i+1; j<=ID_NFEATURES(id); j=j+1) {
+	        if (IS_INDEFD(PIX(id,j)))
+		    next
+		if (abs (PIX(id,i)-PIX(id,j)) < ID_MINSEP(id)) {
+		    if (abs (FIT(id,i)-USER(id,i)) < abs (FIT(id,j)-USER(id,j)))
+			PIX(id,j) = INDEFD
+		    else {
+			PIX(id,i) = INDEFD
+			break
+		    }
+		}
+	    }
+	}
+
+	pix_shift = 0.
+	fit_shift = 0.
+	z_shift = 0.
+	j = 0
+	do i = 1, ID_NFEATURES(id) {
+	    if (IS_INDEFD(PIX(id,i)))
+		next
+
+	    pix_shift = pix_shift + PIX(id,i) - Memd[pix+i-1]
+	    fit_shift = fit_shift + FIT(id,i) - Memd[fit+i-1]
+	    if (Memd[fit+i-1] != 0.)
+		z_shift = z_shift + id_zval (id, FIT(id,i), Memd[fit+i-1])
+
+	    j = j + 1
+	    PIX(id,j) = PIX(id,i)
+	    FIT(id,j) = FIT(id,i)
+	    USER(id,j) = USER(id,i)
+	    WTS(id,j) = WTS(id,i)
+	    FWIDTH(id,j) = FWIDTH(id,i)
+	    FTYPE(id,j) = FTYPE(id,i)
+	}
+	ID_NFEATURES(id) = j
+		    
+	nfeatures2 = j
+	pix_shift = pix_shift / max (1, ID_NFEATURES(id))
+	fit_shift = fit_shift / max (1, ID_NFEATURES(id))
+	z_shift = z_shift / max (1, ID_NFEATURES(id))
+
+	# 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.
+
+	mono = YES
+	switch (ID_REFIT(id)) {
+	case 1:
+	    if (ID_CV(id) != NULL) {
+		if (ID_NFEATURES(id)>1)
+		    call id_dofit (id, NO)
+		else
+		    call id_doshift (id, NO)
+	    }
+	case 2:
+	    if (ID_CV(id) != NULL)
+		call id_doshift (id, NO)
+	case 3:
+	    call id_velocity (id, NO)
+	    v = ID_REDSHIFT(id) * VLIGHT
+	    vrms = ID_RMSRED(id) * VLIGHT
+	}
+	if (ID_NEWCV(id) == YES) {
+	    iferr (call id_fitdata (id))
+		mono = NO
+	    call id_fitfeatures (id)
+	}
+
+	if (clgetb ("addfeatures")) {
+	    ID_FWIDTH(id) = FWIDTH(id,1)
+	    ID_FTYPE(id) = FTYPE(id,1)
+	    call id_linelist (id)
+	    if (ID_NEWFEATURES(id) == YES) {
+		switch (ID_REFIT(id)) {
+		case 1:
+		    if (ID_NFEATURES(id)>1)
+			call id_dofit (id, NO)
+		    else
+			call id_doshift (id, NO)
+		case 2:
+		    call id_doshift (id, NO)
+		case 3:
+		    call id_velocity (id, NO)
+		    v = ID_REDSHIFT(id) * VLIGHT
+		    vrms = ID_RMSRED(id) * VLIGHT
+		}
+		if (ID_NEWCV(id) == YES) {
+	    	    iferr (call id_fitdata (id))
+			mono = NO
+		    call id_fitfeatures (id)
+		}
+	    }
+	}
+
+	# Enter fitting interactively.
+	iden = NO
+	if ((ID_NFEATURES(id)>1) && (ID_CV(id)!=NULL || ID_REFIT(id) == 3)) {
+	    if (ans[1] != 'N') {
+	        if (ans[1] != 'Y') {
+		    nfit = 0
+		    for (j=1; j<=ID_NFEATURES(id); j=j+1)
+		        if (WTS(id,j) > 0.)
+			    nfit = nfit + 1
+	            call printf (
+		    "%s%s%23t%3d/%-3d %3d/%-3d %9.3g  %10.3g  %7.3g  %7.3g\n")
+			call pargstr (Memc[ID_IMAGE(id)])
+			call pargstr (Memc[ID_SECTION(id)])
+			call pargi (nfeatures2)
+			call pargi (nfeatures1)
+			call pargi (nfit)
+			call pargi (ID_NFEATURES(id))
+			call pargd (pix_shift)
+			call pargd (fit_shift)
+			if (ID_REFIT(id) == 3) {
+			    call pargd (v)
+			    call pargd (vrms)
+			} else {
+			    call pargd (z_shift)
+			    call pargd (id_rms(id))
+			}
+	            call flush (STDOUT)
+		    repeat {
+	                ifnoerr (i = clgwrd ("answer", ans, SZ_FNAME,
+			    "|no|yes|NO|YES|"))
+			    break
+		    }
+		    call clpstr ("answer", ans)
+		}
+	        switch (ans[1]) {
+	        case 'y', 'Y':
+		    mono = YES
+		    i = ID_REFIT(id)
+		    call reidentify (id)
+		    ID_REFIT(id) = i
+		    if (ID_REFIT(id) == 3) {
+			call id_velocity (id, NO)
+			v = ID_REDSHIFT(id) * VLIGHT
+			vrms = ID_RMSRED(id) * VLIGHT
+		    }
+		    iden = YES
+		}
+		if (ans[1] != 'Y')
+		    call gdeactivate (ID_GP(id), 0)
+	    }
+	}
+
+	# Record log information if a log file descriptor is given.
+	for (i = 1; i <= nlogfd; i = i + 1) {
+	    if (ans[1] == 'n' && logfd[i] == STDOUT)
+		    next
+	    nfit = 0
+	    for (j=1; j<=ID_NFEATURES(id); j=j+1)
+		if (WTS(id,j) > 0.)
+		    nfit = nfit + 1
+	    call fprintf (logfd[i],
+		"%s%s%23t%3d/%-3d %3d/%-3d %9.3g  %10.3g  %7.3g  %7.3g\n")
+		call pargstr (Memc[ID_IMAGE(id)])
+		call pargstr (Memc[ID_SECTION(id)])
+		call pargi (nfeatures2)
+		call pargi (nfeatures1)
+		call pargi (nfit)
+		call pargi (ID_NFEATURES(id))
+		call pargd (pix_shift)
+		call pargd (fit_shift)
+		if (ID_REFIT(id) == 3) {
+		    call pargd (v)
+		    call pargd (vrms)
+		} else {
+		    call pargd (z_shift)
+		    call pargd (id_rms(id))
+		}
+	    if (ID_REFIT(id) == 3 && logfd[i] != STDOUT)
+		call id_log (id, "", logfd[i])
+	    if (mono == NO)
+		call fprintf (logfd[i], "Non-monotonic dispersion function")
+	    call flush (logfd[i])
+	    if (logfd[i] == STDOUT)
+		iden = NO
+	}
+	# Print log if STDOUT is not used but if the IDENTIFY is done.
+	if (iden == YES) {
+	    call printf (
+		"%s%s%23t%3d/%-3d %3d/%-3d %9.3g  %10.3g  %7.3g  %7.3g\n")
+		call pargstr (Memc[ID_IMAGE(id)])
+		call pargstr (Memc[ID_SECTION(id)])
+		call pargi (nfeatures2)
+		call pargi (nfeatures1)
+		call pargi (nfit)
+		call pargi (ID_NFEATURES(id))
+		call pargd (pix_shift)
+		call pargd (fit_shift)
+		if (ID_REFIT(id) == 3) {
+		    call pargd (v)
+		    call pargd (vrms)
+		} else {
+		    call pargd (z_shift)
+		    call pargd (id_rms(id))
+		}
+	    if (mono == NO)
+		call printf ("Non-monotonic dispersion function")
+	    call flush (STDOUT)
+	}
+		
+	# Make log plot.
+	call ri_plot (id, pd)
+
+	call sfree (sp)
+end
+
+
+# RI_LOGHDR -- Print a log header in the log files.
+
+procedure ri_loghdr (id, reference, logfd, nlogfd, flag)
+
+pointer	id		# Identify structure
+char	reference[ARB]	# Reference image
+int	logfd[ARB]	# Log file descriptors
+int	nlogfd		# Number of log files
+int	flag		# Header type flag (1=banner, 2=Column labels, 3=Error)
+
+int	i
+pointer	str
+
+begin
+	for (i = 1; i <= nlogfd; i = i + 1) {
+	    switch (flag) {
+	    case 1:	# Print ID
+		call malloc (str, SZ_LINE, TY_CHAR)
+	        call sysid (Memc[str], SZ_LINE)
+		switch (ID_REFIT(id)) {
+		case 1, 2:
+		    call fprintf (logfd[i], "\nREIDENTIFY: %s\n")
+			call pargstr (Memc[str])
+		case 3:
+		    call fprintf (logfd[i], "\nRVREIDLINES: %s\n")
+			call pargstr (Memc[str])
+		}
+		call mfree (str, TY_CHAR)
+	    case 2:	# Print labels
+		switch (ID_REFIT(id)) {
+		case 1, 2:
+		    call fprintf (logfd[i],
+			"  Reference image = %s, New image = %s, Refit = %b\n")
+			call pargstr (reference)
+			call pargstr (Memc[ID_IMAGE(id)])
+			call pargi (ID_REFIT(id))
+		    call fprintf (logfd[i],
+			"%20s  %7s %7s %9s  %10s  %7s  %7s\n")
+			call pargstr ("Image Data")
+			call pargstr ("Found")
+			call pargstr ("Fit")
+			call pargstr ("Pix Shift")
+			call pargstr ("User Shift")
+			call pargstr ("Z Shift")
+			call pargstr ("RMS")
+		case 3:
+		    call fprintf (logfd[i],
+			"  Reference image = %s, New image = %s\n")
+			call pargstr (reference)
+			call pargstr (Memc[ID_IMAGE(id)])
+		    call fprintf (logfd[i],
+			"%20s  %7s %7s %9s  %10s %8s  %7s\n")
+			call pargstr ("Image Data")
+			call pargstr ("Found")
+			call pargstr ("Fit")
+			call pargstr ("Pix Shift")
+			call pargstr ("User Shift")
+			call pargstr ("Velocity")
+			call pargstr ("RMS")
+		}
+	    case 3:	# Error
+		call fprintf (logfd[i], "    ** Too many features lost **\n")
+	    }
+	}
+end
+
+
+# RI_PLOT -- Plot residual graph of reidentified lines.
+
+procedure ri_plot (id, pd)
+
+pointer	id				# ID pointer
+pointer	pd				# GIO pointer
+
+int	i, j
+pointer	sp, str, x, y, gt, gt_init()
+double	id_zshiftd()
+
+begin
+	# Check if there is anything to plot.
+	if (pd == NULL || ID_NFEATURES(id) == 0)
+	    return
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (x, ID_NFEATURES(id), TY_REAL)
+	call salloc (y, ID_NFEATURES(id), TY_REAL)
+
+	# Set plot points.
+	j = 0
+	do i = 1, ID_NFEATURES(id) {
+	    if (IS_INDEFD(USER(id,i)))
+		break
+
+	    Memr[x+j] = USER(id,i)
+	    Memr[y+j] = id_zshiftd (id, FIT(id,i), 0)
+	    j = j + 1
+	}
+
+	if (j == 0) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Make the plot.
+	call sprintf (Memc[str], SZ_LINE, "Reidentify: %s")
+	    call pargstr (Memc[ID_IMAGE(id)])
+	gt = gt_init ()
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTXLABEL, "user coordinates")
+	call gt_sets (gt, GTYLABEL, "residuals (fit - user)")
+	call gt_sets (gt, GTTITLE, Memc[str])
+	call gclear (pd)
+	call gascale (pd, Memr[x], j, 1)
+	call gascale (pd, Memr[y], j, 2)
+	call gt_swind (pd, gt)
+	call gt_labax (pd, gt)
+	call gt_plot (pd, gt, Memr[x], Memr[y], j)
+	call gt_free (gt)
+
+	call sfree (sp)
+end