aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/identify/autoid
diff options
context:
space:
mode:
Diffstat (limited to 'noao/onedspec/identify/autoid')
-rw-r--r--noao/onedspec/identify/autoid/aidautoid.x314
-rw-r--r--noao/onedspec/identify/autoid/aidget.x21
-rw-r--r--noao/onedspec/identify/autoid/aidgraph.x240
-rw-r--r--noao/onedspec/identify/autoid/aidinit.x93
-rw-r--r--noao/onedspec/identify/autoid/aidlog.x57
-rw-r--r--noao/onedspec/identify/autoid/aidset.x162
-rw-r--r--noao/onedspec/identify/autoid/aidshift.x67
-rw-r--r--noao/onedspec/identify/autoid/autoid.h90
-rw-r--r--noao/onedspec/identify/autoid/autoid.x1600
-rw-r--r--noao/onedspec/identify/autoid/mkpkg17
10 files changed, 2661 insertions, 0 deletions
diff --git a/noao/onedspec/identify/autoid/aidautoid.x b/noao/onedspec/identify/autoid/aidautoid.x
new file mode 100644
index 00000000..7c213b4a
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidautoid.x
@@ -0,0 +1,314 @@
+include <mach.h>
+include <gset.h>
+include <math/iminterp.h>
+include <smw.h>
+include "../identify.h"
+include "autoid.h"
+
+
+# List of debug key letters.
+# Debug a: Print candidate line assignments.
+# Debug b: Print search limits.
+# Debug c: Print list of line ratios.
+# Debug d: Graph dispersions.
+# Debug f: Print final result.
+# Debug i: Show ICFIT iterations.
+# Debug l: Graph lines and spectra.
+# Debug m: Print miscellaneous debug information
+# Debug n: Show non-linearity constraint
+# Debug r: Print list of reference lines.
+# Debug s: Print search iterations.
+# Debug t: Print list of target lines.
+# Debug v: Print vote array.
+# Debug w: Print wavelength bin limits.
+
+
+# AID_AUTOID -- Automatically identify spectral features.
+# This routine is main entry to the autoidentify algorithms.
+# The return value is true if the ID pointer contains a new solution
+# and false if the ID pointer is the original solution.
+
+bool procedure aid_autoid (id, aid)
+
+pointer id #I ID pointer
+pointer aid #U AID pointer
+
+bool cdflip
+int i, j, k, l, iev, nbins, bin, nbest
+double best, dval1, dval2
+pointer sp, str, idr, ev, evf, sid
+
+bool streq(), strne()
+int stridxs()
+double dcveval(), aid_eval()
+pointer gopen(), aid_evalloc(), id_getid()
+errchk id_mapll, aid_reference, aid_target, aid_autoid1, aid_evalutate
+
+define done_ 10
+define redo_ 20
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Save original data.
+ call id_saveid (id, "autoidentify backup")
+
+ # Initialize.
+ AID_IDT(aid) = id
+ call ic_putr (AID_IC1(aid), "xmin", real (PIXDATA(id,1)))
+ call ic_putr (AID_IC1(aid), "xmax", real (PIXDATA(id,ID_NPTS(id))))
+ AID_IC2(aid) = ID_IC(id)
+
+ if (stridxs ("ild", AID_DEBUG(aid,1)) != 0 && ID_GP(id) == NULL) {
+ call clgstr ("graphics", Memc[str], SZ_LINE)
+ ID_GP(id) = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH)
+ } else if (AID_DEBUG(aid,1) != EOS && ID_GP(id) != NULL)
+ call gdeactivate (ID_GP(id), 0)
+
+ idr = AID_IDR(aid)
+ if (idr == NULL) {
+ call id_init (AID_IDR(aid))
+ idr = AID_IDR(aid)
+ }
+
+ # Set reference coordinate list.
+ if (strne (AID_REFLIST(aid), ID_COORDLIST(idr)) ||
+ streq (AID_REFLIST(aid), "FEATURES")) {
+ call id_unmapll (idr)
+ ID_COORDLIST(idr) = EOS
+
+ if (streq (AID_REFLIST(aid), "FEATURES")) {
+ if (ID_NFEATURES(id) >= 10) {
+ call strcpy (AID_REFLIST(aid), ID_COORDLIST(idr),
+ ID_LENSTRING)
+ i = ID_NFEATURES(id)
+ ID_NLL(idr) = i
+ call calloc (ID_LL(idr), i+1, TY_DOUBLE)
+ call calloc (ID_LLL(idr), i+1, TY_POINTER)
+ call amovd (USER(id,1), Memd[ID_LL(idr)], i)
+ Memd[ID_LL(idr)+i] = INDEFD
+ }
+ } else if (AID_REFLIST(aid) != EOS) {
+ call strcpy (AID_REFLIST(aid), ID_COORDLIST(idr), ID_LENSTRING)
+ call id_mapll (idr)
+ }
+ }
+
+ # Get reference spectrum.
+ if (AID_REFSPEC(aid) != EOS)
+ call strcpy (AID_REFSPEC(aid), ID_COORDSPEC(idr), ID_LENSTRING)
+ else if (ID_COORDSPEC(idr) == EOS)
+ call strcpy (ID_COORDSPEC(id), ID_COORDSPEC(idr), ID_LENSTRING)
+ if (strne (ID_COORDSPEC(idr), ID_IMAGE(idr))) {
+ if (ID_SH(idr) != NULL) {
+ call smw_close (MW(ID_SH(idr)))
+ call imunmap (IM(ID_SH(idr)))
+ call shdr_close (ID_SH(idr))
+ }
+ call strcpy (ID_COORDSPEC(idr), ID_IMAGE(idr), ID_LENSTRING)
+ ifnoerr (call id_map (idr))
+ call id_gdata(idr)
+ else {
+ ID_COORDSPEC(idr) = EOS
+ ID_IMAGE(idr) = EOS
+ }
+ }
+
+ ID_MAXFEATURES(idr) = AID_NRMAX(aid)
+ ID_MINSEP(idr) = ID_MINSEP(id)
+ ID_FTYPE(idr) = ID_FTYPE(id)
+ ID_FWIDTH(idr) = ID_FWIDTH(id)
+ ID_CRADIUS(idr) = ID_CRADIUS(id)
+ ID_THRESHOLD(idr) = ID_THRESHOLD(id)
+ ID_MATCH(idr) = ID_MATCH(id)
+
+ # Use faster, less accurate centering for the search.
+ call c1d_params (II_LINEAR, 0.02)
+
+ # Set target lines and dispersion limits.
+ call aid_target (aid)
+ cdflip = (AID_CDDIR(aid) == CDUNKNOWN ||
+ (IS_INDEFD(AID_CDELT(aid)) && AID_CDDIR(aid) == CDSIGN))
+
+ # Now search for the dispersion function and line identifications.
+ # The reference spectrum is broken up into a varying number of
+ # pieces and each is searched. The order in which the reference
+ # spectrum is divided is from the middle outward and overlapping
+ # bins are tried as a second pass. The hope is to find a
+ # piece that is close enough to the target spectrum as quickly
+ # as possible.
+
+ AID_BEST(aid) = MAX_REAL
+ nbest = 0
+ iev = 0
+redo_
+ do i = 0, 1 {
+ do j = 1, AID_NB(aid) {
+ if (j == 1)
+ nbins = (AID_NB(aid) + 2) / 2
+ else if (mod (j, 2) == 0)
+ nbins = (AID_NB(aid) + 2 - j) / 2
+ else
+ nbins = (AID_NB(aid) + 1 + j) / 2
+ nbins = 2 * nbins - 1
+ do k = 1, nbins {
+ if (k == 1)
+ bin = (nbins + 2) / 2
+ else if (mod (k, 2) == 0)
+ bin = (nbins + 2 - k) / 2
+ else
+ bin = (nbins + 1 + k) / 2
+ if (mod ((nbins-1)/2, 2) == 0) {
+ if (mod (bin, 2) == i)
+ next
+ } else {
+ if (mod (bin, 2) != i)
+ next
+ }
+
+ iferr {
+ iev = iev + 1
+ ev = aid_evalloc (aid, iev)
+ AID_BIN1(ev) = nbins
+ AID_BIN2(ev) = bin
+ call aid_reference (aid, ev, NO)
+ call aid_autoid1 (aid, ev)
+ } then {
+ AID_ND(ev) = 0
+ }
+ if (cdflip) {
+ iferr {
+ iev = iev + 1
+ evf = aid_evalloc (aid, iev)
+ AID_BIN1(evf) = nbins
+ AID_BIN2(evf) = bin
+ call aid_reference (aid, evf, YES)
+ call aid_autoid1 (aid, evf)
+ } then {
+ AID_ND(evf) = 0
+ }
+ }
+
+ # Search the candidates with the highest weights.
+ # Terminate the search if the number of best fit values
+ # less than 1. is equal to the specified number.
+ do l = 1, 5 {
+ best = aid_eval (aid, ev, l)
+ if (!IS_INDEFD(best) && best < 1.) {
+ nbest = nbest + 1
+ if (nbest >= AID_NBEST(aid))
+ goto done_
+ }
+ if (cdflip) {
+ best = aid_eval (aid, evf, l)
+ if (!IS_INDEFD(best) && best < 1.) {
+ nbest = nbest + 1
+ if (nbest >= AID_NBEST(aid))
+ goto done_
+ }
+ }
+ }
+ }
+ }
+ }
+
+ # Go back and evaluate candidates with lower weights.
+ # Terminate the search if the number of best fit values
+ # less than 1. is equal to the specified number.
+ do j = 6, AID_ND(ev) {
+ do i = 1, iev {
+ ev = aid_evalloc (aid, i)
+ best = aid_eval (aid, ev, j)
+ if (!IS_INDEFD(best) && best < 1.) {
+ nbest = nbest + 1
+ if (nbest >= AID_NBEST(aid))
+ goto done_
+ }
+ }
+ }
+
+ # Add other loops here.
+ if (AID_BEST(aid) > 1.) {
+ if (AID_NP(aid) > 3) {
+ AID_NP(aid) = AID_NP(aid) - 1
+ goto redo_
+ }
+ }
+
+done_
+ do i = 1, iev
+ call aid_evfree (aid, i)
+
+ # Evaluate the final solution with the full dispersion function.
+ # Save the final solution. If the final solution has a merit
+ # greater than one restore the original solution.
+
+ sid = id_getid (id, "autoidentify")
+ if (sid != NULL) {
+ call dcvfree (ID_CV(id))
+ iferr (call aid_dofitf (aid, id))
+ ;
+ call id_sid (id, sid)
+ } else {
+ ID_NFEATURES(id) = 0
+ call dcvfree (ID_CV(id))
+ call id_saveid (id, "autoidentify")
+ }
+
+ # Debug f: Print final result.
+ if (stridxs ("f", AID_DEBUG(aid,1)) != 0) {
+ if (AID_BEST(aid) == MAX_REAL) {
+ call eprintf ("%s %8.5g %8.3g No solution found\n")
+ call pargstr (ID_IMAGE(id))
+ call pargd (AID_CRVAL(aid))
+ call pargd (AID_CDELT(aid))
+ } else {
+ call eprintf (
+ "%s %8.5g %8.3g %8.5g %8.3g %3d %3d %6.3f %5.2f\n")
+ call pargstr (ID_IMAGE(id))
+ call pargd (AID_CRVAL(aid))
+ call pargd (AID_CDELT(aid))
+ if (ID_CV(id) == NULL) {
+ dval1 = FITDATA(id,1)
+ dval2 = FITDATA(id,2) - FITDATA(id,1)
+ } else {
+ dval1 = dcveval (ID_CV(id), AID_CRPIX(aid)+1D0)
+ dval2 = dcveval (ID_CV(id), AID_CRPIX(aid)-1D0)
+ dval2 = (dval1 - dval2) / 2D0
+ dval1 = dcveval (ID_CV(id), AID_CRPIX(aid))
+ }
+ call pargd (dval1)
+ call pargd (FITDATA(id,2) - FITDATA(id,1))
+ call pargi (nint(100.*AID_FMATCH(aid)))
+ call pargi (nint(100.*AID_FTMATCH(aid)))
+ call pargr (AID_RMS(aid))
+ call pargr (AID_BEST(aid))
+ call eprintf (
+ " dCRVAL = %.4g (%.3g), dCDELT = %.4g (%.3g)\n")
+ call pargd (dval1 - AID_CRVAL(aid))
+ call pargd (abs((dval1-AID_CRVAL(aid))/
+ (ID_NPTS(id)*AID_CDELT(aid))))
+ call pargd (dval2 - AID_CDELT(aid))
+ call pargd (abs((dval2-AID_CDELT(aid))/AID_CDELT(aid)))
+ }
+ }
+
+ if (AID_BEST(aid) > 1.) {
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ call dcvfree (ID_CV(id))
+ sid = id_getid (id, "autoidentify backup")
+ ID_NEWFEATURES(id) = NO
+ ID_NEWCV(id) = NO
+ ID_NEWGRAPH(id) = NO
+ }
+ call id_fitdata (id)
+
+ # Restore centering.
+ call c1d_params (II_SPLINE3, 0.001)
+
+ call sfree (sp)
+
+ return (AID_BEST(aid) <= 1.)
+end
diff --git a/noao/onedspec/identify/autoid/aidget.x b/noao/onedspec/identify/autoid/aidget.x
new file mode 100644
index 00000000..ba3c9342
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidget.x
@@ -0,0 +1,21 @@
+include "autoid.h"
+
+define AIDGET "|best|"
+
+
+# AID_GETR -- Get AID parameters by name.
+
+real procedure aid_getr (aid, param)
+
+pointer aid #I AID object
+char param[ARB] #I Parameter name
+
+char temp[10]
+int strdic()
+
+begin
+ switch (strdic (param, temp, 10, AIDGET)) {
+ case 1:
+ return (AID_BEST(aid))
+ }
+end
diff --git a/noao/onedspec/identify/autoid/aidgraph.x b/noao/onedspec/identify/autoid/aidgraph.x
new file mode 100644
index 00000000..35494004
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidgraph.x
@@ -0,0 +1,240 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include <smw.h>
+include "../identify.h"
+include "autoid.h"
+
+
+# AID_LGRAPH -- Graph target and reference spectra and associated lines.
+# This is only used for debugging.
+
+procedure aid_lgraph (aid, x1, n1, x2, n2)
+
+pointer aid #I AID pointer
+double x1[n1] #I Reference lines
+int n1 #I Number of reference lines
+double x2[n2] #I Target lines
+int n2 #I Number of target lines
+
+int i, wcs, key, nr, nt, redraw, clgcur(), stridxs()
+real wx, wy, wz, a, b, c, d, dy, ytmin, ytmax
+pointer sp, cmd, id, sht, shr, gp, gt, xr, yr, yt, y, gt_init()
+double shdr_lw()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ id = AID_IDT(aid)
+ sht = ID_SH(id)
+ shr = ID_SH(AID_IDR(aid))
+
+ gp = ID_GP(id)
+ if (gp == NULL)
+ return
+ gt = gt_init()
+ call gt_sets (gt, GTTYPE, "line")
+ call gt_seti (gt, GTSYSID, NO)
+ if (DC(sht) == DCNO || WP(sht) * AID_CDSIGN(aid) < 0.) {
+ call gt_setr (gt, GTXMIN, AID_W1(aid))
+ call gt_setr (gt, GTXMAX, AID_W2(aid))
+ } else {
+ call gt_setr (gt, GTXMIN, W0(sht))
+ call gt_setr (gt, GTXMAX, W1(sht))
+ }
+
+ if (shr != NULL) {
+ xr = SX(shr) + AID_X1R(aid) - 1
+ yr = AID_SPECR(aid)
+ nr = AID_X2R(aid) - AID_X1R(aid) + 1
+ }
+
+ nt = ID_NPTS(id)
+ yt = ID_IMDATA(id)
+ call alimr (Memr[yt], nt, ytmin, ytmax)
+
+ call malloc (y, max(nr,nt), TY_REAL)
+
+ key = 'r'
+ repeat {
+ switch (key) {
+ case ':':
+ call gt_colon (Memc[cmd], gp, gt, redraw)
+ case 'Q':
+ i = stridxs ("l", AID_DEBUG(aid,1))
+ AID_DEBUG(aid,i) = ' '
+ break
+ case 'q':
+ break
+ case 'r':
+ redraw = YES
+ case 'w':
+ call gt_window (gt, gp, "gcur", redraw)
+ }
+
+ if (redraw == YES) {
+ call gclear (gp)
+ call gseti (gp, G_YDRAWTICKS, NO)
+ if (shr != NULL) {
+ call gascale (gp, Memr[xr], nr, 1)
+ call gascale (gp, Memr[yr], nr, 2)
+ } else {
+ call gswind (gp, real(x1[1]), real(x1[n1]), 0., 1.)
+ }
+ call gt_swind (gp, gt)
+ call ggwind (gp, a, b, c, d)
+ dy = 2 * (d - c)
+ call gswind (gp, a, b, c, c + dy)
+ call gt_labax(gp, gt)
+
+ if (shr != NULL) {
+ call aminkr (Memr[yr], c + 0.44 * dy, Memr[y], nr)
+ call gt_plot (gp, gt, Memr[xr], Memr[y], nr)
+ }
+
+ wy = c + 0.46 * dy
+ wz = c + 0.49 * dy
+ do i = 1, n1 {
+ wx = x1[i]
+ if (wx < min (a,b) || wx > max (a,b))
+ next
+ call gline (gp, wx, wy, wx, wz)
+ }
+
+ call amapr (Memr[yt], Memr[y], nt,
+ ytmin, ytmax, c+0.55*dy, c+0.99*dy)
+ wy = c + 0.50 * dy
+ wz = c + 0.53 * dy
+
+ if (DC(sht) == DCNO || WP(sht) * AID_CDSIGN(aid) < 0.) {
+ call gvline (gp, Memr[y], nt, a, b)
+ do i = 1, n2 {
+ wx = a + (x2[i] - 1) / (nt - 1) * (b - a)
+ call gline (gp, wx, wy, wx, wz)
+ }
+ } else {
+ call gpline (gp, Memr[SX(sht)], Memr[y], nt)
+ do i = 1, n2 {
+ wx = shdr_lw (sht, double (x2[i]))
+ call gline (gp, wx, wy, wx, wz)
+ }
+ }
+
+ redraw = NO
+ }
+ } until (clgcur ("gcur", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+
+ call gdeactivate (gp, 0)
+ call mfree (y, TY_REAL)
+ call gt_free (gt)
+ call sfree (sp)
+end
+
+
+# AID_DGRAPH -- Graph candidate dispersions.
+# This routine is only used for debugging.
+
+procedure aid_dgraph (aid, x, y, n, w1, dw, nd)
+
+pointer aid #I AID pointer
+real x[n] #I Array of candidate reference coordinates (sorted)
+real y[n] #I Array of candidate target coordinates
+int n #I Number of candidate pairs
+real w1[nd] #I Dispersion origin
+real dw[nd] #I Dispersion slope
+int nd #I Number of dispersions
+
+int i, ndplot, wcs, key, redraw, clgcur(), stridxs(), ctoi()
+real wx, wy, a, b, c, d, dy, crpix, crval, cdelt
+pointer sp, cmd, sh, gp, gt, gt_init()
+
+begin
+ call smark (sp)
+ call salloc (cmd, SZ_LINE, TY_CHAR)
+
+ sh = ID_SH(AID_IDT(aid))
+ gp = ID_GP(AID_IDT(aid))
+ if (gp == NULL)
+ return
+ gt = gt_init()
+ call gt_seti (gt, GTSYSID, NO)
+ if (DC(sh) != DCNO) {
+ call gt_setr (gt, GTXMIN, W0(sh))
+ call gt_setr (gt, GTXMAX, W1(sh))
+ call gt_setr (gt, GTYMIN, 1.)
+ call gt_setr (gt, GTYMAX, real(SN(sh)))
+ }
+
+ ndplot = nd
+ key = 'r'
+ repeat {
+ switch (key) {
+ case ':':
+ if (Memc[cmd] == '/')
+ call gt_colon (Memc[cmd], gp, gt, redraw)
+ else {
+ i = 1
+ if (ctoi (Memc[cmd], i, ndplot) == 0)
+ ndplot = nd
+ }
+ case 'Q':
+ i = stridxs ("d", AID_DEBUG(aid,1))
+ AID_DEBUG(aid,i) = ' '
+ break
+ case 'q':
+ break
+ case 'r':
+ redraw = YES
+ case 'w':
+ call gt_window (gt, gp, "gcur", redraw)
+ }
+
+ if (redraw == YES) {
+ call gclear (gp)
+ call gascale (gp, x, n, 1)
+ call gascale (gp, y, n, 2)
+ call gt_swind (gp, gt)
+ call gt_labax(gp, gt)
+
+ call gt_plot (gp, gt, x, y, n)
+
+ call ggwind (gp, a, b, c, d)
+ dy = (b - a) / 500.
+ do i = 1, ndplot {
+ crval = w1[i]
+ cdelt = dw[i]
+ wy = c
+ wx = crval + wy * cdelt
+ call gamove (gp, wx, wy)
+ for (wy=wy+dy; wy<d+dy; wy=wy+dy) {
+ wx = crval + wy * cdelt
+ call gadraw (gp, wx, wy)
+ }
+ }
+
+ if (AID_CRMIN(aid) > -MAX_DOUBLE / 10. &&
+ AID_CRMAX(aid) < MAX_DOUBLE / 10.) {
+ crpix = AID_CRPIX(aid)
+ crval = AID_CDSIGN(aid) * AID_CDMIN(aid)
+ cdelt = AID_CDSIGN(aid) * AID_CDMAX(aid)
+ for (wy=c; wy<d+dy; wy=wy+dy) {
+ wx = AID_CRMIN(aid) +
+ min ((wy-crpix)*crval, (wy-crpix)*cdelt)
+ call gmark (gp, wx, wy, GM_POINT, 2, 2)
+ }
+ for (wy=c; wy<d+dy; wy=wy+dy) {
+ wx = AID_CRMAX(aid) +
+ max ((wy-crpix)*crval, (wy-crpix)*cdelt)
+ call gmark (gp, wx, wy, GM_POINT, 2, 2)
+ }
+ }
+
+ redraw = NO
+ }
+ } until (clgcur ("gcur", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+
+ call gdeactivate (gp, 0)
+ call gt_free (gt)
+ call sfree (sp)
+end
diff --git a/noao/onedspec/identify/autoid/aidinit.x b/noao/onedspec/identify/autoid/aidinit.x
new file mode 100644
index 00000000..ac86b34d
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidinit.x
@@ -0,0 +1,93 @@
+include <smw.h>
+include "../identify.h"
+include "autoid.h"
+
+
+# AID_INIT -- Create AID object and initialize algorithm parameters.
+
+procedure aid_init (aid, pset)
+
+pointer aid #O AID object
+char pset[ARB] #I Pset for parameters
+
+pointer pp, clopset()
+int clgpseti(), strdic()
+double clgpsetd()
+
+begin
+ call calloc (aid, AID_LEN, TY_STRUCT)
+
+ # Set default parameters. This can be overridden later by
+ # the application.
+
+ pp = clopset (pset)
+
+ #call clgpseta (pp, "crval", AID_CR(aid), AID_SZLINE)
+ #call clgpseta (pp, "cdelt", AID_CD(aid), AID_SZLINE)
+ call strcpy ("INDEF", AID_CR(aid), AID_SZLINE)
+ call strcpy ("INDEF", AID_CD(aid), AID_SZLINE)
+
+ call clgpseta (pp, "reflist", AID_REFLIST(aid), AID_SZLINE)
+ call clgpseta (pp, "refspec", AID_REFSPEC(aid), AID_SZLINE)
+ call clgpseta (pp, "crpix", AID_CP(aid), AID_SZLINE)
+ call clgpseta (pp, "crquad", AID_CQ(aid), AID_SZLINE)
+ call clgpseta (pp, "cddir", AID_DEBUG(aid,1), AID_SZLINE)
+ AID_CDDIR(aid) = strdic (AID_DEBUG(aid,1), AID_DEBUG(aid,1),
+ AID_SZLINE, CDDIR)
+ call clgpseta (pp, "crsearch", AID_CRS(aid), AID_SZLINE)
+ call clgpseta (pp, "cdsearch", AID_CDS(aid), AID_SZLINE)
+ AID_NTMAX(aid) = clgpseti (pp, "ntarget")
+ #AID_NRMAX(aid) = clgpseti (pp, "nreference")
+ AID_NRMAX(aid) = 2 * AID_NTMAX(aid)
+ AID_ORD(aid) = clgpseti (pp, "aidord")
+ AID_MAXNL(aid) = clgpsetd (pp, "maxnl")
+ AID_NB(aid) = clgpseti (pp, "nbins")
+ AID_NN(aid) = clgpseti (pp, "nneighbors")
+ AID_NP(aid) = clgpseti (pp, "npattern")
+ AID_SIG(aid) = clgpsetd (pp, "sigma")
+ AID_NFOUND(aid) = clgpseti (pp, "nfound")
+ AID_RMSG(aid) = clgpsetd (pp, "rms")
+ AID_FMATCHG(aid) = clgpsetd (pp, "fmatch")
+ AID_FTMATCHG(aid) = clgpsetd (pp, "fmatch")
+ AID_MINRATIO(aid) = clgpsetd (pp, "minratio")
+ AID_NDMAX(aid) = clgpseti (pp, "ndmax")
+ call clgpseta (pp, "debug", AID_DEBUG(aid,1), AID_SZLINE)
+ AID_NBEST(aid) = 3
+ AID_WRMS(aid) = 0.34
+ AID_WFMATCH(aid) = 0.33
+ AID_WFTMATCH(aid) = 0.33
+ call clcpset (pp)
+
+ call ic_open (AID_IC1(aid))
+ call ic_pstr (AID_IC1(aid), "function", "chebyshev")
+ call ic_puti (AID_IC1(aid), "order", AID_ORD(aid))
+ call ic_puti (AID_IC1(aid), "niterate", 5)
+ call ic_putr (AID_IC1(aid), "low", 2.)
+ call ic_putr (AID_IC1(aid), "high", 2.)
+end
+
+
+# AID_FREE -- Free memory associated with the AID algorithms.
+
+procedure aid_free (aid)
+
+pointer aid #U AID object
+
+begin
+ if (AID_IDR(aid) != NULL) {
+ if (ID_SH(AID_IDR(aid)) != NULL) {
+ call smw_close (MW(ID_SH(AID_IDR(aid))))
+ call imunmap (IM(ID_SH(AID_IDR(aid))))
+ call shdr_close (ID_SH(AID_IDR(aid)))
+ }
+ }
+
+ call ic_closed (AID_IC1(aid))
+ call mfree (AID_SPECR(aid), TY_REAL)
+ call mfree (AID_XR(aid), TY_DOUBLE)
+ call mfree (AID_XT(aid), TY_DOUBLE)
+ call mfree (AID_XTF(aid), TY_DOUBLE)
+ call id_free (AID_IDR(aid))
+ call mfree (AID_EVS(aid), TY_POINTER)
+ call mfree (aid, TY_STRUCT)
+end
diff --git a/noao/onedspec/identify/autoid/aidlog.x b/noao/onedspec/identify/autoid/aidlog.x
new file mode 100644
index 00000000..b0247d00
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidlog.x
@@ -0,0 +1,57 @@
+include "../identify.h"
+
+
+# AID_LOG -- Log final solution.
+
+procedure aid_log (id, fd, hdr)
+
+pointer id #I ID object
+int fd #I Log file descriptor
+int hdr #U Print header?
+
+double wc, dw, id_fitpt(), id_rms()
+pointer str
+bool fp_equald()
+
+begin
+ if (fd == NULL)
+ return
+
+ if (fd == STDOUT && ID_GP(id) != NULL)
+ call gdeactivate (ID_GP(id), 0)
+
+ if (hdr == YES) {
+ call malloc (str, SZ_LINE, TY_CHAR)
+ call sysid (Memc[str], SZ_LINE)
+ call fprintf (fd, "\nAUTOIDENTIFY: %s\n")
+ call pargstr (Memc[str])
+ call mfree (str, TY_CHAR)
+
+ call fprintf (fd, " %-20s %10s %10s %10s %10s\n")
+ call pargstr ("Spectrum")
+ call pargstr ("# Found")
+ call pargstr ("Midpoint")
+ call pargstr ("Dispersion")
+ call pargstr ("RMS")
+
+ hdr = NO
+ }
+
+ call fprintf (fd, " %s%s%24t ")
+ call pargstr (ID_IMAGE(id))
+ call pargstr (ID_SECTION(id))
+ if (ID_CV(id) == NULL)
+ call fprintf (fd, " No solution found\n")
+ else {
+ wc = id_fitpt (id, (ID_NPTS(id) + 1D0) / 2D0)
+ dw = wc - id_fitpt (id, (ID_NPTS(id) - 1D0) / 2D0)
+ if (!fp_equald (dw, 0D0)) {
+ call fprintf (fd, "%10d %10.*g %10.3g %10.3g\n")
+ call pargi (ID_NFEATURES(id))
+ call pargi (int (log10 (abs (wc / dw)) + 3))
+ call pargd (wc)
+ call pargd (dw)
+ call pargd (id_rms(id))
+ }
+ }
+end
diff --git a/noao/onedspec/identify/autoid/aidset.x b/noao/onedspec/identify/autoid/aidset.x
new file mode 100644
index 00000000..5905002b
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidset.x
@@ -0,0 +1,162 @@
+include "autoid.h"
+
+define AIDSET "|reflist|refspec|crval|cdelt|crpix|crquad|crsearch|cdsearch\
+ |cddir|ntarget|nreference|aidord|maxnl|nbins|nneighbors\
+ |npattern|sigma|nfound|rms|fmatch|ftmatch|minratio|ndmax\
+ |debug|nbest|wrms|wfmatch|wftmatch|"
+
+
+# AID_SETS -- Set AID parameters by name.
+# If the first word of the value field is "CL" or "ENV" then the second
+# word is the CL parameter name or environment variable name to use
+# for the value.
+
+procedure aid_sets (aid, param, value)
+
+pointer aid #I AID object
+char param[ARB] #I Parameter name
+char value[ARB] #I Value
+
+int i, j, strdic(), strncmp(), envfind(), nowhite(), ctoi(), ctor(), ctod()
+pointer sp, str
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ i = strdic (param, Memc[str], SZ_LINE, AIDSET)
+
+ if (strncmp ("CL ", value, 3) == 0)
+ call clgstr (value[4], Memc[str], SZ_LINE)
+ else if (strncmp ("ENV ", value, 4) == 0) {
+ if (envfind (value[5], Memc[str], SZ_LINE) <= 0)
+ Memc[str] = EOS
+ } else
+ call strcpy (value, Memc[str], SZ_LINE)
+ j = nowhite (Memc[str], Memc[str], SZ_LINE)
+
+ j = 1
+ switch (i) {
+ case 1:
+ call strcpy (Memc[str], AID_REFLIST(aid), AID_SZLINE)
+ case 2:
+ call strcpy (Memc[str], AID_REFSPEC(aid), AID_SZLINE)
+ case 3:
+ call strcpy (Memc[str], AID_CR(aid), AID_SZLINE)
+ case 4:
+ call strcpy (Memc[str], AID_CD(aid), AID_SZLINE)
+ case 5:
+ call strcpy (Memc[str], AID_CP(aid), AID_SZLINE)
+ case 6:
+ i = ctod (Memc[str], j, AID_CRQUAD(aid))
+ case 7:
+ call strcpy (Memc[str], AID_CRS(aid), AID_SZLINE)
+ case 8:
+ call strcpy (Memc[str], AID_CDS(aid), AID_SZLINE)
+ case 9:
+ AID_CDDIR(aid) = strdic (Memc[str], Memc[str], SZ_LINE, CDDIR)
+ if (AID_CDDIR(aid) == 0)
+ AID_CDDIR(aid) = CDUNKNOWN
+ case 10:
+ i = ctoi (Memc[str], j, AID_NTMAX(aid))
+ case 11:
+ i = ctoi (Memc[str], j, AID_NRMAX(aid))
+ case 12:
+ i = ctoi (Memc[str], j, AID_ORD(aid))
+ call ic_puti (AID_IC1(aid), "order", AID_ORD(aid))
+ case 13:
+ i = ctor (Memc[str], j, AID_MAXNL(aid))
+ case 14:
+ i = ctoi (Memc[str], j, AID_NB(aid))
+ case 15:
+ i = ctoi (Memc[str], j, AID_NN(aid))
+ case 16:
+ i = ctoi (Memc[str], j, AID_NP(aid))
+ case 17:
+ i = ctor (Memc[str], j, AID_SIG(aid))
+ case 18:
+ i = ctoi (Memc[str], j, AID_NFOUND(aid))
+ case 19:
+ i = ctor (Memc[str], j, AID_RMSG(aid))
+ case 20:
+ i = ctor (Memc[str], j, AID_FMATCHG(aid))
+ case 21:
+ i = ctor (Memc[str], j, AID_FTMATCHG(aid))
+ case 22:
+ i = ctor (Memc[str], j, AID_MINRATIO(aid))
+ case 23:
+ i = ctoi (Memc[str], j, AID_NDMAX(aid))
+ case 24:
+ call strcpy (Memc[str], AID_DEBUG(aid,1), AID_SZLINE)
+ case 25:
+ i = ctoi (Memc[str], j, AID_NBEST(aid))
+ case 26:
+ i = ctor (Memc[str], j, AID_WRMS(aid))
+ case 27:
+ i = ctor (Memc[str], j, AID_WFMATCH(aid))
+ case 28:
+ i = ctor (Memc[str], j, AID_WFTMATCH(aid))
+ }
+
+ call sfree (sp)
+end
+
+
+# AID_SETI -- Set integer AID parameters.
+
+procedure aid_seti (aid, param, ival)
+
+pointer aid #I AID object
+char param[ARB] #I Parameter name
+int ival #I Value
+
+pointer sp, str
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call sprintf (Memc[str], SZ_FNAME, "%d")
+ call pargi (ival)
+ call aid_sets (aid, param, Memc[str])
+ call sfree (sp)
+end
+
+
+# AID_SETR -- Set real AID parameters.
+
+procedure aid_setr (aid, param, rval)
+
+pointer aid #I AID object
+char param[ARB] #I Parameter name
+real rval #I Value
+
+pointer sp, str
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call sprintf (Memc[str], SZ_FNAME, "%g")
+ call pargr (rval)
+ call aid_sets (aid, param, Memc[str])
+ call sfree (sp)
+end
+
+
+# AID_SETD -- Set double AID parameters.
+
+procedure aid_setd (aid, param, dval)
+
+pointer aid #I AID object
+char param[ARB] #I Parameter name
+double dval #I Value
+
+pointer sp, str
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_FNAME, TY_CHAR)
+ call sprintf (Memc[str], SZ_FNAME, "%g")
+ call pargd (dval)
+ call aid_sets (aid, param, Memc[str])
+ call sfree (sp)
+end
diff --git a/noao/onedspec/identify/autoid/aidshift.x b/noao/onedspec/identify/autoid/aidshift.x
new file mode 100644
index 00000000..1b910338
--- /dev/null
+++ b/noao/onedspec/identify/autoid/aidshift.x
@@ -0,0 +1,67 @@
+include "../identify.h"
+
+
+# AID_SHIFT -- Find a new shifted dispersion solution assuming (nearly) the
+# same dispersion per pixel and the same dispersion direction. The shift is
+# assumed to be less than or equal to the dispersion range of the input
+# dispersion. The input is an ID pointer have the previous dispersion
+# solution and features but with the new spectrum. If there are more than 10
+# features then the list of user feature coordinates is used as the reference
+# list. If there are not enough features or the previous search fails then
+# the the coordinate list is used as the reference. The returned result is a
+# new ID pointer if the algorithm succeeds or the original ID pointer if it
+# fails along with an error status.
+
+procedure aid_shift (id, crsearch, cdsearch)
+
+pointer id #I ID object
+double crsearch #I Search range
+double cdsearch #I Search range
+
+pointer aid
+bool found, aid_autoid()
+double crpix, crval, cdelt, id_fitpt()
+
+begin
+ # Set approximate dispersion from input dispersion solution.
+ crpix = ID_NPTS(id) / 2 + 1
+ crval = id_fitpt (id, crpix)
+ cdelt = (FITDATA(id,ID_NPTS(id)) - FITDATA(id,1)) /
+ (ID_NPTS(id) - 1)
+
+ # Initialize AUTOID.
+ call aid_init (aid, "aidpars")
+ call aid_setd (aid, "crval", crval)
+ call aid_setd (aid, "cdelt", cdelt)
+ call aid_setd (aid, "crpix", crpix)
+ call aid_sets (aid, "cddir", "sign")
+ call aid_setd (aid, "crsearch", crsearch)
+ call aid_setd (aid, "cdsearch", cdsearch)
+ call aid_seti (aid, "nbest", 5)
+
+ found = false
+ if (ID_NFEATURES(id) > 10) {
+ # Try shift using features.
+ call aid_seti (aid, "ntarget", ID_NFEATURES(id))
+ call aid_seti (aid, "nreference", ID_NFEATURES(id))
+ call aid_setr (aid, "wrms", 0.5)
+ call aid_setr (aid, "wfmatch", 0.5)
+ call aid_setr (aid, "wftmatch", 0.)
+ call aid_sets (aid, "refspec", "FEATURES")
+ found = aid_autoid (id, aid)
+ }
+ if (!found) {
+ # Try shift using coordinate list.
+ call aid_seti (aid, "ntarget", max (ID_NFEATURES(id),20))
+ call aid_seti (aid, "nreference", max (ID_NFEATURES(id),40))
+ call aid_setr (aid, "wrms", 0.5)
+ call aid_setr (aid, "wfmatch", 0.25)
+ call aid_setr (aid, "wftmatch", 0.25)
+ call aid_sets (aid, "refspec", "COORDLIST")
+ found = aid_autoid (id, aid)
+ }
+
+ call aid_free (aid)
+ if (!found)
+ call error (1, "No solution not found")
+end
diff --git a/noao/onedspec/identify/autoid/autoid.h b/noao/onedspec/identify/autoid/autoid.h
new file mode 100644
index 00000000..304d675a
--- /dev/null
+++ b/noao/onedspec/identify/autoid/autoid.h
@@ -0,0 +1,90 @@
+# AUTOIDENTIFY data structure.
+
+define AID_SZLINE 99
+define AID_LEN 512
+
+# Algorithm input parameters.
+define AID_REFLIST Memc[P2C($1)] # Reference coordinate list
+define AID_REFSPEC Memc[P2C($1+50)] # Reference spectrum
+define AID_CR Memc[P2C($1+100)] # Coordinate reference value
+define AID_CD Memc[P2C($1+150)] # Coordinate reference value
+define AID_CP Memc[P2C($1+200)] # Coordinate reference value
+define AID_CQ Memc[P2C($1+250)] # Coordinate quad distortion
+define AID_CRS Memc[P2C($1+300)] # Coordinate reference value
+define AID_CDS Memc[P2C($1+350)] # Coordinate reference value
+define AID_DEBUG Memc[P2C($1+400)+ 2-1] # Debug flags (19 chars)
+define AID_CDDIR Memi[$1+450] # Coordinate direction
+define AID_NTMAX Memi[$1+451] # Maximum number of target lines
+define AID_NRMAX Memi[$1+452] # Maximum number of reference lines
+define AID_ORD Memi[$1+453] # Maximum fitting order
+define AID_MAXNL Memr[P2R($1+454)] # Maximum non-linearity
+define AID_NB Memi[$1+455] # Number of sub-bins
+define AID_NN Memi[$1+456] # Number of neighbor lines
+define AID_NP Memi[$1+457] # Number of lines in pattern
+define AID_SIG Memr[P2R($1+458)] # Target line centering sigma
+define AID_NFOUND Memi[$1+459] # Minimum number to be found
+define AID_RMSG Memr[P2R($1+460)] # Pixel RMS (goal)
+define AID_FMATCHG Memr[P2R($1+461)] # Frac of unmatched lines (goal)
+define AID_FTMATCHG Memr[P2R($1+462)] # Frac of unmatched target lines (goal)
+
+define AID_IDT Memi[$1+463] # Target ID pointer
+define AID_IDR Memi[$1+464] # Reference ID pointer
+define AID_IC1 Memi[$1+465] # ICFIT pointer
+define AID_IC2 Memi[$1+466] # ICFIT pointer
+
+define AID_XR Memi[$1+467] # Reference lines (ptr)
+define AID_NR Memi[$1+468] # Number of reference lines
+define AID_XTF Memi[$1+469] # Full target lines sorted by peak
+define AID_NTF Memi[$1+470] # Full number of target lines
+define AID_XT Memi[$1+471] # Target lines to use sorted by pix
+define AID_XTL Memi[$1+472] # Linearized target lines sort by pix
+define AID_NT Memi[$1+473] # Number of target lines to use
+
+define AID_CDSIGN Memi[$1+474] # Sign of coordinate interval
+define AID_CRVAL Memd[P2D($1+476)] # Reference coordinate value
+define AID_CDELT Memd[P2D($1+478)] # Coordinate interval per pixel
+define AID_CRPIX Memd[P2D($1+480)] # Reference pixel
+define AID_CRQUAD Memd[P2D($1+482)] # Quadratic distortion
+define AID_CRSEARCH Memd[P2D($1+484)] # Search radius for ref value
+define AID_CDSEARCH Memd[P2D($1+486)] # Search radius for coord int
+define AID_CRMIN Memd[P2D($1+488)] # Min for central coordinate
+define AID_CRMAX Memd[P2D($1+490)] # Max for central coordinate
+define AID_CDMIN Memd[P2D($1+492)] # Min for coordinate interval
+define AID_CDMAX Memd[P2D($1+494)] # Max for coordinate interval
+
+define AID_MINRATIO Memr[P2R($1+496)] # Minimum ratio
+define AID_NDMAX Memi[$1+497] # Max number of dispersions to check
+define AID_RMS Memr[P2R($1+498)] # Pixel RMS (best)
+define AID_FMATCH Memr[P2R($1+499)] # Fraction of unmatched linelist lines
+define AID_FTMATCH Memr[P2R($1+500)] # Fraction of unmatched target lines
+define AID_WRMS Memr[P2R($1+501)] # Weight for RMS
+define AID_WFMATCH Memr[P2R($1+502)] # Weight for FMATCH
+define AID_WFTMATCH Memr[P2R($1+503)] # Weight for FTMATCH
+define AID_NBEST Memi[$1+504] # Number of best values < 1 to check
+define AID_BEST Memr[P2R($1+505)] # Best fit parameter
+define AID_EVS Memi[$1+506] # Evaluate structure
+
+define AID_SPECR Memi[$1+507] # Reference spectrum (ptr)
+define AID_X1R Memi[$1+508] # First pixel of full ref spectrum
+define AID_X2R Memi[$1+509] # Last pixel of full ref spectrum
+define AID_W1 Memr[P2R($1+510)] # Tentative wavelength of first pixel
+define AID_W2 Memr[P2R($1+511)] # Tentative wavelength of last pixel
+
+
+# Evaluation structure.
+define AID_EVLEN 8
+define AID_BIN1 Memi[$1] # Reference sample bin
+define AID_BIN2 Memi[$1+1] # Reference sample bin
+define AID_X Memi[$1+2] # Pixel coordinates
+define AID_Y Memi[$1+3] # Dispersion coordinates
+define AID_N Memi[$1+4] # Number of coordinate pairs
+define AID_A Memi[$1+5] # Trial dispersion start
+define AID_B Memi[$1+6] # Trial dispersion step
+define AID_ND Memi[$1+7] # Number of trial dispersions
+
+# Dispersion direction options.
+define CDDIR "|sign|increasing|decreasing|unknown|"
+define CDSIGN 1
+define CDINC 2
+define CDDEC 3
+define CDUNKNOWN 4
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
diff --git a/noao/onedspec/identify/autoid/mkpkg b/noao/onedspec/identify/autoid/mkpkg
new file mode 100644
index 00000000..7d46d183
--- /dev/null
+++ b/noao/onedspec/identify/autoid/mkpkg
@@ -0,0 +1,17 @@
+$checkout libpkg.a ../../
+$update libpkg.a
+$checkin libpkg.a ../../
+$exit
+
+libpkg.a:
+ aidautoid.x autoid.h <gset.h> ../identify.h <mach.h>\
+ <math/iminterp.h> <smw.h>
+ aidget.x autoid.h
+ aidgraph.x autoid.h <gset.h> ../identify.h <mach.h> <pkg/gtools.h>\
+ <smw.h>
+ aidinit.x autoid.h ../identify.h <smw.h>
+ aidlog.x ../identify.h
+ aidset.x autoid.h
+ aidshift.x ../identify.h
+ autoid.x autoid.h <error.h> ../identify.h <mach.h> <smw.h>
+ ;
generated by cgit (git 2.34.1) at 2025-09-20 13:48:10 -0400