Initial commit

1 files changed, 1600 insertions, 0 deletions

diff --git a/noao/onedspec/identify/autoid/autoid.x b/noao/onedspec/identify/autoid/autoid.x
new file mode 100644
index 00000000..3f169ca7
--- /dev/null
+++ b/noao/onedspec/identify/autoid/autoid.x
@@ -0,0 +1,1600 @@
+include	<mach.h>
+include	<error.h>
+include	<smw.h>
+include	"../identify.h"
+include	"autoid.h"
+
+
+# AID_TARGET -- Select target lines and the dispersion limits to be searched.
+# The dispersion limits may be based on header parameters.
+
+procedure aid_target (aid)
+
+pointer	aid		#I AID pointer
+
+int	i, j, l, nt, n
+double	dw, dwmin, dwmax, pix, aid_imgd(), id_center()
+pointer	sp, x, y, idt, idr, im, xt, xtl
+int	id_upeaks(), stridxs()
+errchk	id_upeaks, id_center
+
+begin
+        call smark (sp)
+        call salloc (x, ID_NPTS(AID_IDT(aid)), TY_REAL)
+
+	idt = AID_IDT(aid)
+	idr = AID_IDR(aid)
+	im = IM(ID_SH(idt))
+	nt = ID_NPTS(idt)
+
+	# Set the approximate coordinate information.
+	AID_CRVAL(aid) = aid_imgd (im, AID_CR(aid))
+	AID_CDELT(aid) = aid_imgd (im, AID_CD(aid))
+	AID_CRPIX(aid) = aid_imgd (im, AID_CP(aid))
+	AID_CRQUAD(aid) = aid_imgd (im, AID_CQ(aid))
+	AID_CRSEARCH(aid) = aid_imgd (im, AID_CRS(aid))
+	AID_CDSEARCH(aid) = aid_imgd (im, AID_CDS(aid))
+
+	if (IS_INDEFD(AID_CRPIX(aid)))
+	    AID_CRPIX(aid) = (nt+1) / 2.
+
+	if (IS_INDEFD(AID_CRQUAD(aid)))
+	    AID_CRQUAD(aid) = 0D0
+
+	if (!IS_INDEFD(AID_CRVAL(aid)) && !IS_INDEFD(AID_CDELT(aid))) {
+	    dw = nt * AID_CDELT(aid)
+	    if (IS_INDEFD(AID_CRSEARCH(aid)))
+		AID_CRSEARCH(aid) = abs (0.1 * dw)
+	    else if (AID_CRSEARCH(aid) < 0.)
+		AID_CRSEARCH(aid) = abs (AID_CRSEARCH(aid) * dw)
+	    if (IS_INDEFD(AID_CDSEARCH(aid)))
+		AID_CDSEARCH(aid) = abs (0.1 * AID_CDELT(aid))
+	    else if (AID_CDSEARCH(aid) < 0.)
+		AID_CDSEARCH(aid) = abs (AID_CDSEARCH(aid) * AID_CDELT(aid))
+	    AID_CRSEARCH(aid) = max (abs (0.0001 * dw),
+		AID_CRSEARCH(aid))
+	    AID_CDSEARCH(aid) = max (abs (0.0001 * AID_CDELT(aid)),
+		AID_CDSEARCH(aid))
+	    dwmax = 2 * AID_CRSEARCH(aid) + (nt - 1) *
+		(abs (AID_CDELT(aid)) + AID_CDSEARCH(aid))
+	    dwmin = (abs (AID_CDELT(aid)) - AID_CDSEARCH(aid)) * (nt - 1)
+	    dwmin = max (1.0D-1, dwmin / dwmax)
+	    AID_NB(aid) = nint (1. / dwmin)
+	}
+
+        # Find the peaks in the target spectrum.
+	if (ID_FTYPE(idt) == ABSORPTION) {
+	    call anegr (IMDATA(idt,1), IMDATA(idt,1), nt)
+	    n = id_upeaks (idt, IMDATA(idt,1), Memr[x], nt, INDEF,
+		int (ID_MINSEP(idt)), 0, AID_NTMAX(aid), 5, INDEF, false)
+	    call anegr (IMDATA(idt,1), IMDATA(idt,1), nt)
+	} else {
+	    n = id_upeaks (idt, IMDATA(idt,1), Memr[x], nt, INDEF,
+		int (ID_MINSEP(idt)), 0, AID_NTMAX(aid), 5, INDEF, false)
+	}
+	call salloc (y, n, TY_REAL)
+	do i = 1, n
+	    Memr[y+i-1] = -IMDATA(idt,nint(Memr[x+i-1]))
+	call xt_sort2 (Memr[y], Memr[x], n)
+
+        # Center and sort the lines.
+        if (AID_XTF(aid) == NULL)
+            call malloc (AID_XTF(aid), n, TY_DOUBLE)
+        else
+            call realloc (AID_XTF(aid), n, TY_DOUBLE)
+        xt = AID_XTF(aid)
+
+        j = 0
+        do i = 1, n {
+	    pix = Memr[x+i-1]
+            pix = id_center (idt, pix, ID_FWIDTH(idt), ID_FTYPE(idt))
+            if (IS_INDEFD(pix))
+                next
+	    if (IS_INDEFD(pix))
+		next
+	    do l = 1, j {
+	        if (abs (pix-Memd[xt+l-1]) < 1.)
+		    break
+	    }
+	    if (l <= j)
+		next
+	    Memd[xt+j] = pix
+            j = j + 1
+        }
+        AID_NTF(aid) = j
+
+        # Sort the lines.
+        if (AID_XT(aid) == NULL)
+            call malloc (AID_XT(aid), j, TY_DOUBLE)
+        else
+            call realloc (AID_XT(aid), j, TY_DOUBLE)
+        xt = AID_XT(aid)
+        if (j > 0)
+            call asrtd (Memd[AID_XTF(aid)], Memd[xt], j)
+        else {
+            call salloc (x, SZ_LINE, TY_CHAR)
+            call sprintf (Memc[x], SZ_LINE, "No target lines found in `%s'")
+                call pargstr (ID_IMAGE(idt))
+            call error (1, Memc[x])
+        }
+
+	# Linearize the lines.
+        if (AID_XTL(aid) == NULL)
+            call malloc (AID_XTL(aid), j, TY_DOUBLE)
+        else
+            call realloc (AID_XTL(aid), j, TY_DOUBLE)
+	xt = AID_XT(aid)
+	xtl = AID_XTL(aid)
+	do i = 1, j
+	    Memd[xtl+i-1] = Memd[xt+i-1] +
+	         AID_CRQUAD(aid) * (Memd[xt+i-1]-AID_CRPIX(aid))**2
+
+        # Debug t: Print list of target lines.
+	if (stridxs ("t", AID_DEBUG(aid,1)) != 0) {
+	    call eprintf ("# Selected target lines:\n")
+	    call eprintf ("#%10s %11s\n")
+	        call pargstr ("Measured")
+		call pargstr ("Undistorted")
+            do i = 1, j {
+                call eprintf ("%11.6g %11.6g\n")
+                    call pargd (Memd[xt+i-1])
+                    call pargd (Memd[xtl+i-1])
+            }
+	    call eprintf ("\n")
+        }
+
+	call sfree (sp)
+end
+
+
+# AID_REFERENCE -- Set reference lines from spectrum or line list.  
+
+procedure aid_reference (aid, ev, flip)
+
+pointer	aid		#I AID pointer
+pointer	ev		#I EV pointer
+int	flip		#I Flip dispersion?
+
+int	i, j, i1, i2, npts, nr, nt, nll, id_peaks(), stridxs()
+double	w, w0, w1, wp, cdelt, wa, wb
+real	sig, wt, center1d()
+pointer	sp, x, idt, idr, specr, xr, sh, label, ll
+double	shdr_wl(), shdr_lw()
+errchk	id_peaks, center1d
+
+begin
+	call smark (sp)
+
+	idr = AID_IDR(aid)
+	npts = ID_NPTS(idr)
+	sh = ID_SH(idr)
+	specr = AID_SPECR(aid)
+	idt = AID_IDT(aid)
+	nt = ID_NPTS(idt)
+
+	# Set reference parameters.
+	if (sh != NULL) {
+	    w0 = min (W0(sh), W1(sh))
+	    w1 = max (W0(sh), W1(sh))
+	    wp = abs (WP(sh))
+	} else {
+	    ll = ID_LL(idr)
+	    nll = ID_NLL(idr)
+	    if (ll == NULL) {
+		ll = ID_LL(idt)
+		nll = ID_NLL(idt)
+	    }
+	    x = ll
+	    w0 = Memd[x]
+	    w1 = Memd[x+nll-1]
+	    wp = INDEFD
+	}
+
+	# Set limits for reference coordinate and dispersion values.
+	AID_CRMIN(aid) = -MAX_DOUBLE
+	AID_CRMAX(aid) = MAX_DOUBLE
+	AID_CDMIN(aid) = 0D0
+	AID_CDMAX(aid) = MAX_DOUBLE
+	if (IS_INDEFD(AID_CDELT(aid))) {
+	    switch (AID_CDDIR(aid)) {
+	    case CDINC:
+		AID_CDSIGN(aid) = 1
+	    case CDDEC:
+		AID_CDSIGN(aid) = -1
+	    default:
+		if (flip == YES)
+		    AID_CDSIGN(aid) = -1
+		else
+		    AID_CDSIGN(aid) = 1
+	    }
+
+	    if (!IS_INDEFD(AID_CRVAL(aid))) {
+		AID_CRMIN(aid) = AID_CRVAL(aid) - AID_CRSEARCH(aid)
+		AID_CRMAX(aid) = AID_CRVAL(aid) + AID_CRSEARCH(aid)
+	    }
+
+	    if (sh != NULL) {
+		i1 = 1
+		i2 = npts
+		sig = 0.
+	    } else {
+		wa = -MAX_DOUBLE
+		wb = MAX_DOUBLE
+	    }
+
+	    AID_W1(aid) = INDEF
+	    AID_W2(aid) = INDEF
+	} else if (IS_INDEFD(AID_CRVAL(aid))) {
+	    switch (AID_CDDIR(aid)) {
+	    case CDINC:
+		cdelt = abs (AID_CDELT(aid))
+		AID_CDSIGN(aid) = 1
+	    case CDDEC:
+		cdelt = -abs (AID_CDELT(aid))
+		AID_CDSIGN(aid) = -1
+	    default:
+		if (flip == YES)
+		    cdelt = -AID_CDELT(aid)
+		else
+		    cdelt = AID_CDELT(aid)
+		if (cdelt < 0.)
+		    AID_CDSIGN(aid) = -1
+		else
+		    AID_CDSIGN(aid) = 1
+	    }
+
+	    AID_CDMIN(aid) = abs (cdelt) - AID_CDSEARCH(aid)
+	    AID_CDMAX(aid) = abs (cdelt) + AID_CDSEARCH(aid)
+
+	    if (sh != NULL) {
+		i1 = 1
+		i2 = npts
+		sig = abs (AID_CDELT(aid)) / wp
+	    } else {
+		wa = -MAX_DOUBLE
+		wb = MAX_DOUBLE
+	    }
+
+	    AID_W1(aid) = INDEF
+	    AID_W2(aid) = INDEF
+	} else {
+	    switch (AID_CDDIR(aid)) {
+	    case CDINC:
+		cdelt = abs (AID_CDELT(aid))
+		AID_CDSIGN(aid) = 1
+	    case CDDEC:
+		cdelt = -abs (AID_CDELT(aid))
+		AID_CDSIGN(aid) = -1
+	    default:
+		if (flip == YES)
+		    cdelt = -AID_CDELT(aid)
+		else
+		    cdelt = AID_CDELT(aid)
+		if (cdelt < 0.)
+		    AID_CDSIGN(aid) = -1
+		else
+		    AID_CDSIGN(aid) = 1
+	    }
+
+	    AID_CRMIN(aid) = AID_CRVAL(aid) - AID_CRSEARCH(aid)
+	    AID_CRMAX(aid) = AID_CRVAL(aid) + AID_CRSEARCH(aid)
+	    AID_CDMIN(aid) = abs (cdelt) - AID_CDSEARCH(aid)
+	    AID_CDMAX(aid) = abs (cdelt) + AID_CDSEARCH(aid)
+
+	    if (cdelt > 0.) {
+		wa = AID_CRMIN(aid) + (cdelt + AID_CDSEARCH(aid)) *
+		    (1 - AID_CRPIX(aid))
+		wb = AID_CRMAX(aid) + (cdelt + AID_CDSEARCH(aid)) *
+		    (nt - AID_CRPIX(aid))
+	    } else {
+		wa = AID_CRMIN(aid) + (cdelt - AID_CDSEARCH(aid)) *
+		    (nt - AID_CRPIX(aid))
+		wb = AID_CRMAX(aid) + (cdelt - AID_CDSEARCH(aid)) *
+		    (1 - AID_CRPIX(aid))
+	    }
+
+	    if (stridxs ("m", AID_DEBUG(aid,1)) != 0) {
+		call eprintf ("wa=%g wb=%g\n")
+		    call pargd (wa)
+		    call pargd (wb)
+	    }
+
+	    if (sh != NULL) {
+		i1 = max (1, min (npts, nint (shdr_wl (sh, wa))))
+		i2 = max (1, min (npts, nint (shdr_wl (sh, wb))))
+		sig = abs (AID_CDELT(aid)) / wp
+	    }
+
+	    AID_W1(aid) = AID_CRVAL(aid) + (1-AID_CRPIX(aid)) * cdelt
+	    AID_W2(aid) = AID_CRVAL(aid) + (nt-AID_CRPIX(aid)) * cdelt
+	}
+
+	# Select lines from line list.
+	if (ID_IMAGE(idr) == EOS) {
+	    ll = ID_LL(idr)
+	    if (ll == NULL)
+		ll = ID_LL(idt)
+	    x = ll
+	    npts = 0
+	    while (!IS_INDEFD(Memd[x])) {
+		if (Memd[x] > wb)
+		    break
+		if (Memd[x] >= wa)
+		    npts = npts + 1
+		x = x + 1
+	    }
+	    x = x - npts
+	    if (npts == 0) {
+		call salloc (x, SZ_LINE, TY_CHAR)
+		call sprintf (Memc[x], SZ_LINE, "No reference lines found")
+		call error (1, Memc[x])
+	    }
+
+	    wa = Memd[x]
+	    wb = Memd[x+npts-1] - Memd[x]
+	    wb = wb / ((AID_BIN1(ev) + 1) / 2)
+	    wa = wa + wb / 2 * (AID_BIN2(ev) - 1)
+	    wb = wa + wb
+
+	    x = ll
+	    npts = 0
+	    while (!IS_INDEFD(Memd[x])) {
+		if (Memd[x] > wb)
+		    break
+		if (Memd[x] >= wa)
+		    npts = npts + 1
+		x = x + 1
+	    }
+	    x = x - npts
+	    if (npts == 0) {
+		call salloc (x, SZ_LINE, TY_CHAR)
+		call sprintf (Memc[x], SZ_LINE, "No reference lines found")
+		call error (1, Memc[x])
+	    }
+
+	    AID_NRMAX(aid) = npts
+	    nr = AID_NRMAX(aid)
+	    AID_NR(aid) = nr
+	    if (AID_XR(aid) == NULL)
+		call malloc (AID_XR(aid), nr, TY_DOUBLE)
+	    else
+		call realloc (AID_XR(aid), nr, TY_DOUBLE)
+	    xr = AID_XR(aid)
+
+	    if (nr < npts) {
+		w = real (npts) / nr
+		do i = 0, nr {
+		    j = i * w
+		    Memd[xr+i] = Memd[x+j]
+		}
+	    } else {
+		do i = 0, nr-1
+		    Memd[xr+i] = Memd[x+i]
+	    }
+
+	# Select lines using reference spectrum.
+	} else {
+	    wb = (i2 - i1) / ((AID_BIN1(ev) + 1) / 2)
+	    i1 = max (i1, nint (i1 + wb / 2 * (AID_BIN2(ev) - 1)))
+	    i2 = min (i2, nint (i1 + wb))
+
+	    if (i2 - i1 + 1 < 100) {
+		i1 = 1
+		i2 = npts
+	    }
+	    npts = i2 - i1 + 1
+
+	    if (specr == NULL)
+		call malloc (specr, npts, TY_REAL)
+	    else
+		call realloc (specr, npts, TY_REAL)
+	    AID_SPECR(aid) = specr
+	    AID_X1R(aid) = i1
+	    AID_X2R(aid) = i2
+	    wa = Memr[SX(sh)+i1-1]
+	    wb = Memr[SX(sh)+i2-1]
+	    call amovr (IMDATA(idr,i1), Memr[specr], npts)
+
+	    if (sig > 1.) {
+		ID_MINSEP(idr) = sig * ID_MINSEP(idt)
+		ID_FWIDTH(idr) = sig * ID_FWIDTH(idt)
+		sig = sig / 1.1774
+		j = nint (3 * sig)
+		call malloc (x, npts, TY_REAL)
+		call malloc (xr, npts+2*j+1, TY_REAL)
+		xr = xr + j
+		call amovr (Memr[specr], Memr[xr], npts)
+		do i = 1, j {
+		    wt = exp (-0.5 * (i / sig) ** 2)
+		    call amulkr (Memr[specr], wt, Memr[x], npts)
+		    call aaddr (Memr[x], Memr[xr+i], Memr[xr+i], npts)
+		    call aaddr (Memr[x], Memr[xr-i], Memr[xr-i], npts)
+		}
+		call amovr (Memr[xr], Memr[specr], npts)
+		call mfree (x, TY_REAL)
+		call mfree (xr-j, TY_REAL)
+	    }
+
+	    call salloc (x, npts, TY_REAL)
+
+	    # Find the peaks in the reference spectrum.
+	    AID_NRMAX(aid) = 2 * AID_NTF(aid)
+	    if (ID_FTYPE(idr) == ABSORPTION) {
+		call anegr (Memr[specr], Memr[specr], nt)
+		nr = id_peaks (idr, Memr[specr], Memr[x], npts, INDEF,
+		    int (ID_MINSEP(idr)), 0, AID_NRMAX(aid), INDEF, false)
+		call anegr (Memr[specr], Memr[specr], nt)
+	    } else {
+		nr = id_peaks (idr, Memr[specr], Memr[x], npts, INDEF,
+		    int (ID_MINSEP(idr)), 0, AID_NRMAX(aid), INDEF, false)
+	    }
+
+	    # Center and sort the lines.
+	    if (AID_XR(aid) == NULL)
+		call malloc (AID_XR(aid), nr, TY_DOUBLE)
+	    else
+		call realloc (AID_XR(aid), nr, TY_DOUBLE)
+	    xr = AID_XR(aid)
+
+	    j = 0
+	    label = NULL
+	    do i = 1, nr {
+		wt = center1d (Memr[x+i-1], Memr[specr], npts, ID_FWIDTH(idr),
+		    ID_FTYPE(idr), ID_CRADIUS(idt), 0.)
+		if (IS_INDEF(wt))
+		    next
+		w = shdr_lw (sh, double(wt+i1-1))
+		Memd[xr+j] = w
+		call id_match (idt, w, Memd[xr+j], label, -2.)
+		if (IS_INDEFD(Memd[xr+j]) || (j>0 && Memd[xr+j]==Memd[xr+j-1]))
+		    next
+		j = j + 1
+	    }
+	    call mfree (label, TY_CHAR)
+	    nr = j
+	    AID_NR(aid) = nr
+
+	    # Sort the lines.
+	    if (j > 0)
+		call asrtd (Memd[xr], Memd[xr], nr)
+	    else {
+		call salloc (x, SZ_LINE, TY_CHAR)
+		call sprintf (Memc[x], SZ_LINE,
+		    "No reference lines found in `%s'")
+		    call pargstr (ID_IMAGE(idr))
+		call error (1, Memc[x])
+	    }
+	}
+
+	#AID_NT(aid) = min (2 * AID_NR(aid), AID_NTF(aid))
+	AID_NT(aid) = AID_NTF(aid)
+	call asrtd (Memd[AID_XTF(aid)], Memd[AID_XT(aid)], AID_NT(aid))
+
+	# Debug w: Print wavelength bin limits.
+	if (stridxs ("w", AID_DEBUG(aid,1)) != 0) {
+	    call eprintf ("%2d/%-2d %g %g\n")
+		call pargi (AID_BIN1(ev))
+		call pargi (AID_BIN2(ev))
+		call pargd (wa)
+		call pargd (wb)
+	}
+
+	# Debug b: Print search limits.
+	if (stridxs ("b", AID_DEBUG(aid,1)) != 0) {
+	    if (ev == AID_EVS(aid)) {
+		call eprintf ("Search: CRVAL = %.8g - %.8g,  CDELT = %.5g - %.5g\n\n")
+		    call pargd (AID_CRMIN(aid))
+		    call pargd (AID_CRMAX(aid))
+		    call pargd (AID_CDMIN(aid))
+		    call pargd (AID_CDMAX(aid))
+	    }
+	}
+
+	# Debug r: Print list of reference lines.
+	if (stridxs ("r", AID_DEBUG(aid,1)) != 0) {
+	    call eprintf ("# Selected reference lines:\n")
+	    do i = 1, nr {
+		call eprintf ("%10.6g\n")
+		    call pargd (Memd[xr+i-1])
+	    }
+	    call eprintf ("\n")
+	}
+
+	call sfree (sp)
+end
+
+
+# AID_AUTOID1 -- Automatically identify lines.
+# This routine takes preset target and reference line lists and tries to
+# find correspondences.  It returns lists of possible correspondences
+# and dispersions.
+
+procedure aid_autoid1 (aid, ev)
+
+pointer	aid		#I AID pointer
+pointer	ev		#I EV pointer
+
+int	i, nn, n1, n2, nr1, nr2, n, nd
+pointer	sp, idt, x1, x2, x3, r1, s1, r2, s2, votes, svotes
+pointer	x, y, w, w1, dw, nw, nv
+
+int	aid_rsort(), aid_vsort(), stridxs()
+extern	aid_rsort, aid_vsort
+errchk	aid_select, aid_disp
+
+begin
+	call smark (sp)
+
+	idt = AID_IDT(aid)
+	nn = AID_NN(aid)
+	x1 = AID_XR(aid)
+	n1 = AID_NR(aid)
+	x2 = AID_XTL(aid)
+	x3 = AID_XT(aid)
+	n2 = AID_NT(aid)
+
+        # Debug l: Graph lines and spectra.
+	if (stridxs ("l", AID_DEBUG(aid,1)) != 0)
+	    call aid_lgraph (aid, Memd[x1], n1, Memd[x2], n2)
+
+	# Make ratio lists.
+	i = min (nn, n1-1)
+	nr1 = (n1-i) * i * (i - 1) / 2 + i * (i - 1) * (i - 2) / 6
+	call salloc (r1, nr1, TY_REAL)
+	call aid_ratios (aid, Memd[x1], n1, 1, Memr[r1], nr1, 1)
+	call salloc (s1, nr1, TY_INT)
+	do i = 1, nr1
+	    Memi[s1+i-1] = i
+	call gqsort (Memi[s1], nr1, aid_rsort, r1)
+
+	i = min (nn, n2-1)
+	nr2 = (n2-i) * i * (i - 1) / 2 + i * (i - 1) * (i - 2) / 6
+	call salloc (r2, 2*nr2, TY_REAL)
+	call aid_ratios (aid, Memd[x2], n2, AID_CDSIGN(aid), Memr[r2], nr2, 2)
+	call salloc (s2, nr2, TY_INT)
+	do i = 1, nr2
+	    Memi[s2+i-1] = i
+	call gqsort (Memi[s2], nr2, aid_rsort, r2)
+
+	call salloc (votes, n1 * n2, TY_INT)
+	call aid_votes (aid, Memr[r1], Memi[s1], nr1, Memr[r2], Memi[s2],
+	    nr2, Memd[x1], Memd[x2], Memi[votes], n1, n2)
+
+	call salloc (svotes, n1 * n2, TY_INT)
+	do i = 1, n1 * n2
+	    Memi[svotes+i-1] = i
+	call gqsort (Memi[svotes], n1*n2, aid_vsort, votes)
+
+	do n = 1, n1 * n2
+	    if (Memi[votes+Memi[svotes+n-1]-1] < 1)
+	        break
+	n = max (3 * n2, n-1)
+
+	call malloc (x, n, TY_REAL)
+	call malloc (y, n, TY_REAL)
+	call salloc (w, n, TY_REAL)
+	iferr (call aid_select (aid, Memd[x1], Memd[x2], Memd[x3], Memi[votes],
+	    Memi[svotes], n1, n2, Memr[x], Memr[y], Memr[w], n)) {
+	    call sfree (sp)
+	    call erract (EA_ERROR)
+	}
+
+	nd = AID_NDMAX(aid)
+	call malloc (w1, nd, TY_REAL)
+	call malloc (dw, nd, TY_REAL)
+	call salloc (nw, nd, TY_INT)
+	call salloc (nv, nd, TY_INT)
+	call aid_disp (aid, Memr[y], Memr[x], Memr[w], n, Memr[w1], Memr[dw],
+	    Memi[nw], Memi[nv], nd)
+
+	AID_X(ev) = x
+	AID_Y(ev) = y
+	AID_N(ev) = n
+	AID_A(ev) = w1
+	AID_B(ev) = dw
+	AID_ND(ev) = nd
+
+	call sfree (sp)
+end
+
+
+# AID_RATIOS -- Generate list of spacing ratios from list of lines.
+
+procedure aid_ratios (aid, x, n, cdsign, r, nr, nv)
+
+pointer	aid		#I AID pointer
+double	x[n]		#I Line positions (sorted)
+int	n		#I Number of lines
+int	cdsign		#I Sign of dispersion
+real	r[nr,nv]	#O Ratios
+int	nr		#O Number of ratios
+int	nv		#I Number of values
+
+int	i, j, k, l, nn, stridxs()
+real	minr, maxr, xi, xj, xk, xij, xjk, sig, ratio, err
+
+begin
+	nn = AID_NN(aid)
+	sig = AID_SIG(aid)
+	minr = AID_MINRATIO(aid)
+	maxr = 1 / AID_MINRATIO(aid)
+
+	# Compute ratios.
+	l = 0
+	if (cdsign == 1) {
+	    do i = 1, n-2 {
+		xi = x[i]
+		do j = i+1, min (i+nn-1, n-1) {
+		    xj = x[j]
+		    xij = xj - xi
+		    do k = j+1, min (i+nn, n) {
+			xk = x[k]
+			xjk = xk - xj
+			ratio = xij / xjk
+
+			l = l + 1
+			if (nv == 1) {
+			    if (ratio < minr || ratio > maxr)
+				r[l,1] = 1000.
+			    else
+				r[l,1] = ratio
+			} else {
+			    if (ratio < minr || ratio > maxr) {
+				r[l,1] = 1000.
+				r[l,2] = 1000.
+			    } else {
+				err = sig * sqrt (2*(1+ratio+ratio**2)) / xjk
+				r[l,1] = ratio - err
+				r[l,2] = ratio + err
+			    }
+			}
+		    }
+		}
+	    }
+	} else {
+	    do i = n, 3, -1 {
+		xi = x[i]
+		do j = i-1, max (i-nn+1, 2), -1 {
+		    xj = x[j]
+		    xij = xi - xj
+		    do k = j-1, max (i-nn, 1), -1 {
+			xk = x[k]
+			xjk = xj - xk
+			ratio = xij / xjk
+
+			l = l + 1
+			if (nv == 1) {
+			    if (ratio < minr || ratio > maxr)
+				r[l,1] = 1000.
+			    else
+				r[l,1] = ratio
+			} else {
+			    if (ratio < minr || ratio > maxr) {
+				r[l,1] = 1000.
+				r[l,2] = 1000.
+			    } else {
+				err = sig * sqrt (2*(1+ratio+ratio**2)) / xjk
+				r[l,1] = ratio - err
+				r[l,2] = ratio + err
+			    }
+			}
+		    }
+		}
+	    }
+	}
+	nr = l
+
+	# Debug c: Print list of line ratios.
+	if (stridxs ("c", AID_DEBUG(aid,1)) != 0) {
+	    do l = 1, nr {
+		call aid_lines (l, n, nn, i, j, k)
+		if (nv == 1)
+		    call printf ("%2d %2d %2d %8.2f %8.2f %8.2f %6.4f\n")
+		else
+		    call printf ("%2d %2d %2d %8.2f %8.2f %8.2f %6.4f %6.4f\n")
+		    call pargi (i)
+		    call pargi (j)
+		    call pargi (k)
+		    if (cdsign == 1) {
+			call pargd (x[i])
+			call pargd (x[j])
+			call pargd (x[k])
+		    } else {
+			call pargd (x[n-i+1])
+			call pargd (x[n-j+1])
+			call pargd (x[n-k+1])
+		    }
+		    call pargr (r[l,1])
+		    if (nv == 2)
+			call pargr (r[l,2])
+	    }
+	}
+end
+
+
+# AID_LINES --  Convert ratio index to line indices.
+
+procedure aid_lines (s, n, nn, i, j, k)
+
+int	s		# Index into ratio array
+int	n		# Number of lines
+int	nn		# Number of neigbhors
+int	i		#O Index of first line
+int	j		#O Index of second line
+int	k		#O Index of third line
+
+int	l
+
+begin
+	k = s
+	for (i=1;;i=i+1) {
+	    l = min (nn, n-i)
+	    l = l * (l-1) / 2
+	    if (k <= l)
+		break
+	    k = k - l
+	}
+	for (j=i+1;;j=j+1) {
+	    l = min (nn-1, n-j)
+	    if (k <= l)
+		break
+	    k = k - l
+	}
+	k = k + j
+end
+
+
+# AID_RSORT -- Compare ratio array with smallest first.
+
+int procedure aid_rsort (ptr, i, j)
+
+pointer	ptr		#I Pointer to array to be sorted.
+int	i		#I Index 1
+int	j		#I Index 2
+
+real	a, b
+
+begin
+	a = Memr[ptr+i-1]
+	b = Memr[ptr+j-1]
+
+	if (a < b)
+	    return (-1)
+	else if (b < a)
+	    return (1)
+	else
+	    return (0)
+end
+
+
+# AID_VSORT -- Compare vote array with biggest first.
+
+int procedure aid_vsort (ptr, i, j)
+
+pointer	ptr		#I Pointer to array to be sorted.
+int	i		#I Index 1
+int	j		#I Index 2
+
+int	a, b
+
+begin
+	a = Memi[ptr+i-1]
+	b = Memi[ptr+j-1]
+
+	if (a < b)
+	    return (1)
+	else if (b < a)
+	    return (-1)
+	else
+	    return (0)
+end
+
+
+# AID_VOTES -- Find ratio matches and increment the vote array.
+
+procedure aid_votes (aid, r1, s1, nr1, r2, s2, nr2, x1, x2, votes, n1, n2)
+
+pointer	aid		#I AID pointer
+real	r1[nr1]		#I Ratio array (reference)
+int	s1[nr1]		#I Sort array
+int	nr1		#I Number of ratios
+real	r2[nr2,2]	#I Ratio array (target)
+int	s2[nr2]		#I Sort array
+int	nr2		#I Number of ratios
+double	x1[n1]		#I Reference lines
+double	x2[n2]		#I Target lines
+int	votes[n1,n2]	#O Votes
+int	n1, n2		#I Size of votes array
+
+int	i, j, nn, np, start, stridxs()
+real	maxr, ra, rb1, rb2
+pointer	sp, a, b
+
+begin
+	nn = AID_NN(aid)
+	np = max (3, min (AID_NP(aid), n1 - 5))
+	maxr = 1. / AID_MINRATIO(aid)
+
+	call smark (sp)
+	call salloc (a, np, TY_INT)
+	call salloc (b, np, TY_INT)
+
+	call aclri (votes, n1*n2)
+
+	start = 1
+	do j = 1, nr2 {
+	    rb1 = r2[s2[j],1]
+	    if (rb1 > maxr)
+		break
+	    rb2 = r2[s2[j],2]
+	    do i = start, nr1 {
+		ra = r1[s1[i]]
+		if (ra > rb2)
+		    break
+		if (ra < rb1) {
+		    start = i + 1
+		    next
+		}
+		call aid_lines (s1[i], n1, nn, Memi[a], Memi[a+1], Memi[a+2])
+		call aid_lines (s2[j], n2, nn, Memi[b], Memi[b+1], Memi[b+2])
+		call aid_addlines (aid, r1, nr1, s1[i], r2, nr2, s2[j], nn,
+		    Memi[a], Memi[b], np, votes, n1, n2)
+	    }
+	}
+
+        # Debug v: Print vote array.
+	if (stridxs ("v", AID_DEBUG(aid,1)) != 0) {
+	    call printf ("%4w")
+	    do i = 1, n2 {
+		call printf (" %3d")
+		    call pargi (nint (x2[i]))
+	    }
+	    call printf ("\n")
+	    do i = 1, n1 {
+		call printf ("%4d")
+		    call pargi (nint (x1[i]))
+		do j = 1, n2 {
+		    call printf (" %3d")
+			call pargi (votes[i,j])
+		}
+		call printf ("\n")
+	    }
+	    call printf ("\n")
+	    call flush (STDOUT)
+	}
+
+	call sfree (sp)
+end
+
+
+# AID_ADDLINES -- Starting with a matching triplets add more lines.
+# The lines are added recursively.  To avoid recursive calls this
+# routine is repeated to a maximum depth.  The indentation is intentionally
+# non-standard.
+
+procedure aid_addlines (aid, r1, nr1, s1, r2, nr2, s2, nn, a, b, npattern,
+	votes, n1, n2)
+
+pointer	aid		#I AID pointer
+real	r1[nr1]		#I Reference ratios
+int	nr1		#I Number of ratios
+int	s1		#I Ratio index
+real	r2[nr2,2]	#I Target ratios
+int	nr2		#I Number of ratios
+int	s2		#I Ratio index
+int	nn		#I Number of neighbors
+int	a[npattern]	#I Reference lines (indices)
+int	b[npattern]	#I Target lines (indices)
+int	npattern	#I Number of lines in pattern
+int	votes[n1,n2]	#O Vote array
+int	n1, n2		#O Number of reference and target lines
+
+int	i, j, i1, j1, na, nb
+
+begin
+	na = min (a[1] + nn, n1)
+	nb = min (b[1] + nn, n2)
+	i1 = s1 - a[3]
+	j1 = s2 - b[3]
+
+	if (npattern > 3) {
+	    for (a[4]=a[3]+1; a[4]<=na; a[4]=a[4]+1) {
+		for (b[4]=b[3]+1; b[4]<=nb; b[4]=b[4]+1) {
+		    i = i1 + a[4]
+		    j = j1 + b[4]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 4) {
+	    for (a[5]=a[4]+1; a[5]<=na; a[5]=a[5]+1) {
+		for (b[5]=b[4]+1; b[5]<=nb; b[5]=b[5]+1) {
+		    i = i1 + a[5]
+		    j = j1 + b[5]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 5) {
+	    for (a[6]=a[5]+1; a[6]<=na; a[6]=a[6]+1) {
+		for (b[6]=b[5]+1; b[6]<=nb; b[6]=b[6]+1) {
+		    i = i1 + a[6]
+		    j = j1 + b[6]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 6) {
+	    for (a[7]=a[6]+1; a[7]<=na; a[7]=a[7]+1) {
+		for (b[7]=b[6]+1; b[7]<=nb; b[7]=b[7]+1) {
+		    i = i1 + a[7]
+		    j = j1 + b[7]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 7) {
+	    for (a[8]=a[7]+1; a[8]<=na; a[8]=a[8]+1) {
+		for (b[8]=b[7]+1; b[8]<=nb; b[8]=b[8]+1) {
+		    i = i1 + a[8]
+		    j = j1 + b[8]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 8) {
+	    for (a[9]=a[8]+1; a[9]<=na; a[9]=a[9]+1) {
+		for (b[9]=b[8]+1; b[9]<=nb; b[9]=b[9]+1) {
+		    i = i1 + a[9]
+		    j = j1 + b[9]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+	if (npattern > 9) {
+	    for (a[10]=a[9]+1; a[10]<=na; a[10]=a[10]+1) {
+		for (b[10]=b[9]+1; b[10]<=nb; b[10]=b[10]+1) {
+		    i = i1 + a[10]
+		    j = j1 + b[10]
+		    if (r1[i] < r2[j,1] || r1[i] > r2[j,2])
+			next
+		    call aid_vote (aid, a, b, 10, votes, n1, n2)
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+	}
+	}
+	} else {
+	    call aid_vote (aid, a, b, npattern, votes, n1, n2)
+	}
+end
+
+
+# AID_VOTE -- Add votes for the lines in the pattern to the vote array.
+
+procedure aid_vote (aid, a, b, npattern, votes, n1, n2)
+
+pointer	aid		#I AID pointer
+int	a[npattern]	#I Reference lines (indices)
+int	b[npattern]	#I Target lines (indices)
+int	npattern	#I Number of lines in pattern
+int	votes[n1,n2]	#O Vote array
+int	n1, n2		#O Number of reference and target lines
+
+int	i, stridxs()
+pointer	xr, xt
+
+begin
+	if (AID_CDSIGN(aid) == 1) {
+	    do i = 1, npattern
+		votes[a[i],b[i]] = votes[a[i],b[i]] + 1
+	} else {
+	    do i = 1, npattern
+		votes[a[i],n2-b[i]+1] = votes[a[i],n2-b[i]+1] + 1
+	}
+
+        # Debug a: Print line assignments.
+	if (stridxs ("a", AID_DEBUG(aid,1)) != 0) {
+	    xr = AID_XR(aid)-1
+	    xt = AID_XT(aid)-1
+	    if (AID_CDSIGN(aid) == 1) {
+		do i = 1, npattern {
+		    call eprintf (" %6g %6g %5d")
+			call pargd (Memd[xr+a[i]])
+			call pargd (Memd[xt+b[i]])
+			call pargi (b[i])
+		}
+	    } else {
+		xt = xt+n2+1
+		do i = 1, npattern {
+		    call eprintf (" %6g %6g %5d")
+			call pargd (Memd[xr+a[i]])
+			call pargd (Memd[xt-b[i]])
+			call pargi (n2-b[i]+1)
+		}
+	    }
+	    call eprintf ("\n")
+	}
+end
+
+
+# AID_SELECT -- Select top vote getters.
+
+procedure aid_select (aid, x1, x2, x3, votes, svotes, n1, n2, x, y, w, ns)
+
+pointer	aid		#I AID pointer
+double	x1[n1]		#I Reference lines
+double	x2[n2]		#I Linearized target lines
+double	x3[n2]		#I Target lines
+int	votes[n1,n2]	#I Vote array
+int	svotes[ARB]	#I Sort indices for vote array
+int	n1, n2		#I Number of lines
+real	x[ns]		#O Selected target coordinates
+real	y[ns]		#O Selected reference coordinates
+real	w[ns]		#O Weight (votes)
+int	ns		#U Maximum number on input, number selected on output
+
+int	i, j, k, n
+double	a, b
+bool	check
+
+begin
+	check = (AID_CRMIN(aid) > -MAX_DOUBLE / 10. &&
+	    AID_CRMAX(aid) < MAX_DOUBLE / 10.)
+
+	# Select the highest votes.
+	n = 0
+	for (k=1; k<=n1*n2 && n<ns; k=k+1) {
+	    i = mod (svotes[k]-1, n1) + 1 
+	    j = (svotes[k]-1) / n1 + 1 
+	    if (votes[i,j] < 1)
+		break
+	    if (check) {
+		a = (x2[j] - AID_CRPIX(aid)) * AID_CDSIGN(aid) * AID_CDMIN(aid)
+		b = (x2[j] - AID_CRPIX(aid)) * AID_CDSIGN(aid) * AID_CDMAX(aid)
+		if (x1[i] < AID_CRMIN(aid) + min (a,b))
+		    next
+		if (x1[i] > AID_CRMAX(aid) + max (a,b))
+		    next
+	    }
+	    n = n + 1
+	    x[n] = x3[j]
+	    y[n] = x1[i]
+	    w[n] = votes[i,j]
+	}
+	ns = n
+
+	if (ns < 1)
+	    call error (1, "No matches found")
+end
+
+
+# AID_DISP -- Given a set of candidate identifications (pixel, wavelength)
+# find all linear dispersions between two or more identifications which
+# satisfy the dispersion constraints.  The list of ranked dispersions with
+# higher rankings for higher number of points the dispersion goes through
+# higher total votes for the points.  Hopefully the true dispersion will be
+# in the highest ranked dispersions.
+
+procedure aid_disp (aid, x, y, v, n, w1, dw, nw, nv, nd)
+
+pointer	aid		#I AID pointer
+real	x[n]		#I Array of candidate reference coordinates
+real	y[n]		#I Array of candidate target coordinates
+real	v[n]		#I Votes
+int	n		#I Number of candidate pairs
+real	w1[nd]		#O Dispersion origin
+real	dw[nd]		#O Dispersion slope
+int	nw[nd]		#O Number of points
+int	nv[nd]		#O Sum of votes
+int	nd		#U Number of dispersions
+
+bool	debug, skip
+int	i, j, k, l, m, ii, sumn, sumv, stridxs()
+double	aw, bw, cw, sumx, sumy, sumyy, sumxy
+pointer	iii
+
+begin
+	# Sort the candidates by reference coordinate.
+	call xt_sort2 (x, y, n)
+
+	debug = (stridxs ("m", AID_DEBUG(aid,1)) != 0)
+	if (debug) {
+	    call eprintf ("# Selected pairs with votes.\n")
+	    do i = 1, n {
+		call eprintf ("%4d %8.6g %8.6g %d\n")
+		call pargi (i)
+		call pargr (x[i])
+		call pargr (y[i])
+		call pargr (v[i])
+	    }
+	    call eprintf ("# Dispersions to check up to %d.\n")
+	        call pargi (nd)
+	}
+
+	m = 0
+	ii = 0
+	call malloc (iii, nd, TY_INT)
+	do i = 1, n-2 {
+	    do j = i+1, n-1 {
+		if (x[j] == x[i] || y[j] == y[i])
+		    next
+
+		bw = (x[j] - x[i]) / (y[j] - y[i])
+		aw = x[i] - bw * y[i]
+		cw = aw + bw * AID_CRPIX(aid)
+
+		# Check dispersion ranges.
+		skip = false
+		if (abs (bw) < AID_CDMIN(aid) || abs (bw) > AID_CDMAX(aid))
+		    skip = true
+		else if (cw < AID_CRMIN(aid) || cw > AID_CRMAX(aid))
+		    skip = true
+		if (AID_CDSIGN(aid) * bw < 0.)
+		    skip = true
+		if (skip)
+		    next
+
+		sumn = 2
+		sumv = v[i] + v[j]
+		sumx = x[i] + x[j]
+		sumy = y[i] + y[j]
+		sumyy = y[i]*y[i] + y[j]*y[j]
+		sumxy = x[i]*y[i] + x[j]*y[j]
+
+		do k = j+1, n {
+		    if (abs ((x[k] - aw - bw * y[k]) / bw) > 2.) 
+			next
+
+		    sumn = sumn + 1
+		    sumv = sumv + v[k]
+		    sumx = sumx + x[k]
+		    sumy = sumy + y[k]
+		    sumyy = sumyy + y[k]*y[k]
+		    sumxy = sumxy + x[k]*y[k]
+		}
+
+		aw = (sumx*sumyy - sumy*sumxy) / (sumn * sumyy - sumy * sumy)
+		bw = (sumn*sumxy - sumx*sumy) / (sumn * sumyy - sumy * sumy)
+		cw = aw + bw * AID_CRPIX(aid)
+		ii = ii + 1
+
+		if (debug) {
+		    call eprintf ("     %4d %4d %4d %8.5g  %8.3g  %8d %8d")
+		    call pargi (ii)
+		    call pargi (i)
+		    call pargi (j)
+		    call pargd (aw+bw*(ID_NPTS(AID_IDT(aid))/2.+1))
+		    call pargd (bw)
+		    call pargi (sumn)
+		    call pargi (sumv)
+		}
+
+		# Check if already found.
+		for (k = 1; k <= m; k = k + 1)
+		    if (abs ((x[1]-aw)/bw - (x[1]-w1[k])/dw[k]) < 2. &&
+		        abs ((x[n]-aw)/bw - (x[n]-w1[k])/dw[k]) < 2.)
+			break
+		if (k <= m) {
+		    if (sumn > nw[k] || (sumn == nw[k] && sumv > nv[k])) {
+			for (l = k; l > 1; l = l - 1) {
+			    if (sumn<nw[l-1] || (sumn==nw[l-1] && sumv<nv[l-1]))
+				break
+			    w1[l] = w1[l-1]
+			    dw[l] = dw[l-1]
+			    nw[l] = nw[l-1]
+			    nv[l] = nv[l-1]
+			    Memi[iii+l-1] = Memi[iii+l-2]
+			}
+			if (debug) {
+			    call eprintf (" replace %4d\n")
+			    call pargi (Memi[iii+l-1])
+			}
+			w1[l] = aw
+			dw[l] = bw
+			nw[l] = sumn
+			nv[l] = sumv
+			Memi[iii+l-1] = ii
+		    } else if (debug) {
+			call eprintf (" use %4d\n")
+			call pargi (Memi[iii+k-1])
+		    }
+		    next
+		}
+
+		# Check dispersion ranges.
+		if (abs (bw) < AID_CDMIN(aid) || abs (bw) > AID_CDMAX(aid))
+		    skip = true
+		else if (cw < AID_CRMIN(aid) || cw > AID_CRMAX(aid))
+		    skip = true
+		if (AID_CDSIGN(aid) * bw < 0.)
+		    skip = true
+		if (skip) {
+		    if (debug)
+			call eprintf (" out of range\n")
+		    next
+		}
+
+		# Add to ordered list.
+		for (k = 1; k <= m; k = k + 1)
+		    if (sumn > nw[k] || (sumn == nw[k] && sumv > nv[k]))
+			break
+		if (k <= nd) {
+		    if (m < nd) {
+			m = m + 1
+			if (debug)
+			    call eprintf (" add\n")
+		    } else if (debug) {
+			call eprintf (" bump %4d\n")
+			    call pargi (Memi[iii+m-1])
+		    }
+		    for (l = m; l > k; l = l - 1) {
+			w1[l] = w1[l-1]
+			dw[l] = dw[l-1]
+			nw[l] = nw[l-1]
+			nv[l] = nv[l-1]
+			Memi[iii+l-1] = Memi[iii+l-2]
+		    }
+		    w1[k] = aw
+		    dw[k] = bw
+		    nw[k] = sumn
+		    nv[k] = sumv
+		    Memi[iii+k-1] = ii
+		} else if (debug)
+		    call eprintf (" failed\n")
+	    }
+	}
+
+	nd = m
+
+	if (debug) {
+	    call eprintf ("# Final ordered dispersions to try.\n")
+	    do i = 1, nd {
+		call eprintf ("     %4d %8.5g  %8.3g  %8d %8d\n")
+		call pargi (Memi[iii+i-1])
+		call pargr (w1[i]+dw[i]*(ID_NPTS(AID_IDT(aid))/2.+1))
+		call pargr (dw[i])
+		call pargi (nw[i])
+		call pargi (nv[i])
+	    }
+	}
+	call mfree (iii, TY_INT)
+
+        # Debug d: Graph dispersions.
+	if (stridxs ("d", AID_DEBUG(aid,1)) != 0)
+	    call aid_dgraph (aid, x, y, n, w1, dw, nd)
+end
+
+
+# AID_EVAL -- Evaluate possible solutions.
+
+double procedure aid_eval (aid, ev, nd)
+
+pointer	aid		#I AID pointer
+pointer	ev		#I EV pointer
+int	nd		#I Dispersion candidate to evaluate
+double	best		#O Best statistic
+
+int	i, n
+pointer	idt, x, y
+double	a, b, c, d, rms, fmatch, ftmatch
+int	stridxs()
+
+int	ncandidate, nmatch1, nmatch2
+common	/llstat/ ncandidate, nmatch1, nmatch2
+
+define	done_	90
+
+begin
+	best = INDEFD
+	if (nd > AID_ND(ev))
+	    return (best)
+
+	idt = AID_IDT(aid)
+	x = AID_X(ev) - 1
+	y = AID_Y(ev) - 1
+	n = AID_N(ev)
+
+	a = Memr[AID_A(ev)+nd-1]
+	b = Memr[AID_B(ev)+nd-1]
+	c = ID_NPTS(AID_IDT(aid)) / 2. + 1
+	if (IS_INDEFD(AID_CDELT(aid)))
+	    d = b
+	else
+	    d = AID_CDELT(aid)
+
+	ID_IC(idt) = AID_IC1(aid)
+	ID_NFEATURES(idt) = 0
+	do i = 1, n {
+	    if (abs ((Memr[y+i] - a - b * Memr[x+i]) / b) < 2.)
+		call id_newfeature (idt, double(Memr[x+i]), double(Memr[x+i]),
+		    double(Memr[y+i]), 1.0D0, ID_FWIDTH(idt), ID_FTYPE(idt),
+		    NULL)
+	}
+	if (ID_NFEATURES(idt) <= 1)
+	    goto done_
+
+	call dcvfree (ID_CV(idt))
+	iferr (call aid_dofit (aid, idt, d, rms, fmatch, ftmatch, best))
+	    goto done_
+
+	# Debug s: Print search iterations.
+	if (stridxs ("s", AID_DEBUG(aid,1)) != 0) {
+	    call eprintf (
+	"%2d/%-2d %8.2f %8.3f %3d %3d/%-3d %3d/%-3d %3d %3d %6.3f %5.2f\n")
+		call pargi (AID_BIN1(ev))
+		call pargi (AID_BIN2(ev))
+		call pargd (a+c*b)
+		call pargd (b)
+		call pargi (ID_NFEATURES(idt))
+		call pargi (nmatch2)
+		call pargi (ncandidate)
+		call pargi (nint(min (ncandidate, AID_NT(aid))*(1-ftmatch)))
+		call pargi (min (ncandidate, AID_NT(aid)))
+		call pargi (nint(100.*fmatch))
+		call pargi (nint(100.*ftmatch))
+		call pargd (rms)
+		call pargd (best)
+	}
+
+	if (best < AID_BEST(aid)) {
+	    AID_FMATCH(aid) = fmatch
+	    AID_FTMATCH(aid) = ftmatch
+	    AID_RMS(aid) = rms
+	    AID_BEST(aid) = best
+	    ID_IC(idt) = AID_IC2(aid)
+	    call id_saveid (idt, "autoidentify")
+	}
+
+done_
+	ID_IC(idt) = AID_IC2(aid)
+	return (best)
+end
+
+
+# AID_DOFIT -- From a set of candidate identifications fit and evaluate
+# a dispersion solution.
+
+procedure aid_dofit (aid, id, cdelt, rms, fmatch, ftmatch, best)
+
+pointer	aid		#I AID pointer
+pointer	id		#I ID pointer
+double	cdelt		#I Dispersion to use in pixel rms conversion
+double	rms		#O Final RMS in pixels
+double	fmatch		#O Line list non-matching fraction
+double	ftmatch		#O Target line non-matching fraction
+double	best		#O Best fit parameter
+
+int	i, j, k, l, nmin, nfound, nt, ntmatch, maxfeatures, stridxs()
+double	fit, user, dcveval(), id_fitpt()
+pointer	cv, xt, label
+
+int	ncandidate, nmatch1, nmatch2
+common	/llstat/ ncandidate, nmatch1, nmatch2
+
+errchk	id_dofit, id_fitdata, id_fitfeatures, id_linelist, id_match
+
+begin
+	maxfeatures = ID_MAXFEATURES(id)
+	ID_MAXFEATURES(id) = 1000
+	iferr {
+	    do k = 1, 3 {
+		if (ID_NFEATURES(id) < 2)
+		    call error (0, "aid_dofit: not enough features")
+		if (k > 1)
+		    call id_linelist (id)
+
+		if (stridxs ("i", AID_DEBUG(aid,1)) != 0)
+		    call id_dofit (id, YES)
+		else
+		    call id_dofit (id, NO)
+		do l = AID_ORD(aid)-1, 2, -1 {
+		    cv = ID_CV(id)
+		    user = dcveval (cv, 1D0)
+		    fit = (dcveval (cv, double (ID_NPTS(id)/2)) - user) /
+			(dcveval (cv, double (ID_NPTS(id))) - user)
+		    if (abs (fit - 0.5) <= AID_MAXNL(aid))
+		        break
+		    if (stridxs ("n", AID_DEBUG(aid,1)) != 0) {
+			call eprintf (
+			    "order %d: non-linearity of %.1f%% > %.1f%%\n")
+			call pargi (l+1)
+			call pargd (100*abs(fit-0.5))
+			call pargr (100*AID_MAXNL(aid))
+		    }
+		    call ic_puti (ID_IC(id), "order", l)
+		    if (stridxs ("i", AID_DEBUG(aid,1)) != 0)
+			call id_dofit (id, YES)
+		    else
+			call id_dofit (id, NO)
+		    call ic_puti (ID_IC(id), "order", AID_ORD(aid))
+		}
+		call id_fitdata (id)
+		call id_fitfeatures (id)
+		
+		j = 0
+		do i = 1, ID_NFEATURES(id) {
+		    if (IS_INDEFD(USER(id,i)) || WTS(id,i) != 0.) {
+			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
+	    }
+	    ID_MAXFEATURES(id) = maxfeatures
+	} then {
+	    ID_MAXFEATURES(id) = maxfeatures
+	    call erract (EA_ERROR)
+	}
+	if (IS_INDEFD(cdelt))
+	    return
+
+	nmin = 2
+        nfound = AID_NFOUND(aid)
+	if (ID_NFEATURES(id) < nfound)
+	    call error (0, "aid_dofit: not enough features")
+
+	# Compute fwidth rms.
+	rms = 0.
+	for (i=1; i<=ID_NFEATURES(id); i=i+1)
+	    rms = rms + (FIT(id,i) - USER(id,i)) ** 2
+	rms = sqrt (rms/ max (1, ID_NFEATURES(id)-nmin)) / abs (cdelt)
+	rms = rms / ID_FWIDTH(id)
+
+	# Compute line list matching fraction.
+	ncandidate = max (nfound, (ncandidate-(nmatch1-nmatch2)))
+	fmatch = 1 - real (nmatch2) / ncandidate
+
+	# Compute target line matching fraction.
+	xt = AID_XT(aid)
+	nt = AID_NT(aid)
+	label = NULL
+	ntmatch = 0
+	do i = 1, nt {
+	    fit = id_fitpt (id, Memd[xt+i-1])
+	    user = INDEFD
+	    call id_match (id, fit, user, label, ID_MATCH(id))
+	    if (!IS_INDEFD(user))
+		ntmatch = ntmatch + 1
+	}
+	ftmatch = 1 - real (ntmatch) / min (nt, ncandidate)
+	call mfree (label, TY_CHAR)
+
+	if (AID_RMSG(aid) > 0. && AID_FMATCHG(aid) > 0.) {
+	    best = AID_WRMS(aid) *  rms / AID_RMSG(aid)
+	    best = best + AID_WFMATCH(aid) * fmatch / AID_FMATCHG(aid)
+	    best = best + AID_WFTMATCH(aid) * ftmatch / AID_FMATCHG(aid)
+	} else
+	    best = MAX_DOUBLE
+end
+
+
+# AID_DOFITF -- From a set of candidate identifications fit and evaluate
+# a final dispersion solution.
+
+procedure aid_dofitf (aid, id)
+
+pointer	aid		#I AID pointer
+pointer	id		#I ID pointer
+
+int	i, j, k, maxfeatures
+
+errchk	id_dofit, id_fitdata, id_fitfeatures, id_linelist
+
+begin
+	maxfeatures = ID_MAXFEATURES(id)
+	ID_MAXFEATURES(id) = 1000
+	iferr {
+	    do k = 1, 3 {
+		if (ID_NFEATURES(id) < 2)
+		    call error (0, "aid_dofit: not enough features")
+		if (k > 1)
+		    call id_linelist (id)
+
+		call id_dofit (id, NO)
+		call id_fitdata (id)
+		call id_fitfeatures (id)
+		if (k < 3) {
+		    j = 0
+		    do i = 1, ID_NFEATURES(id) {
+			if (IS_INDEFD(USER(id,i)) || WTS(id,i) != 0.) {
+			    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
+		}
+	    }
+	    ID_MAXFEATURES(id) = maxfeatures
+	} then {
+	    ID_MAXFEATURES(id) = maxfeatures
+	    call erract (EA_ERROR)
+	}
+end
+
+
+# AID_EVALLOC -- Allocate memory to save the candidate identifications
+# and dispersions to be evaluated.
+
+pointer procedure aid_evalloc (aid, index)
+
+pointer	aid		#I AID pointer
+int	index		#I Reference sample index
+
+begin
+	if (AID_EVS(aid) == NULL)
+	    call calloc (AID_EVS(aid), (index+49)*AID_EVLEN, TY_POINTER)
+	else if (index > 1 && mod (index-1, 50) == 0) {
+	    call realloc (AID_EVS(aid), (index+49)*AID_EVLEN, TY_POINTER)
+	    call aclri (Memi[AID_EVS(aid)+(index-1)*AID_EVLEN], 50*AID_EVLEN)
+	}
+	return (AID_EVS(aid)+(index-1)*AID_EVLEN)
+end
+
+
+# AID_EVFREE -- Free memory from the evaluation step.
+
+procedure aid_evfree (aid, index)
+
+pointer	aid		#I AID pointer
+int	index		#I Reference sample index
+
+pointer	ev, aid_evalloc()
+
+begin
+	ev = aid_evalloc (aid, index)
+	call mfree (AID_X(ev), TY_REAL)
+	call mfree (AID_Y(ev), TY_REAL)
+	call mfree (AID_A(ev), TY_REAL)
+	call mfree (AID_B(ev), TY_REAL)
+end
+
+
+# AID_IMGD -- Get value from image header or parameter string.
+
+double procedure aid_imgd (im, param)
+
+pointer	im			#I IMIO pointer
+char	param[ARB]		#I Parameter
+
+int	i, ctod()
+double	dval, imgetd()
+
+begin
+	if (param[1] == '!') {
+	    iferr (dval = imgetd (im, param[2]))
+		dval = INDEFD
+	} else {
+	    iferr (dval = imgetd (im, param)) {
+		i = 1
+		if (ctod (param, i, dval) == 0)
+		    dval = INDEFD
+	    }
+	}
+	return (dval)
+end