aboutsummaryrefslogtreecommitdiff
path: root/noao/rv/rvidlines
diff options
context:
space:
mode:
authorJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
committerJoseph Hunkeler <jhunkeler@gmail.com>2015-07-08 20:46:52 -0400
commitfa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch)
treebdda434976bc09c864f2e4fa6f16ba1952b1e555 /noao/rv/rvidlines
downloadiraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz
Initial commit
Diffstat (limited to 'noao/rv/rvidlines')
-rw-r--r--noao/rv/rvidlines/Revisions52
-rw-r--r--noao/rv/rvidlines/idcenter.x257
-rw-r--r--noao/rv/rvidlines/idcolon.x288
-rw-r--r--noao/rv/rvidlines/iddb.x436
-rw-r--r--noao/rv/rvidlines/iddeblend.x413
-rw-r--r--noao/rv/rvidlines/iddelete.x26
-rw-r--r--noao/rv/rvidlines/iddofit.x101
-rw-r--r--noao/rv/rvidlines/iddoshift.x42
-rw-r--r--noao/rv/rvidlines/identify.h97
-rw-r--r--noao/rv/rvidlines/identify.key104
-rw-r--r--noao/rv/rvidlines/idfitdata.x140
-rw-r--r--noao/rv/rvidlines/idfixx.x27
-rw-r--r--noao/rv/rvidlines/idgdata.x74
-rw-r--r--noao/rv/rvidlines/idgraph.x168
-rw-r--r--noao/rv/rvidlines/ididentify.x795
-rw-r--r--noao/rv/rvidlines/idinit.x352
-rw-r--r--noao/rv/rvidlines/idlabel.x30
-rw-r--r--noao/rv/rvidlines/idlinelist.x250
-rw-r--r--noao/rv/rvidlines/idlog.x190
-rw-r--r--noao/rv/rvidlines/idmap.x379
-rw-r--r--noao/rv/rvidlines/idmark.x97
-rw-r--r--noao/rv/rvidlines/idnearest.x29
-rw-r--r--noao/rv/rvidlines/idnewfeature.x87
-rw-r--r--noao/rv/rvidlines/idnoextn.x11
-rw-r--r--noao/rv/rvidlines/idpeak.x23
-rw-r--r--noao/rv/rvidlines/idrms.x28
-rw-r--r--noao/rv/rvidlines/idshift.x65
-rw-r--r--noao/rv/rvidlines/idshow.x83
-rw-r--r--noao/rv/rvidlines/idvelocity.x188
-rw-r--r--noao/rv/rvidlines/idvhelio.x102
-rw-r--r--noao/rv/rvidlines/mkpkg47
-rw-r--r--noao/rv/rvidlines/peaks.gx447
-rw-r--r--noao/rv/rvidlines/peaks.x446
-rw-r--r--noao/rv/rvidlines/reidentify.x609
-rw-r--r--noao/rv/rvidlines/rvidlines.key100
-rw-r--r--noao/rv/rvidlines/t_identify.x108
-rw-r--r--noao/rv/rvidlines/t_reidentify.x1092
37 files changed, 7783 insertions, 0 deletions
diff --git a/noao/rv/rvidlines/Revisions b/noao/rv/rvidlines/Revisions
new file mode 100644
index 00000000..9f93511a
--- /dev/null
+++ b/noao/rv/rvidlines/Revisions
@@ -0,0 +1,52 @@
+idinit.x
+reidentify.x
+ Changed aclrr to aclri and updated the 'b' key in reidentify.x.
+ (12/1/93, Valdes)
+
+idlog.x
+ Added the spectrum title and heliocentric Julian date to the log output.
+ (11/17/93, Valdes)
+
+iddblend.x +
+ Needed copy of ONEDSPEC deblend with names changes to protect the
+ innocent and avoid conflict with the version in RV. (8/25/93, Valdes)
+
+t_identify.x
+t_reidentify.x
+idcenter.x
+idcolon.x
+iddb.x
+iddofit.x
+iddoshift.x
+idgdata.x
+idgraph.x
+ididentify.x
+idinit.x
+idlinelist.x
+idlog.x
+idmark.x
+idrms.x
+idshift.x
+idshow.x
+idvelocity.x +
+idvhelio.x
+peaks.gx
+reidentify.x
+mkpkg
+identify.h
+identify.par
+reidentify.par
+rvidlines.par +
+rvreidlines.par +
+identify.key
+rvidlines.key +
+doc/rvidlines.hlp +
+doc/rvreidlines.hlp +
+ 1. A new cursor key was added for marking and deblending features.
+ 2. Two new feature types (ftype) were added:
+ gabsorption Absorption features fit by gaussian
+ gemission Emission features fit by gaussian
+ 3. The 'u' key was modified to allow setting a feature weight
+ 4. Two new tasks were added RVIDLINES and RVREIDLINES. These
+ are implemented as special entry points and a task flag.
+ (8/23/93, Valdes)
diff --git a/noao/rv/rvidlines/idcenter.x b/noao/rv/rvidlines/idcenter.x
new file mode 100644
index 00000000..05c636cb
--- /dev/null
+++ b/noao/rv/rvidlines/idcenter.x
@@ -0,0 +1,257 @@
+include <gset.h>
+include <smw.h>
+include "identify.h"
+
+# ID_CENTER -- Locate the center of a feature.
+
+double procedure id_center (id, x, n, width, type, interactive)
+
+pointer id # ID pointer
+double x[n] # Initial guess
+int n # Number of features
+real width # Feature width
+int type # Feature type
+int interactive # Interactive?
+
+int np1
+real value
+
+real center1d()
+double smw_c1trand()
+
+begin
+ switch (type) {
+ case EMISSION, ABSORPTION:
+ np1 = NP1(ID_SH(id)) - 1
+ value = smw_c1trand (ID_PL(id), x[1]) - np1
+ value = center1d (value, IMDATA(id,1), ID_NPTS(id),
+ width, type, ID_CRADIUS(id), ID_THRESHOLD(id))
+ if (IS_INDEF(value))
+ return (INDEFD)
+ else
+ return (smw_c1trand (ID_LP(id), double(value+np1)))
+ case GEMISSION, GABSORPTION:
+ iferr (call id_gcenter (id, x, n, INDEF, INDEF, INDEF, INDEF,
+ width, type, interactive))
+ return (INDEFD)
+ return (x[1])
+ }
+end
+
+
+# ID_GCENTER -- Locate the center of a feature.
+
+procedure id_gcenter (id, x, ng, xa, ya, xb, yb, width, type, interactive)
+
+pointer id # ID pointer
+double x[ng] # Initial guess
+int ng # Number of features
+real xa, ya, xb, yb # Background points
+real width # Feature width
+int type # Feature type
+int interactive # Draw gaussian fit?
+
+int i, np1, x1, x2
+real ag, bg, pix, w, y
+pointer sp, xr, xg, yg, sg, gp
+
+bool fp_equalr()
+real id_model()
+double smw_c1trand(), id_fitpt()
+
+errchk gcenter1d
+
+begin
+ call smark (sp)
+ call salloc (xr, ng, TY_REAL)
+ call salloc (xg, ng, TY_REAL)
+ call salloc (yg, ng, TY_REAL)
+ call salloc (sg, ng, TY_REAL)
+
+ np1 = NP1(ID_SH(id)) - 1
+
+ # Compute background in logical units.
+ if (IS_INDEF(xa) || IS_INDEF(ya) || IS_INDEF(xb) || IS_INDEF(yb)) {
+ ag = INDEF
+ bg = INDEF
+ } else {
+ ag = smw_c1trand (ID_PL(id), double(xa)) - np1
+ bg = smw_c1trand (ID_PL(id), double(xb)) - np1
+ if (!fp_equalr (ag, bg)) {
+ bg = (yb - ya) / (bg - ag)
+ ag = ya - bg * ag
+ } else {
+ ag = INDEF
+ bg = INDEF
+ }
+ }
+
+ do i = 1, ng
+ Memr[xr+i-1] = smw_c1trand (ID_PL(id), x[i]) - np1
+ call gcenter1d (Memr[xr], ng, IMDATA(id,1), ID_NPTS(id),
+ width, type, ID_CRADIUS(id), ID_THRESHOLD(id),
+ x1, x2, Memr[xg], Memr[yg], Memr[sg], ag, bg)
+ do i = 1, ng
+ x[i] = smw_c1trand (ID_LP(id), double(Memr[xg+i-1]+np1))
+
+ if (interactive == YES) {
+ gp = ID_GP(id)
+ call gseti (gp, G_PLTYPE, 2)
+ call gseti (gp, G_PLCOLOR, 2)
+ pix = x1
+ w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1)))
+ y = id_model (pix, Memr[xg], Memr[yg], Memr[sg], ng) + ag +
+ bg * pix
+ call gamove (gp, w, y)
+ for (pix = x1; pix <= x2; pix = pix + .1) {
+ w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1)))
+ y = id_model (pix, Memr[xg], Memr[yg], Memr[sg], ng) +
+ ag + bg * pix
+ call gadraw (gp, w, y)
+ }
+ call gseti (gp, G_PLTYPE, 3)
+ call gseti (gp, G_PLCOLOR, 3)
+ pix = x1
+ w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1)))
+ y = ag + bg * pix
+ call gamove (gp, w, y)
+ pix = x2
+ w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1)))
+ y = ag + bg * pix
+ call gadraw (gp, w, y)
+ call gseti (gp, G_PLTYPE, 1)
+ call gseti (gp, G_PLCOLOR, 1)
+ call gflush (gp)
+ }
+
+ call sfree (sp)
+end
+
+
+define MIN_WIDTH 3. # Minimum centering width
+
+
+# GCENTER1D -- Locate the center of a one dimensional feature by guassian fit.
+# A value of INDEF is returned in the centering fails for any reason.
+# This procedure just sets up the data and adjusts for emission or
+# absorption features. The actual centering is done by GFIT.
+# If width <= 1 return the nearest minima or maxima.
+
+procedure gcenter1d (x, ng, data, npts, width, type, radius, threshold,
+ x1, x2, xg, yg, sg, ag, bg)
+
+real x[ng] # Initial guess
+int ng # Number of gaussians
+real data[npts] # Data points
+int npts # Number of data points
+real width # Feature width
+int type # Feature type
+real radius # Centering radius
+real threshold # Minimum range in feature
+int x1, x2 # Fitting region
+real xg[ng], yg[ng], sg[ng], ag, bg # Gaussian parameters
+
+int i, nx, xa, xb
+real a, b, c, d, rad, wid, ya, yb, chisq
+pointer xfit
+
+errchk id_mr_dofit
+
+begin
+ # Check starting values.
+ do i = 1, ng
+ if (IS_INDEF(x[i]) || (x[i] < 1) || (x[i] > npts))
+ call error (1, "Invalid starting values")
+
+ # Set minimum width and error radius. The minimum in the error radius
+ # is for defining the data window. The user error radius is used to
+ # check for an error in the derived center at the end of the centering.
+
+ wid = max (width, MIN_WIDTH)
+ rad = max (2., radius)
+
+ # Determine the pixel value range around the initial center, including
+ # the width and error radius buffer. Check for a minimum range.
+
+ call alimr (x, ng, c, d)
+ x1 = max (1., c - wid / 2 - rad - wid)
+ x2 = min (real (npts), d + wid / 2 + rad + wid + 1)
+ nx = x2 - x1 + 1
+ call alimr (data[x1], nx, a, b)
+ if (b - a < threshold)
+ call error (1, "Data range below threshold")
+
+ # Allocate memory for the continuum subtracted data vector. The X
+ # range is just large enough to include the error radius and the
+ # half width.
+
+ x1 = max (1., c - wid / 2 - rad)
+ x2 = min (real (npts), d + wid / 2 + rad + 1)
+ nx = x2 - x1 + 1
+
+ # Make the centering data positive, subtract the continuum, and
+ # apply a threshold to eliminate noise spikes.
+
+ xa = nint(c)
+ ya = data[xa]
+ xb = nint(d)
+ yb = data[xb]
+ switch (type) {
+ case GEMISSION:
+ for (i = xa; i >= x1; i=i-1)
+ if (data[i] < ya) {
+ xa = i
+ ya = data[i]
+ }
+ for (i = xb; i <= x2; i=i+1)
+ if (data[i] < yb) {
+ xb = i
+ yb = data[i]
+ }
+ case GABSORPTION:
+ for (i = xa; i >= x1; i=i-1)
+ if (data[i] > ya) {
+ xa = i
+ ya = data[i]
+ }
+ for (i = xb; i <= x2; i=i+1)
+ if (data[i] > yb) {
+ xb = i
+ yb = data[i]
+ }
+ default:
+ call error (0, "Unknown feature type")
+ }
+
+ # Set initial gaussian parameters.
+ if (IS_INDEF(ag) || IS_INDEF(bg)) {
+ if (xa == xb)
+ call error (1, "Can't determine background")
+ bg = (yb-ya) / (xb-xa)
+ ag = ya - bg * xa
+ }
+ do i = 1, ng {
+ xg[i] = x[i]
+ yg[i] = data[nint(x[i])] - ag - bg * x[i]
+ sg[i] = width / 6.
+ }
+
+ # Determine the center.
+ call malloc (xfit, nx, TY_REAL)
+ do i = x1, x2
+ Memr[xfit+i-x1] = i
+
+ call id_mr_dofit (0, 3, 3, Memr[xfit], data[x1], nx,
+ ag, bg, xg, yg, sg, ng, chisq)
+
+ # Check user centering error radius.
+ do i = 1, ng {
+ if (!IS_INDEF(xg[i])) {
+ if (abs (x[i] - xg[i]) > radius)
+ call error (2, "Error radius exceeded")
+ }
+ }
+
+ # Free memory and return the center position.
+ call mfree (xfit, TY_REAL)
+end
diff --git a/noao/rv/rvidlines/idcolon.x b/noao/rv/rvidlines/idcolon.x
new file mode 100644
index 00000000..c976c66b
--- /dev/null
+++ b/noao/rv/rvidlines/idcolon.x
@@ -0,0 +1,288 @@
+include <gset.h>
+include <error.h>
+include <smw.h>
+include "identify.h"
+
+# List of colon commands.
+define CMDS "|show|features|image|nsum|database|read|write|add|coordlist|match\
+ |maxfeatures|minsep|zwidth|labels|fwidth|ftype|cradius|threshold|"
+
+define SHOW 1 # Show parameters
+define FEATURES 2 # Show list of features
+define IMAGE 3 # Set new image
+define NSUM 4 # Set the number of lines or columns to sum
+define DATABASE 5 # Set new database
+define READ 6 # Read database entry
+define WRITE 7 # Write database entry
+define ADD 8 # Add features from database
+define COORDLIST 9 # Set new coordinate list
+define MATCH 10 # Set coordinate list matching distance
+define MAXFEATURES 11 # Set maximum number of features for auto find
+define MINSEP 12 # Set minimum separation distance
+define ZWIDTH 13 # Set zoom window width
+define LABEL 14 # Set label type
+define WIDTH 15 # Set centering width
+define TYPE 16 # Set centering type
+define RADIUS 17 # Set centering radius
+define THRESHOLD 18 # Set the centering threshold
+
+# ID_COLON -- Respond to colon command.
+
+procedure id_colon (id, cmdstr, newimage, prfeature)
+
+pointer id # ID pointer
+char cmdstr[ARB] # Colon command
+char newimage[ARB] # New image name
+int prfeature # Print current feature on status line
+
+char cmd[SZ_LINE]
+int i, ncmd, ival[2]
+real rval[2]
+pointer im
+
+int nscan(), strdic()
+pointer immap()
+errchk immap, id_dbread, id_dbwrite, id_log
+
+begin
+ # Scan the command string and get the first word.
+ call sscan (cmdstr)
+ call gargwrd (cmd, SZ_LINE)
+ ncmd = strdic (cmd, cmd, SZ_LINE, CMDS)
+
+ switch (ncmd) {
+ case SHOW: # :show - show values of parameters
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call gdeactivate (ID_GP(id), AW_CLEAR)
+ call id_show (id, "STDOUT")
+ call greactivate (ID_GP(id), AW_PAUSE)
+ } else {
+ iferr (call id_show (id, cmd)) {
+ call erract (EA_WARN)
+ prfeature = NO
+ }
+ }
+ case FEATURES: # :features - list features
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call gdeactivate (ID_GP(id), AW_CLEAR)
+ call id_log (id, "STDOUT", NULL)
+ call greactivate (ID_GP(id), AW_PAUSE)
+ } else {
+ iferr (call id_log (id, cmd, NULL)) {
+ call erract (EA_WARN)
+ prfeature = NO
+ }
+ }
+ case IMAGE: # :image - set image to identify
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call printf ("image %s\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ prfeature = NO
+ } else {
+ call strcpy (cmd, newimage, SZ_FNAME)
+ iferr {
+ im = immap (newimage, READ_ONLY, 0)
+ call imunmap (im)
+ } then {
+ newimage[1] = EOS
+ call erract (EA_WARN)
+ prfeature = NO
+ }
+ }
+ case NSUM: # :nsum - set number of lines or columns to sum in image
+ call gargi (ival[1])
+ if (nscan() == 1) {
+ call printf ("nsum %d %d\n")
+ call pargi (ID_NSUM(id,1))
+ call pargi (ID_NSUM(id,2))
+ prfeature = NO
+ } else {
+ ID_NSUM(id,1) = ival[1]
+ call gargi (ival[2])
+ if (nscan() == 3)
+ ID_NSUM(id,2) = ival[2]
+ call smw_daxis (NULL, NULL, SMW_PAXIS(MW(ID_SH(id)),1),
+ ID_NSUM(id,1), ID_NSUM(id,2))
+ }
+ case DATABASE: # :database - set database
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call printf ("database %s\n")
+ call pargstr (Memc[ID_DATABASE(id)])
+ prfeature = NO
+ } else {
+ call strcpy (cmd, Memc[ID_DATABASE(id)], SZ_FNAME)
+ ID_NEWDBENTRY(id) = YES
+ }
+ case READ: # :read - read database entry
+ prfeature = NO
+ iferr {
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1)
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ NO, YES)
+ else {
+ call gargi (ival[1])
+ if (nscan() < 3)
+ ival[1] = ID_AP(id,1)
+ call gargi (ival[2])
+ if (nscan() < 4)
+ ival[2] = ID_AP(id,2)
+ call id_dbread (id, cmd, ival, NO, YES)
+ }
+ } then
+ call erract (EA_WARN)
+ case WRITE: # :write - write database entry
+ prfeature = NO
+ iferr {
+ ival[1] = ID_AP(id,1)
+ ival[2] = ID_AP(id,2)
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1)
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ival, YES)
+ else {
+ call gargi (ival[1])
+ if (nscan() < 3)
+ ival[1] = ID_AP(id,1)
+ call gargi (ival[2])
+ if (nscan() < 4)
+ ival[2] = ID_AP(id,2)
+ call id_dbwrite (id, cmd, ival, YES)
+ }
+ } then
+ call erract (EA_WARN)
+ case ADD: # :add - add features from database entry
+ prfeature = NO
+ iferr {
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1)
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ YES, YES)
+ else {
+ call gargi (ival[1])
+ if (nscan() < 3)
+ ival[1] = ID_AP(id,1)
+ call gargi (ival[2])
+ if (nscan() < 4)
+ ival[2] = ID_AP(id,2)
+ call id_dbread (id, cmd, ival, YES, YES)
+ }
+ } then
+ call erract (EA_WARN)
+ case COORDLIST: # :coordlist - set coordinate list
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ call printf ("coordlist %s\n")
+ call pargstr (Memc[ID_COORDLIST(id)])
+ prfeature = NO
+ } else {
+ call strcpy (cmd, Memc[ID_COORDLIST(id)], SZ_FNAME)
+ call id_unmapll (id)
+ call id_mapll (id)
+ }
+ case MATCH: # :match - set matching distance for coordinate list
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("match %g\n")
+ call pargr (ID_MATCH(id))
+ prfeature = NO
+ } else
+ ID_MATCH(id) = rval[1]
+ case MAXFEATURES: # :maxfeatures - set max num features for auto find
+ call gargi (ival[1])
+ if (nscan() == 1) {
+ call printf ("maxfeatures %d\n")
+ call pargi (ID_MAXFEATURES(id))
+ prfeature = NO
+ } else
+ ID_MAXFEATURES(id) = ival[1]
+ case MINSEP: # :minsep - set minimum feature separation allowed
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("minsep %g\n")
+ call pargr (ID_MINSEP(id))
+ prfeature = NO
+ } else
+ ID_MINSEP(id) = rval[1]
+ case ZWIDTH: # :zwidth - set zoom window width
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("zwidth %g\n")
+ call pargr (ID_ZWIDTH(id))
+ prfeature = NO
+ } else {
+ ID_ZWIDTH(id) = rval[1]
+ if (ID_GTYPE(id) == 2)
+ ID_NEWGRAPH(id) = YES
+ }
+ case LABEL: # :labels - set label type
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ switch (ID_LABELS(id)) {
+ case 2:
+ call printf ("labels index\n")
+ case 3:
+ call printf ("labels pixel\n")
+ case 4:
+ call printf ("labels coords\n")
+ case 5:
+ call printf ("labels user\n")
+ case 6:
+ call printf ("labels both\n")
+ default:
+ call printf ("labels none\n")
+ }
+ prfeature = NO
+ } else {
+ ID_LABELS(id) = strdic (cmd, cmd, SZ_LINE, LABELS)
+ do i = 1, ID_NFEATURES(id)
+ call id_mark (id, i)
+ }
+ case WIDTH: # :fwidth - set centering width
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("fwidth %g\n")
+ call pargr (ID_FWIDTH(id))
+ prfeature = NO
+ } else
+ ID_FWIDTH(id) = rval[1]
+ case TYPE: # :ftype - set centering type
+ call gargwrd (cmd, SZ_LINE)
+ if (nscan() == 1) {
+ switch (ID_FTYPE(id)) {
+ case EMISSION:
+ call printf ("ftype emission\n")
+ case ABSORPTION:
+ call printf ("ftype absorption\n")
+ case GEMISSION:
+ call printf ("ftype gemission\n")
+ case GABSORPTION:
+ call printf ("ftype gabsorption\n")
+ }
+ prfeature = NO
+ } else
+ ID_FTYPE(id) = strdic (cmd, cmd, SZ_LINE, FTYPES)
+ case RADIUS: # :cradius - set centering radius
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("cradius %g\n")
+ call pargr (ID_CRADIUS(id))
+ prfeature = NO
+ } else
+ ID_CRADIUS(id) = rval[1]
+ case THRESHOLD: # :threshold - set centering threshold
+ call gargr (rval[1])
+ if (nscan() == 1) {
+ call printf ("threshold %g\n")
+ call pargr (ID_THRESHOLD(id))
+ prfeature = NO
+ } else
+ ID_THRESHOLD(id) = rval[1]
+ default:
+ call printf ("Unrecognized or ambiguous command\007")
+ prfeature = NO
+ }
+end
diff --git a/noao/rv/rvidlines/iddb.x b/noao/rv/rvidlines/iddb.x
new file mode 100644
index 00000000..90bafea9
--- /dev/null
+++ b/noao/rv/rvidlines/iddb.x
@@ -0,0 +1,436 @@
+include <imset.h>
+include <math/curfit.h>
+include <smw.h>
+include "identify.h"
+
+# ID_DBREAD -- Read features data from the database.
+
+procedure id_dbread (id, name, ap, add, verbose)
+
+pointer id # ID pointer
+char name[SZ_LINE] # Image name
+int ap[2] # Aperture number
+int add # Add features?
+int verbose # Verbose flag
+
+double pix
+int i, j, k, ncoeffs, rec
+pointer dt, sp, coeffs, line, str, sh
+
+int dtgeti(), dcvstati(), dtlocate(), dtscan(), nscan()
+real dtgetr()
+double dcvstatd()
+
+errchk dtremap(), dtlocate(), dtgeti(), dtgad()
+
+begin
+ call smark (sp)
+ call salloc (line, SZ_LINE, TY_CHAR)
+
+ call strcpy ("id", Memc[line], SZ_LINE)
+ call imgcluster (name, Memc[line+2], SZ_LINE)
+ call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[line], READ_ONLY)
+
+ call id_dbsection (id, name, ap, Memc[ID_SECTION(id)], SZ_FNAME)
+ call sprintf (Memc[line], SZ_LINE, "identify %s%s")
+ call pargstr (name)
+ call pargstr (Memc[ID_SECTION(id)])
+
+ dt = ID_DT(id)
+ sh = ID_SH(id)
+ iferr (rec = dtlocate (dt, Memc[line])) {
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call sprintf (Memc[str], SZ_LINE, "Entry not found: %s")
+ call pargstr (Memc[line])
+ call error (0, Memc[str])
+ }
+
+ if (add == YES) {
+ j = dtgeti (dt, rec, "features")
+ k = j + ID_NFEATURES(id)
+
+ call realloc (ID_PIX(id), k, TY_DOUBLE)
+ call realloc (ID_FIT(id), k, TY_DOUBLE)
+ call realloc (ID_USER(id), k, TY_DOUBLE)
+ call realloc (ID_WTS(id), k, TY_DOUBLE)
+ call realloc (ID_FWIDTHS(id), k, TY_REAL)
+ call realloc (ID_FTYPES(id), k, TY_INT)
+ call realloc (ID_LABEL(id), k, TY_POINTER)
+
+ do i = 1, j {
+ k = dtscan (dt)
+ call gargd (pix)
+
+ ID_NFEATURES(id) = ID_NFEATURES(id) + 1
+ for (k=ID_NFEATURES(id); (k>1)&&(pix<PIX(id,k-1)); k=k-1) {
+ PIX(id,k) = PIX(id,k-1)
+ FIT(id,k) = FIT(id,k-1)
+ USER(id,k) = USER(id,k-1)
+ WTS(id,k) = WTS(id,k-1)
+ FWIDTH(id,k) = FWIDTH(id,k-1)
+ FTYPE(id,k) = FTYPE(id,k-1)
+ Memi[ID_LABEL(id)+k-1] = Memi[ID_LABEL(id)+k-2]
+ }
+ PIX(id,k) = pix
+ call gargd (FIT(id,k))
+ call gargd (USER(id,k))
+ call gargr (FWIDTH(id,k))
+ call gargi (FTYPE(id,k))
+ call gargd (WTS(id,k))
+ call gargstr (Memc[line], SZ_LINE)
+ Memi[ID_LABEL(id)+k-1] = NULL
+ call id_label (Memc[line], Memi[ID_LABEL(id)+k-1])
+
+ # The following initialization is for backwards compatibility.
+ if (nscan() < 5) {
+ FWIDTH(id,k) = ID_FWIDTH(id)
+ FTYPE(id,k) = ID_FTYPE(id)
+ } else if (nscan() < 6)
+ WTS(id,k) = 1.
+ }
+
+ } else {
+ if (SMW_FORMAT(MW(sh)) == SMW_ES || SMW_FORMAT(MW(sh)) == SMW_MS) {
+ iferr (APLOW(sh,1) = dtgetr (dt, rec, "aplow"))
+ APLOW(sh,1) = INDEF
+ iferr (APHIGH(sh,1) = dtgetr (dt, rec, "aphigh"))
+ APHIGH(sh,1) = INDEF
+ }
+
+ do i = 1, ID_NFEATURES(id)
+ call mfree (Memi[ID_LABEL(id)+i-1], TY_CHAR)
+
+ k = dtgeti (dt, rec, "features")
+ ID_NFEATURES(id) = k
+ ID_NALLOC(id) = k
+ call realloc (ID_PIX(id), k, TY_DOUBLE)
+ call realloc (ID_FIT(id), k, TY_DOUBLE)
+ call realloc (ID_USER(id), k, TY_DOUBLE)
+ call realloc (ID_WTS(id), k, TY_DOUBLE)
+ call realloc (ID_FWIDTHS(id), k, TY_REAL)
+ call realloc (ID_FTYPES(id), k, TY_INT)
+ call realloc (ID_LABEL(id), k, TY_POINTER)
+
+ do i = 1, ID_NFEATURES(id) {
+ k = dtscan (dt)
+ call gargd (PIX(id,i))
+ call gargd (FIT(id,i))
+ call gargd (USER(id,i))
+ call gargr (FWIDTH(id,i))
+ call gargi (FTYPE(id,i))
+ call gargd (WTS(id,i))
+ call gargstr (Memc[line], SZ_LINE)
+ Memi[ID_LABEL(id)+i-1] = NULL
+ call id_label (Memc[line], Memi[ID_LABEL(id)+i-1])
+
+ # The following initialization is for backwards compatibility.
+ if (nscan() < 5) {
+ FWIDTH(id,i) = ID_FWIDTH(id)
+ FTYPE(id,i) = ID_FTYPE(id)
+ } else if (nscan() < 6)
+ WTS(id,i) = 1.
+ }
+
+ iferr (ID_SHIFT(id) = dtgetr (dt, rec, "shift"))
+ ID_SHIFT(id) = 0.
+ iferr (ID_REDSHIFT(id) = dtgetr (dt, rec, "redshift"))
+ ID_REDSHIFT(id) = 0.
+ iferr (ID_RMSRED(id) = dtgetr (dt, rec, "redshift_rms"))
+ ID_RMSRED(id) = 0.
+
+ iferr {
+ ncoeffs = dtgeti (dt, rec, "coefficients")
+ call salloc (coeffs, ncoeffs, TY_DOUBLE)
+ call dtgad (dt, rec, "coefficients", Memd[coeffs], ncoeffs,
+ ncoeffs)
+
+ if (ID_CV(id) != NULL)
+ call dcvfree (ID_CV(id))
+ call dcvrestore (ID_CV(id), Memd[coeffs])
+
+ call ic_putr (ID_IC(id), "xmin", real (dcvstatd(ID_CV(id),
+ CVXMIN)))
+ call ic_putr (ID_IC(id), "xmax", real (dcvstatd(ID_CV(id),
+ CVXMAX)))
+ ifnoerr (call dtgstr (dt,rec,"function",Memc[line],SZ_LINE)) {
+ call ic_pstr (ID_IC(id), "function", Memc[line])
+ call ic_puti (ID_IC(id), "order", dtgeti (dt, rec, "order"))
+ call dtgstr (dt, rec, "sample", Memc[line], SZ_LINE)
+ call ic_pstr (ID_IC(id), "sample", Memc[line])
+ call ic_puti (ID_IC(id), "naverage",
+ dtgeti (dt, rec, "naverage"))
+ call ic_puti (ID_IC(id), "niterate",
+ dtgeti (dt, rec, "niterate"))
+ call ic_putr (ID_IC(id), "low",
+ dtgetr (dt, rec, "low_reject"))
+ call ic_putr (ID_IC(id), "high",
+ dtgetr (dt, rec, "high_reject"))
+ call ic_putr (ID_IC(id), "grow", dtgetr (dt, rec, "grow"))
+ } else {
+ call ic_puti (ID_IC(id), "order", dcvstati (ID_CV(id),
+ CVORDER))
+ switch (dcvstati (ID_CV(id), CVTYPE)) {
+ case LEGENDRE:
+ call ic_pstr (ID_IC(id), "function", "legendre")
+ case CHEBYSHEV:
+ call ic_pstr (ID_IC(id), "function", "chebyshev")
+ case SPLINE1:
+ call ic_pstr (ID_IC(id), "function", "spline1")
+ case SPLINE3:
+ call ic_pstr (ID_IC(id), "function", "spline3")
+ }
+ }
+
+ ID_NEWCV(id) = YES
+ ID_CURRENT(id) = min (1, ID_NFEATURES(id))
+ } then
+ ;
+ }
+
+ call sfree (sp)
+
+ if (ID_NFEATURES(id) > 0) {
+ ID_NEWGRAPH(id) = YES
+ ID_NEWFEATURES(id) = YES
+ ID_CURRENT(id) = 1
+ } else
+ ID_CURRENT(id) = 0
+
+ if (verbose == YES) {
+ call printf ("identify %s%s\n")
+ call pargstr (name)
+ call pargstr (Memc[ID_SECTION(id)])
+ }
+end
+
+
+# ID_DBWRITE -- Write features data to the database.
+
+procedure id_dbwrite (id, name, ap, verbose)
+
+pointer id # ID pointer
+char name[ARB] # Image name
+int ap[2] # Aperture number
+int verbose # Verbose flag
+
+int i, ncoeffs
+pointer dt, sp, coeffs, root, sh, im
+
+int dcvstati(), ic_geti()
+real ic_getr()
+
+errchk dtremap
+
+begin
+ call smark (sp)
+ call salloc (root, SZ_FNAME, TY_CHAR)
+
+ call strcpy ("id", Memc[root], SZ_FNAME)
+ call imgcluster (name, Memc[root+2], SZ_FNAME)
+ call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[root], APPEND)
+
+ call id_dbsection (id, name, ap, Memc[ID_SECTION(id)], SZ_FNAME)
+
+ sh = ID_SH(id)
+ dt = ID_DT(id)
+ call dtptime (dt)
+ call dtput (dt, "begin\tidentify %s%s\n")
+ call pargstr (name)
+ call pargstr (Memc[ID_SECTION(id)])
+ call dtput (dt, "\tid\t%s\n")
+ call pargstr (name)
+ call dtput (dt, "\ttask\tidentify\n")
+ call dtput (dt, "\timage\t%s%s\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ if (SMW_FORMAT(MW(sh)) == SMW_ES || SMW_FORMAT(MW(sh)) == SMW_MS) {
+ call dtput (dt, "\taperture\t%d\n")
+ call pargi (ID_AP(id,1))
+ call dtput (dt, "\taplow\t%g\n")
+ call pargr (APLOW(sh,1))
+ call dtput (dt, "\taphigh\t%g\n")
+ call pargr (APHIGH(sh,1))
+ }
+
+ call dtput (dt, "\tfeatures\t%d\n")
+ call pargi (ID_NFEATURES(id))
+ do i = 1, ID_NFEATURES(id) {
+ call dtput (dt, "\t %10.2f %10.8g %10.8g %5.1f %d %d %s\n")
+ call pargd (PIX(id,i))
+ call pargd (FIT(id,i))
+ call pargd (USER(id,i))
+ call pargr (FWIDTH(id,i))
+ call pargi (FTYPE(id,i))
+ call pargd (WTS(id,i))
+ if (Memi[ID_LABEL(id)+i-1] != NULL)
+ call pargstr (Memc[Memi[ID_LABEL(id)+i-1]])
+ else
+ call pargstr ("")
+ }
+
+ if (ID_SHIFT(id) != 0.) {
+ call dtput (dt, "\tshift\t%g\n")
+ call pargd (ID_SHIFT(id))
+ }
+ if (ID_REDSHIFT(id) != 0.) {
+ call dtput (dt, "\tredshift\t%g\n")
+ call pargd (ID_REDSHIFT(id))
+ call dtput (dt, "\tredshift_rms\t%g\n")
+ call pargd (ID_RMSRED(id))
+ }
+
+ if (ID_CV(id) != NULL) {
+ call dtput (dt, "\tfunction %s\n")
+ call ic_gstr (ID_IC(id), "function", Memc[root], SZ_FNAME)
+ call pargstr (Memc[root])
+ call dtput (dt, "\torder %d\n")
+ call pargi (ic_geti (ID_IC(id), "order"))
+ call dtput (dt, "\tsample %s\n")
+ call ic_gstr (ID_IC(id), "sample", Memc[root], SZ_FNAME)
+ call pargstr (Memc[root])
+ call dtput (dt, "\tnaverage %d\n")
+ call pargi (ic_geti (ID_IC(id), "naverage"))
+ call dtput (dt, "\tniterate %d\n")
+ call pargi (ic_geti (ID_IC(id), "niterate"))
+ call dtput (dt, "\tlow_reject %g\n")
+ call pargr (ic_getr (ID_IC(id), "low"))
+ call dtput (dt, "\thigh_reject %g\n")
+ call pargr (ic_getr (ID_IC(id), "high"))
+ call dtput (dt, "\tgrow %g\n")
+ call pargr (ic_getr (ID_IC(id), "grow"))
+
+ ncoeffs = dcvstati (ID_CV(id), CVNSAVE)
+ call salloc (coeffs, ncoeffs, TY_DOUBLE)
+ call dcvsave (ID_CV(id), Memd[coeffs])
+ call dtput (dt, "\tcoefficients\t%d\n")
+ call pargi (ncoeffs)
+ do i = 1, ncoeffs {
+ call dtput (dt, "\t\t%g\n")
+ call pargd (Memd[coeffs+i-1])
+ }
+ }
+
+ call dtput (dt, "\n")
+
+ ID_NEWFEATURES(id) = NO
+ ID_NEWCV(id) = NO
+ ID_NEWDBENTRY(id) = NO
+
+ if (ID_TASK(id) == IDENTIFY) {
+ if (verbose == YES) {
+ call printf ("identify %s%s\n")
+ call pargstr (name)
+ call pargstr (Memc[ID_SECTION(id)])
+ }
+
+ # Enter reference spectrum name in image header.
+ im = IM(sh)
+ call imseti (im, IM_WHEADER, YES)
+ call imastr (im, "REFSPEC1", Memc[ID_IMAGE(id)])
+ iferr (call imdelf (im, "REFSPEC2"))
+ ;
+ }
+
+ call sfree (sp)
+end
+
+
+# ID_DBCHECK -- Check if there is an entry in the database.
+# This does not actually read the database entry. It also assumes that
+# if a database is already open it is for the same image (the image
+# names are not checked) and the database has been scanned.
+
+int procedure id_dbcheck (id, name, ap)
+
+pointer id # ID pointer
+char name[SZ_LINE] # Image name
+int ap[2] # Aperture number
+
+int rec, stat
+pointer sp, line, sec
+
+int dtlocate()
+
+errchk dtremap(), dtlocate()
+
+begin
+ call smark (sp)
+ call salloc (line, SZ_LINE, TY_CHAR)
+ call salloc (sec, SZ_LINE, TY_CHAR)
+
+ if (ID_DT(id) == NULL) {
+ call strcpy ("id", Memc[line], SZ_LINE)
+ call imgcluster (name, Memc[line+2], SZ_LINE)
+ iferr (call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[line],
+ READ_ONLY)) {
+ call sfree (sp)
+ return (NO)
+ }
+ }
+
+ call id_dbsection (id, name, ap, Memc[sec], SZ_LINE)
+ call sprintf (Memc[line], SZ_LINE, "identify %s%s")
+ call pargstr (name)
+ call pargstr (Memc[sec])
+
+ iferr (rec = dtlocate (ID_DT(id), Memc[line]))
+ stat = NO
+ else
+ stat = YES
+
+ call sfree (sp)
+ return (stat)
+end
+
+
+# ID_DBSECTION -- Make the IDENTIFY section.
+
+procedure id_dbsection (id, name, ap, section, sz_section)
+
+pointer id #I ID pointer
+char name[SZ_LINE] #I Image name
+int ap[2] #I Aperture number
+char section[sz_section] #O IDENTIFY section
+int sz_section #I Size of section string
+
+pointer sh, smw
+bool streq()
+
+begin
+ sh = ID_SH(id)
+ smw = MW(sh)
+
+ switch (SMW_FORMAT(smw)) {
+ case SMW_ND:
+ section[1] = EOS
+ if (streq (name, Memc[ID_IMAGE(id)])) {
+ switch (SMW_LDIM(smw)) {
+ case 2:
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (section, sz_section, "[*,%d]")
+ case 2:
+ call sprintf (section, sz_section, "[%d,*]")
+ }
+ #call pargi (LINDEX(sh,1))
+ call pargi (ap[1])
+ case 3:
+ switch (SMW_LAXIS(smw,1)) {
+ case 1:
+ call sprintf (section, sz_section, "[*,%d,%d]")
+ case 2:
+ call sprintf (section, sz_section, "[%d,*,%d]")
+ case 3:
+ call sprintf (section, sz_section, "[%d,%d,*]")
+ }
+ #call pargi (LINDEX(sh,1))
+ #call pargi (LINDEX(sh,2))
+ call pargi (ap[1])
+ call pargi (ap[2])
+ }
+ }
+ case SMW_ES, SMW_MS:
+ call sprintf (section, sz_section, " - Ap %d")
+ call pargi (ap[1])
+ }
+end
diff --git a/noao/rv/rvidlines/iddeblend.x b/noao/rv/rvidlines/iddeblend.x
new file mode 100644
index 00000000..75e32a78
--- /dev/null
+++ b/noao/rv/rvidlines/iddeblend.x
@@ -0,0 +1,413 @@
+include <mach.h>
+
+
+# ID_MR_DOFIT -- Fit gaussian components. This is an interface to ID_DOFIT1
+# which puts parameters into the form required by ID_DOFIT1 and vice-versa.
+# It also implements a constrained approach to the solution.
+
+procedure id_mr_dofit (bkgfit, posfit, sigfit, x, y, npts, y1, dy, xg, yg, sg,
+ ng, chisq)
+
+int bkgfit # Fit background (0=no, 1=yes)
+int posfit # Position fitting flag (1=fixed, 2=single, 3=all)
+int sigfit # Sigma fitting flag (1=fixed, 2=single, 3=all)
+real x[npts] # X data
+real y[npts] # Y data
+int npts # Number of points
+real y1 # Continuum offset
+real dy # Continuum slope
+real xg[ng] # Initial and final x coordinates of gaussians
+real yg[ng] # Initial and final y coordinates of gaussians
+real sg[ng] # Initial and final sigmas of gaussians
+int ng # Number of gaussians
+real chisq # Chi squared
+
+int i
+pointer sp, a, j
+errchk id_dofit1
+
+begin
+ call smark (sp)
+ call salloc (a, 4 + 3 * ng, TY_REAL)
+
+ # Convert positions and widths relative to first component.
+ Memr[a] = y1
+ Memr[a+1] = dy
+ Memr[a+2] = xg[1]
+ Memr[a+3] = sg[1]
+ do i = 1, ng {
+ j = a + 3 * i + 1
+ Memr[j] = yg[i]
+ Memr[j+1] = xg[i] - Memr[a+2]
+ Memr[j+2] = sg[i] / Memr[a+3]
+ }
+
+ # Do fit.
+ do i = 0, bkgfit {
+ switch (10*posfit+sigfit) {
+ case 11:
+ call id_dofit1 (i, 1, 1, x, y, npts, Memr[a], ng, chisq)
+ case 12:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ case 13:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 1, 3, x, y, npts, Memr[a], ng, chisq)
+ case 21:
+ call id_dofit1 (i, 2, 1, x, y, npts, Memr[a], ng, chisq)
+ case 22:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq)
+ case 23:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 2, 3, x, y, npts, Memr[a], ng, chisq)
+ case 31:
+ call id_dofit1 (i, 2, 1, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 3, 1, x, y, npts, Memr[a], ng, chisq)
+ case 32:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 3, 2, x, y, npts, Memr[a], ng, chisq)
+ case 33:
+ call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 3, 2, x, y, npts, Memr[a], ng, chisq)
+ call id_dofit1 (i, 3, 3, x, y, npts, Memr[a], ng, chisq)
+ }
+ }
+
+ y1 = Memr[a]
+ dy = Memr[a+1]
+ do i = 1, ng {
+ j = a + 3 * i + 1
+ yg[i] = Memr[j]
+ xg[i] = Memr[j+1] + Memr[a+2]
+ sg[i] = abs (Memr[j+2] * Memr[a+3])
+ }
+
+ call sfree (sp)
+end
+
+
+# ID_MODEL -- Compute model.
+#
+# I(x) = I(i) exp {-[(x-xg(i)) / sg(i)]**2 / 2.}
+#
+# where the params are I1, I2, xg, yg, and sg.
+
+real procedure id_model (x, xg, yg, sg, ng)
+
+real x # X value to be evaluated
+real xg[ng] # X coordinates of gaussians
+real yg[ng] # Y coordinates of gaussians
+real sg[ng] # Sigmas of gaussians
+int ng # Number of gaussians
+
+int i
+real y, arg
+
+begin
+ y = 0.
+ do i = 1, ng {
+ arg = (x - xg[i]) / sg[i]
+ if (abs (arg) < 7.)
+ y = y + yg[i] * exp (-arg**2 / 2.)
+ }
+ return (y)
+end
+
+
+# ID_DOFIT1 -- Perform nonlinear iterative fit for the specified parameters.
+# This uses the Levenberg-Marquardt method from NUMERICAL RECIPES.
+
+procedure id_dofit1 (bkgfit, posfit, sigfit, x, y, npts, a, nlines, chisq)
+
+int bkgfit # Background fit (0=no, 1=yes)
+int posfit # Position fitting flag (1=fixed, 2=one, 3=all)
+int sigfit # Sigma fitting flag (1=fixed, 2=one, 3=all)
+real x[npts] # X data
+real y[npts] # Y data
+int npts # Number of points
+real a[ARB] # Fitting parameters
+int nlines # Number of lines
+real chisq # Chi squared
+
+int i, np, nfit
+real mr, chi2
+pointer sp, flags, ptr
+errchk id_mr_solve
+
+begin
+ # Number of terms is 3 for each line plus common background, center
+ # and sigma.
+
+ np = 3 * nlines + 4
+
+ call smark (sp)
+ call salloc (flags, np, TY_INT)
+ ptr = flags
+
+ # Background.
+ if (bkgfit == 1) {
+ Memi[ptr] = 1
+ Memi[ptr+1] = 2
+ ptr = ptr + 2
+ }
+
+ # Peaks are always fit.
+ do i = 1, nlines {
+ Memi[ptr] = 3 * i + 2
+ ptr = ptr + 1
+ }
+
+ # Positions.
+ switch (posfit) {
+ case 2:
+ Memi[ptr] = 3
+ ptr = ptr + 1
+ case 3:
+ Memi[ptr] = 3
+ ptr = ptr + 1
+ do i = 1, nlines {
+ Memi[ptr] = 3 * i + 3
+ ptr = ptr + 1
+ }
+ }
+
+ # Sigmas.
+ switch (sigfit) {
+ case 2:
+ Memi[ptr] = 4
+ ptr = ptr + 1
+ case 3:
+ Memi[ptr] = 4
+ ptr = ptr + 1
+ do i = 1, nlines {
+ Memi[ptr] = 3 * i + 4
+ ptr = ptr + 1
+ }
+ }
+
+ nfit = ptr - flags
+ mr = -1.
+ i = 0
+ chi2 = MAX_REAL
+ repeat {
+ call id_mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq)
+ if (chi2 - chisq > 1.)
+ i = 0
+ else
+ i = i + 1
+ chi2 = chisq
+ } until (i == 3)
+
+ mr = 0.
+ call id_mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq)
+
+ call sfree (sp)
+end
+
+
+# ID_DERIVS -- Compute model and derivatives for MR_SOLVE procedure.
+#
+# I(x) = I1 + I2 * x + I(i) exp {-[(x-xc-dx(i)) / (sig * sig(i))]**2 / 2.}
+#
+# where the params are I1, I2, xc, sig, I(i), dx(i), and sig(i) (i=1,nlines).
+
+procedure id_derivs (x, a, y, dyda, na)
+
+real x # X value to be evaluated
+real a[na] # Parameters
+real y # Function value
+real dyda[na] # Derivatives
+int na # Number of parameters
+
+int i
+real sig, arg, ex, fac
+
+begin
+ y = a[1] + a[2] * x
+ dyda[1] = 1.
+ dyda[2] = x
+ dyda[3] = 0.
+ dyda[4] = 0.
+ do i = 5, na, 3 {
+ sig = a[4] * a[i+2]
+ arg = (x - a[3] - a[i+1]) / sig
+ if (abs (arg) < 7.)
+ ex = exp (-arg**2 / 2.)
+ else
+ ex = 0.
+ fac = a[i] * ex * arg
+
+ y = y + a[i] * ex
+ dyda[3] = dyda[3] + fac / sig
+ dyda[4] = dyda[4] + fac * arg / a[4]
+ dyda[i] = ex
+ dyda[i+1] = fac / sig
+ dyda[i+2] = fac * arg / a[i+2]
+ }
+end
+
+
+# ID_MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization.
+#
+# Use the Levenberg-Marquardt method to minimize the chi squared of a set
+# of paraemters. The parameters being fit are indexed by the flag array.
+# To initialize the Marquardt parameter, MR, is less than zero. After that
+# the parameter is adjusted as needed. To finish set the parameter to zero
+# to free memory. This procedure requires a subroutine, DERIVS, which
+# takes the derivatives of the function being fit with respect to the
+# parameters. There is no limitation on the number of parameters or
+# data points. For a description of the method see NUMERICAL RECIPES
+# by Press, Flannery, Teukolsky, and Vetterling, p523.
+
+procedure id_mr_solve (x, y, npts, params, flags, np, nfit, mr, chisq)
+
+real x[npts] # X data array
+real y[npts] # Y data array
+int npts # Number of data points
+real params[np] # Parameter array
+int flags[np] # Flag array indexing parameters to fit
+int np # Number of parameters
+int nfit # Number of parameters to fit
+real mr # MR parameter
+real chisq # Chi square of fit
+
+int i
+real chisq1
+pointer new, a1, a2, delta1, delta2
+
+errchk id_mr_invert
+
+begin
+ # Allocate memory and initialize.
+ if (mr < 0.) {
+ call mfree (new, TY_REAL)
+ call mfree (a1, TY_REAL)
+ call mfree (a2, TY_REAL)
+ call mfree (delta1, TY_REAL)
+ call mfree (delta2, TY_REAL)
+
+ call malloc (new, np, TY_REAL)
+ call malloc (a1, nfit*nfit, TY_REAL)
+ call malloc (a2, nfit*nfit, TY_REAL)
+ call malloc (delta1, nfit, TY_REAL)
+ call malloc (delta2, nfit, TY_REAL)
+
+ call amovr (params, Memr[new], np)
+ call id_mr_eval (x, y, npts, Memr[new], flags, np, Memr[a2],
+ Memr[delta2], nfit, chisq)
+ mr = 0.001
+ }
+
+ # Restore last good fit and apply the Marquardt parameter.
+ call amovr (Memr[a2], Memr[a1], nfit * nfit)
+ call amovr (Memr[delta2], Memr[delta1], nfit)
+ do i = 1, nfit
+ Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr)
+
+ # Matrix solution.
+ call id_mr_invert (Memr[a1], Memr[delta1], nfit)
+
+ # Compute the new values and curvature matrix.
+ do i = 1, nfit
+ Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1]
+ call id_mr_eval (x, y, npts, Memr[new], flags, np, Memr[a1],
+ Memr[delta1], nfit, chisq1)
+
+ # Check if chisq has improved.
+ if (chisq1 < chisq) {
+ mr = max (EPSILONR, 0.1 * mr)
+ chisq = chisq1
+ call amovr (Memr[a1], Memr[a2], nfit * nfit)
+ call amovr (Memr[delta1], Memr[delta2], nfit)
+ call amovr (Memr[new], params, np)
+ } else
+ mr = 10. * mr
+
+ if (mr == 0.) {
+ call mfree (new, TY_REAL)
+ call mfree (a1, TY_REAL)
+ call mfree (a2, TY_REAL)
+ call mfree (delta1, TY_REAL)
+ call mfree (delta2, TY_REAL)
+ }
+end
+
+
+# ID_MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS.
+
+procedure id_mr_eval (x, y, npts, params, flags, np, a, delta, nfit, chisq)
+
+real x[npts] # X data array
+real y[npts] # Y data array
+int npts # Number of data points
+real params[np] # Parameter array
+int flags[np] # Flag array indexing parameters to fit
+int np # Number of parameters
+real a[nfit,nfit] # Curvature matrix
+real delta[nfit] # Delta array
+int nfit # Number of parameters to fit
+real chisq # Chi square of fit
+
+int i, j, k
+real ymod, dy, dydpj, dydpk
+pointer sp, dydp
+
+begin
+ call smark (sp)
+ call salloc (dydp, np, TY_REAL)
+
+ do j = 1, nfit {
+ do k = 1, j
+ a[j,k] = 0.
+ delta[j] = 0.
+ }
+
+ chisq = 0.
+ do i = 1, npts {
+ call id_derivs (x[i], params, ymod, Memr[dydp], np)
+ dy = y[i] - ymod
+ do j = 1, nfit {
+ dydpj = Memr[dydp+flags[j]-1]
+ delta[j] = delta[j] + dy * dydpj
+ do k = 1, j {
+ dydpk = Memr[dydp+flags[k]-1]
+ a[j,k] = a[j,k] + dydpj * dydpk
+ }
+ }
+ chisq = chisq + dy * dy
+ }
+
+ do j = 2, nfit
+ do k = 1, j-1
+ a[k,j] = a[j,k]
+
+ call sfree (sp)
+end
+
+
+# MR_INVERT -- Solve a set of linear equations using Householder transforms.
+
+procedure id_mr_invert (a, b, n)
+
+real a[n,n] # Input matrix and returned inverse
+real b[n] # Input RHS vector and returned solution
+int n # Dimension of input matrices
+
+int krank
+real rnorm
+pointer sp, h, g, ip
+
+begin
+ call smark (sp)
+ call salloc (h, n, TY_REAL)
+ call salloc (g, n, TY_REAL)
+ call salloc (ip, n, TY_INT)
+
+ call hfti (a, n, n, n, b, n, 1, 1E-10, krank, rnorm,
+ Memr[h], Memr[g], Memi[ip])
+
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/iddelete.x b/noao/rv/rvidlines/iddelete.x
new file mode 100644
index 00000000..cd96abb1
--- /dev/null
+++ b/noao/rv/rvidlines/iddelete.x
@@ -0,0 +1,26 @@
+include "identify.h"
+
+# ID_DELETE -- Delete a feature.
+
+procedure id_delete (id, feature)
+
+pointer id # ID pointer
+int feature # Feature to be deleted
+
+int i
+
+begin
+ call mfree (Memi[ID_LABEL(id)+feature-1], TY_CHAR)
+ do i = feature + 1, ID_NFEATURES(id) {
+ PIX(id,i-1) = PIX(id,i)
+ FIT(id,i-1) = FIT(id,i)
+ USER(id,i-1) = USER(id,i)
+ WTS(id,i-1) = WTS(id,i)
+ FWIDTH(id,i-1) = FWIDTH(id,i)
+ FTYPE(id,i-1) = FTYPE(id,i)
+ Memi[ID_LABEL(id)+i-2] = Memi[ID_LABEL(id)+i-1]
+ }
+ Memi[ID_LABEL(id)+ID_NFEATURES(id)-1] = NULL
+ ID_NFEATURES(id) = ID_NFEATURES(id) - 1
+ ID_NEWFEATURES(id) = YES
+end
diff --git a/noao/rv/rvidlines/iddofit.x b/noao/rv/rvidlines/iddofit.x
new file mode 100644
index 00000000..691e7350
--- /dev/null
+++ b/noao/rv/rvidlines/iddofit.x
@@ -0,0 +1,101 @@
+include "identify.h"
+
+# ID_DOFIT -- Fit a function to the features. Eliminate INDEF points.
+
+procedure id_dofit (id, interactive)
+
+pointer id # ID pointer
+int interactive # Interactive fit?
+
+int i, j, k, nfit, ic_geti()
+pointer gt1, sp, x, y, wts, rejpts, str, gt_init()
+
+begin
+ if (ID_NFEATURES(id) == 0) {
+ if (ID_CV(id) != NULL) {
+ call dcvfree (ID_CV(id))
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+ ID_NEWGRAPH(id) = YES
+ ID_NEWCV(id) = YES
+ }
+ return
+ }
+
+ call smark (sp)
+ call salloc (x, ID_NFEATURES(id), TY_DOUBLE)
+ call salloc (y, ID_NFEATURES(id), TY_DOUBLE)
+ call salloc (wts, ID_NFEATURES(id), TY_DOUBLE)
+
+ nfit = 0
+ do i = 1, ID_NFEATURES(id) {
+ if (IS_INDEFD (PIX(id,i)) || IS_INDEFD (USER(id,i)))
+ next
+ Memd[x+nfit] = PIX(id,i)
+ Memd[y+nfit] = USER(id,i)
+ Memd[wts+nfit] = 1.
+ nfit = nfit + 1
+ }
+
+ if (nfit > 1) {
+ if (interactive == YES) {
+ call salloc (str, SZ_LINE, TY_CHAR)
+ gt1 = gt_init()
+ call icg_fitd (ID_IC(id), ID_GP(id), "cursor", gt1, ID_CV(id),
+ Memd[x], Memd[y], Memd[wts], nfit)
+ call gt_free (gt1)
+ } else
+ call ic_fitd (ID_IC(id), ID_CV(id), Memd[x], Memd[y], Memd[wts],
+ nfit, YES, YES, YES, YES)
+
+ if (ic_geti (ID_IC(id), "nreject") > 0 &&
+ ic_geti (ID_IC(id), "nfit") == nfit)
+ rejpts = ic_geti (ID_IC(id), "rejpts")
+ else
+ rejpts = NULL
+
+ j = 0
+ k = 0
+ do i = 1, ID_NFEATURES(id) {
+ if (IS_INDEFD (PIX(id,i)) || IS_INDEFD (USER(id,i))) {
+ 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)
+ } else {
+ if (Memd[wts+k] != 0.) {
+ j = j + 1
+ PIX(id,j) = Memd[x+k]
+ FIT(id,j) = FIT(id,i)
+ USER(id,j) = Memd[y+k]
+ WTS(id,j) = Memd[wts+k]
+ if (rejpts != NULL)
+ if (Memi[rejpts+k] == YES)
+ WTS(id,j) = 0.
+ FWIDTH(id,j) = FWIDTH(id,i)
+ FTYPE(id,j) = FTYPE(id,i)
+ }
+ k = k + 1
+ }
+ }
+ ID_NFEATURES(id) = j
+
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ } else {
+ if (ID_CV(id) != NULL) {
+ call dcvfree (ID_CV(id))
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ }
+ }
+
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/iddoshift.x b/noao/rv/rvidlines/iddoshift.x
new file mode 100644
index 00000000..6e397455
--- /dev/null
+++ b/noao/rv/rvidlines/iddoshift.x
@@ -0,0 +1,42 @@
+include "identify.h"
+
+# ID_DOSHIFT -- Minimize residuals by constant shift.
+
+procedure id_doshift (id, interactive)
+
+pointer id # ID pointer
+int interactive # Called interactively?
+
+int i, j
+double shft, delta, rms, id_fitpt()
+
+begin
+ shft = 0.
+ rms = 0.
+ j = 0
+ for (i=1; i <= ID_NFEATURES(id); i = i + 1) {
+ if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.)
+ next
+ delta = USER(id,i) - id_fitpt (id, PIX(id,i))
+ shft = shft + delta
+ rms = rms + delta * delta
+ j = j + 1
+ }
+
+ if (j > 0) {
+ shft = shft / j
+ rms = rms / j
+ if (interactive == YES) {
+ call printf ("%s%s: Coordinate shift=%5f, rms=%5f, npts=%3d\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargd (shft)
+ call pargd (sqrt (rms - shft ** 2))
+ call pargi (j)
+ }
+ ID_SHIFT(id) = ID_SHIFT(id) + shft
+ ID_REDSHIFT(id) = 0.
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ }
+end
diff --git a/noao/rv/rvidlines/identify.h b/noao/rv/rvidlines/identify.h
new file mode 100644
index 00000000..ad34b767
--- /dev/null
+++ b/noao/rv/rvidlines/identify.h
@@ -0,0 +1,97 @@
+# Task parameters
+
+define LEN_IDSTRUCT 64 # Length ID structure
+
+define ID_TASK Memi[$1] # Task ID
+define ID_IMAGE Memi[$1+1] # Image (pointer)
+define ID_SECTION Memi[$1+2] # Section for 2D and 3D images (pointer)
+define ID_LINE Memi[$1+$2+2] # Image line or column [2]
+define ID_MAXLINE Memi[$1+$2+4] # Maximum line or column [2]
+define ID_AP Memi[$1+$2+6] # Aperture if appropriate [2]
+define ID_APS Memi[$1+9] # Array of apertures (pointer)
+define ID_NSUM Memi[$1+$2+10] # Number of lines to sum [2]
+define ID_MAXFEATURES Memi[$1+13] # Maximum number of features
+define ID_FTYPE Memi[$1+14] # Feature type
+define ID_MINSEP Memr[P2R($1+15)] # Minimum pixel separation
+define ID_MATCH Memr[P2R($1+16)] # Maximum matching separation
+define ID_FWIDTH Memr[P2R($1+17)] # Feature width in pixels
+define ID_CRADIUS Memr[P2R($1+18)] # Centering radius in pixels
+define ID_THRESHOLD Memr[P2R($1+19)] # Centering threshold
+define ID_ZWIDTH Memr[P2R($1+20)] # Zoom window width in fit units
+define ID_DATABASE Memi[$1+21] # Name of database (pointer)
+define ID_COORDLIST Memi[$1+22] # Name of coordinate list (pointer)
+define ID_LL Memi[$1+23] # Pointer to lines in coordinate list
+define ID_LABELS Memi[$1+24] # Type of feature labels
+define ID_LOGFILES Memi[$1+25] # List of logfiles
+
+# Common image data
+
+define ID_SHIFT Memd[P2D($1+26)]# Wavelength shift
+define ID_REDSHIFT Memd[P2D($1+28)]# Redshift of spectrum
+define ID_RMSRED Memd[P2D($1+30)]# Redshift of spectrum
+define ID_ZHELIO Memd[P2D($1+32)]# Heliocentric correction in redshift
+define ID_IMDATA Memi[$1+34] # Image data (pointer)
+define ID_PIXDATA Memi[$1+35] # Pixel coordinates (pointer)
+define ID_FITDATA Memi[$1+36] # Fit coordinates (pointer)
+define ID_NPTS Memi[$1+37] # Number of points
+
+# Features
+
+define ID_NFEATURES Memi[$1+38] # Number of features
+define ID_NALLOC Memi[$1+39] # Length of allocated feature arrays
+define ID_PIX Memi[$1+40] # Feature pixel coordinates (pointer)
+define ID_FIT Memi[$1+41] # Feature fit coordinates (pointer)
+define ID_USER Memi[$1+42] # Feature user coordinates (pointer)
+define ID_WTS Memi[$1+43] # Feature weights (pointer)
+define ID_FWIDTHS Memi[$1+44] # Feature width (pointer)
+define ID_FTYPES Memi[$1+45] # Feature type (pointer)
+define ID_LABEL Memi[$1+46] # Feature label (pointer)
+define ID_CURRENT Memi[$1+47] # Current feature
+
+# Pointers for other packages and to save data
+
+define ID_SH Memi[$1+48] # SHDR pointer
+define ID_LP Memi[$1+49] # Logical to physical transformation
+define ID_PL Memi[$1+50] # Physical to logical transformation
+define ID_IC Memi[$1+51] # ICFIT pointer
+define ID_CV Memi[$1+52] # Curfit pointer
+define ID_GP Memi[$1+53] # GIO pointer
+define ID_GT Memi[$1+54] # Gtools pointer
+define ID_ID Memi[$1+55] # Array of structure pointers (pointer)
+define ID_NID Memi[$1+56] # Number of saved structure
+define ID_DT Memi[$1+57] # Database pointer
+
+# Flags
+
+define ID_NEWFEATURES Memi[$1+58] # Has feature list changed?
+define ID_NEWCV Memi[$1+59] # Has fitting function changed?
+define ID_NEWGRAPH Memi[$1+60] # Has graph changed?
+define ID_NEWDBENTRY Memi[$1+61] # Has database entry changed?
+define ID_REFIT Memi[$1+62] # Refit feature data?
+define ID_GTYPE Memi[$1+63] # Graph type
+
+# End of structure ----------------------------------------------------------
+
+# Task ID
+define IDENTIFY 1 # Standard identify
+define RVIDLINES 2 # Line radial velocities
+
+define LABELS "|none|index|pixel|coord|user|both|"
+define FTYPES "|emission|absorption|gemission|gabsorption|"
+define EMISSION 1 # Emission feature (center1d)
+define ABSORPTION 2 # Absorption feature (center1d)
+define GEMISSION 3 # Emission feature (center1d)
+define GABSORPTION 4 # Absorption feature (center1d)
+
+define IMDATA Memr[ID_IMDATA($1)+$2-1]
+define PIXDATA Memd[ID_PIXDATA($1)+$2-1]
+define FITDATA Memd[ID_FITDATA($1)+$2-1]
+
+define PIX Memd[ID_PIX($1)+$2-1]
+define FIT Memd[ID_FIT($1)+$2-1]
+define USER Memd[ID_USER($1)+$2-1]
+define WTS Memd[ID_WTS($1)+$2-1]
+define FWIDTH Memr[ID_FWIDTHS($1)+$2-1]
+define FTYPE Memi[ID_FTYPES($1)+$2-1]
+
+define VLIGHT 2.997925e5 # Speed of light, Km/sec
diff --git a/noao/rv/rvidlines/identify.key b/noao/rv/rvidlines/identify.key
new file mode 100644
index 00000000..ae2f9339
--- /dev/null
+++ b/noao/rv/rvidlines/identify.key
@@ -0,0 +1,104 @@
+ IDENTIFY HELP
+
+
+STATUS LINE
+
+The status line gives the pixel position, fitted wavelength, user wavelength,
+wavelength residual, and optional line identification:
+
+ pixel fitted user residual [identification]
+
+
+CURSOR KEY SUMMARY
+
+? Help l Match list (fit) w Window graph
+a Affect all features m Mark feature x Crosscorrelate peaks
+b Deblend n Next feature y Find peaks
+c Center feature(s) o Go to line z Zoom graph
+d Delete feature(s) p Pan graph + Next feature
+f Fit positions q Quit - Previous feature
+g Fit zero point shift r Redraw graph . Nearest feature
+i Initialize s Shift feature I Interrupt
+j Preceding line t Reset position
+k Next line u Enter coordinate
+
+
+COLON COMMAND SUMMARY
+
+:show [file] :features [file] :coordlist [file]
+:cradius [value] :threshold [value] :database [file]
+:ftype [type] :fwidth [value] :image [image]
+:labels [type] :match [value] :maxfeatures [value]
+:minsep [value] :read [image [ap]] :write [image [ap]]
+:zwidth [value] :add [image [ap]] :vlog [file]
+
+
+CURSOR KEYS
+
+? Clear the screen and print menu of options
+a Apply next (c)enter or (d)elete operation to (a)ll features
+b Mark and de(b)lend features by Gaussian fitting
+c (C)enter the feature nearest the cursor
+d (D)elete the feature nearest the cursor
+f (F)it a function of pixel coordinate to the user coordinates
+g Fit a zero point shift to the user coordinates
+h Match coordinates in the coordinate list without modifying the fit
+i (I)nitialize (delete features and coordinate fit)
+j Go to the preceding image line or column in a 2D or multispec image
+k Go to the next image line or column in a 2D or multispec image
+l Match coordinates in the coordinate (l)ist and fit/refit the dispersion
+m (M)ark a new feature near the cursor and enter coordinate and label
+n Move the cursor or zoom to the (n)ext feature (same as +)
+o Go to the specified image line or column in a 2D or multispec image
+p (P)an to user defined window after (z)ooming on a feature
+q (Q)uit and continue with next image (also carriage return)
+r (R)edraw the graph
+s (S)hift the current feature to the position of the cursor
+t Reset the position of a feature without centering
+u Enter a new (u)ser coordinate and label for the current feature
+v Compute a redshift and velocity from the fitted and user coordinates
+w (W)indow the graph. Use '?' to window prompt for more help.
+x Crosscorrelate features with the data peaks and reregister
+y Automatically find "maxfeatures" strongest peaks and identify them
+z (Z)oom on the feature nearest the cursor
+. Move the cursor or zoom to the feature nearest the cursor
++ Move the cursor or zoom to the next feature
+- Move the cursor or zoom to the previous feature
+I Interrupt task and exit immediately. Database information is not saved.
+
+
+COLON COMMANDS
+
+The parameters are listed or set with the following commands which may be
+abbreviated. To list the value of a parameter type the command alone.
+
+:show file Show the values of all the parameters
+:features file Write feature list to file (default is STDOUT)
+:vlog file Write velocity information to file (default is STDOUT)
+
+:coordlist file Coordinate list file
+:cradius value Centering radius in pixels
+:threshold value Detection threshold for feature centering
+:database name Database for recording feature records
+:ftype value Feature type (emission or absorption)
+:fwidth value Feature width in pixels
+:image imagename Set a new image or show the current image
+:labels value Feature label type (none|index|pixel|coords|user|both)
+:match value Coordinate list matching distance
+:maxfeatures value Maximum number of features automatically found
+:minsep value Minimum separation allowed between features
+:read name ap Read a record from the database
+ (name and ap default to the current spectrum)
+:write name ap Write a record to the database
+ (name and ap default to the current spectrum)
+:add name ap Add features from the database
+ (name and ap default to the current spectrum)
+:zwidth value Zoom width in user units
+
+Labels:
+ none - No labels
+ index - Sequential numbers in order of increasing pixel position
+ pixel - Pixel coordinates
+ coords - User coordinates such as wavelength
+ user - User labels
+ both - Combination of coords and user
diff --git a/noao/rv/rvidlines/idfitdata.x b/noao/rv/rvidlines/idfitdata.x
new file mode 100644
index 00000000..48c36984
--- /dev/null
+++ b/noao/rv/rvidlines/idfitdata.x
@@ -0,0 +1,140 @@
+include <smw.h>
+include "identify.h"
+
+# ID_FITDATA -- Compute fit coordinates from pixel coordinates.
+
+procedure id_fitdata (id)
+
+pointer id # ID pointer
+int i
+
+begin
+ call mfree (ID_FITDATA(id), TY_DOUBLE)
+ call malloc (ID_FITDATA(id), ID_NPTS(id), TY_DOUBLE)
+
+ if (ID_CV(id) == NULL)
+ call achtrd (Memr[SX(ID_SH(id))], FITDATA(id,1), ID_NPTS(id))
+ else {
+ call dcvvector (ID_CV(id), PIXDATA(id,1), FITDATA(id,1),
+ ID_NPTS(id))
+ if (FITDATA(id,2) > FITDATA(id,1)) {
+ do i = 3, ID_NPTS(id)
+ if (FITDATA(id,i) < FITDATA(id,i-1))
+ call error (1, "Coordinate solution is not monotonic")
+ } else {
+ do i = 3, ID_NPTS(id)
+ if (FITDATA(id,i) > FITDATA(id,i-1))
+ call error (1, "Coordinate solution is not monotonic")
+ }
+ }
+ if (ID_SHIFT(id) != 0.)
+ call aaddkd (FITDATA(id,1), ID_SHIFT(id), FITDATA(id,1),ID_NPTS(id))
+end
+
+
+# ID_FITFEATURES -- Compute fit coordinates for features.
+
+procedure id_fitfeatures (id)
+
+pointer id # ID pointer
+int i
+
+double id_fitpt()
+
+begin
+ if (ID_NFEATURES(id) < 1)
+ return
+
+ if (ID_CV(id) == NULL)
+ do i = 1, ID_NFEATURES(id)
+ FIT(id,i) = id_fitpt (id, PIX(id,i))
+ else {
+ call dcvvector (ID_CV(id), PIX(id,1), FIT(id,1), ID_NFEATURES(id))
+ if (ID_SHIFT(id) != 0.)
+ call aaddkd (FIT(id,1), ID_SHIFT(id), FIT(id,1), ID_NFEATURES(id))
+ }
+end
+
+
+# ID_FITPT -- Compute fit coordinates from pixel coordinates.
+
+double procedure id_fitpt (id, pix)
+
+pointer id # ID pointer
+double pix # Pixel coordinate
+
+double fit
+
+double smw_c1trand(), shdr_lw(), dcveval()
+
+begin
+ if (ID_CV(id) == NULL) {
+ fit = smw_c1trand (ID_PL(id), pix)
+ fit = shdr_lw (ID_SH(id), fit)
+ } else
+ fit = dcveval (ID_CV(id), pix)
+ fit = fit + ID_SHIFT(id)
+
+ return (fit)
+end
+
+
+# FIT_TO_PIX -- Transform fit coordinate to pixel coordinate.
+
+define DXMIN .01
+
+double procedure fit_to_pix (id, fitcoord)
+
+pointer id # ID pointer
+double fitcoord # Fit coordinate to be transformed
+double pixcoord # Pixel coordinate returned
+
+int i, np1
+double dx
+
+double smw_c1trand(), id_fitpt()
+
+begin
+ np1 = NP1(ID_SH(id)) - 1
+ if (FITDATA(id,1) < FITDATA(id,ID_NPTS(id))) {
+ if ((fitcoord<FITDATA(id,1)) || (fitcoord>FITDATA(id,ID_NPTS(id))))
+ return (INDEFD)
+
+ for (i = 1; fitcoord > FITDATA(id,i); i = i + 1)
+ ;
+
+ if (FITDATA(id,i) == fitcoord)
+ return (PIXDATA(id,i))
+
+ pixcoord = smw_c1trand (ID_LP(id), double(i+np1-.5))
+ dx = smw_c1trand (ID_LP(id), double(i+np1+.5)) - pixcoord
+ while (dx > DXMIN) {
+ dx = dx / 2
+ if (id_fitpt (id, pixcoord) < fitcoord)
+ pixcoord = pixcoord + dx
+ else
+ pixcoord = pixcoord - dx
+ }
+ } else {
+ if ((fitcoord<FITDATA(id,ID_NPTS(id))) || (fitcoord>FITDATA(id,1)))
+ return (INDEFD)
+
+ for (i = 1; fitcoord < FITDATA(id,i); i = i + 1)
+ ;
+
+ if (FITDATA(id,i) == fitcoord)
+ return (PIXDATA(id,i))
+
+ pixcoord = smw_c1trand (ID_LP(id), double(i+np1-.5))
+ dx = smw_c1trand (ID_LP(id), double(i+np1+.5)) - pixcoord
+ while (dx > DXMIN) {
+ dx = dx / 2
+ if (id_fitpt (id, pixcoord) < fitcoord)
+ pixcoord = pixcoord - dx
+ else
+ pixcoord = pixcoord + dx
+ }
+ }
+
+ return (pixcoord)
+end
diff --git a/noao/rv/rvidlines/idfixx.x b/noao/rv/rvidlines/idfixx.x
new file mode 100644
index 00000000..ae74fcdf
--- /dev/null
+++ b/noao/rv/rvidlines/idfixx.x
@@ -0,0 +1,27 @@
+include <smw.h>
+
+# ID_FIXX - Adjust so that pixel indices are increasing.
+
+procedure id_fixx (sh, x1, x2, y1, y2, i1, i2)
+
+pointer sh
+real x1, x2, y1, y2
+int i1, i2
+
+double z, z1, z2, shdr_wl(), shdr_lw()
+
+begin
+ z1 = x1
+ z2 = x2
+ z1 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z1)))
+ z2 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z2)))
+ if (z1 > z2) {
+ z = y1; y1 = y2; y2 = z
+ z = z1; z1 = z2; z2 = z
+ }
+
+ x1 = shdr_lw (sh, z1)
+ x2 = shdr_lw (sh, z2)
+ i1 = nint (z1)
+ i2 = nint (z2)
+end
diff --git a/noao/rv/rvidlines/idgdata.x b/noao/rv/rvidlines/idgdata.x
new file mode 100644
index 00000000..caddce02
--- /dev/null
+++ b/noao/rv/rvidlines/idgdata.x
@@ -0,0 +1,74 @@
+include <imhdr.h>
+include <imio.h>
+include <pkg/gtools.h>
+include <smw.h>
+include "identify.h"
+
+define SZ_TITLE 320 # Size of long string for title.
+
+# ID_GDATA -- Get image data.
+
+procedure id_gdata (id)
+
+pointer id # ID pointer
+
+int i, np1
+double hjd
+pointer sp, str, im, mw, sh
+
+double smw_c1trand()
+errchk shdr_open, id_vhelio
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_TITLE, TY_CHAR)
+
+ sh = ID_SH(id)
+ im = IM(sh)
+ mw = MW(sh)
+
+ # If format is multispec then header info depends on line.
+ if (SMW_FORMAT(mw) == SMW_ES || SMW_FORMAT(mw) == SMW_MS)
+ ID_LINE(id,2) = 1
+ call shdr_open (im, mw, ID_LINE(id,1), ID_LINE(id,2),
+ INDEFI, SHDATA, sh)
+ ID_AP(id,1) = AP(sh)
+ ID_AP(id,2) = ID_LINE(id,2)
+ ID_NPTS(id) = SN(sh)
+ call id_dbsection (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ Memc[ID_SECTION(id)], SZ_FNAME)
+ call sprintf (Memc[str], SZ_TITLE, "%s %s%s\n%s")
+ if (ID_TASK(id) == IDENTIFY)
+ call pargstr ("identify")
+ else
+ call pargstr ("rvidlines")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargstr (TITLE(sh))
+ call gt_sets (ID_GT(id), GTTITLE, Memc[str])
+
+ # Free previous vectors and allocate new vectors.
+ call mfree (ID_PIXDATA(id), TY_DOUBLE)
+
+ call malloc (ID_PIXDATA(id), ID_NPTS(id), TY_DOUBLE)
+
+ # Set the physical coordinates.
+ np1 = NP1(sh) - 1
+ do i = 1, ID_NPTS(id)
+ PIXDATA(id,i) = smw_c1trand (ID_LP(id), double(i+np1))
+
+ # Set the image data
+ ID_IMDATA(id) = SY(sh)
+
+ # Set the heliocentric correction.
+ if (ID_TASK(id) == RVIDLINES) {
+ call id_vhelio (im, ID_ZHELIO(id), hjd, NULL)
+ ID_ZHELIO(id) = ID_ZHELIO(id) / VLIGHT
+ } else
+ ID_ZHELIO(id) = 0
+
+ ID_NEWGRAPH(id) = YES
+ ID_NEWCV(id) = YES
+
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/idgraph.x b/noao/rv/rvidlines/idgraph.x
new file mode 100644
index 00000000..c8d92b79
--- /dev/null
+++ b/noao/rv/rvidlines/idgraph.x
@@ -0,0 +1,168 @@
+include <gset.h>
+include <pkg/gtools.h>
+include <smw.h>
+include "identify.h"
+
+# ID_GRAPH -- Graph image vector in which features are to be identified.
+
+procedure id_graph (id, gtype)
+
+pointer id # ID pointer
+int gtype # Graph type
+
+begin
+ switch (gtype) {
+ case 1:
+ call id_graph1 (id)
+ case 2:
+ call id_graph2 (id)
+ default:
+ call id_graph1 (id)
+ }
+end
+
+
+procedure id_graph1 (id)
+
+pointer id # ID pointer
+
+int i, n
+real xmin, xmax, ymin, ymax, dy, xminz, xmaxz, id_zshiftr()
+pointer gp, sh, sp, x, y, str
+
+begin
+ gp = ID_GP(id)
+ sh = ID_SH(id)
+
+ call smark (sp)
+ call salloc (x, SN(sh), TY_REAL)
+ y = SY(sh)
+ n = SN(sh)
+
+ call achtdr (FITDATA(id,1), Memr[x], n)
+
+ call gclear (gp)
+ xmin = min (Memr[x], Memr[x+n-1])
+ xmax = max (Memr[x], Memr[x+n-1])
+ call alimr (Memr[y], n, ymin, ymax)
+ dy = ymax - ymin
+ call gswind (gp, xmin, xmax, ymin - .2 * dy, ymax + .2 * dy)
+ call gt_swind (gp, ID_GT(id))
+
+ if (ID_TASK(id) == RVIDLINES && ID_REDSHIFT(id) != 0D0) {
+ call salloc (str, SZ_LINE, TY_CHAR)
+ if (ID_ZHELIO(id) == 0D0)
+ call sprintf (Memc[str], SZ_LINE,
+ "Vobs = %.5g, Zobs = %.5g\n\n\n")
+ else
+ call sprintf (Memc[str], SZ_LINE,
+ "Vhelio = %.5g, Zhelio = %.5g\n\n\n")
+ call pargd ((ID_REDSHIFT(id)+ID_ZHELIO(id)) * VLIGHT)
+ call pargd (ID_REDSHIFT(id)+ID_ZHELIO(id))
+ call gt_sets (ID_GT(id), GTSUBTITLE, Memc[str])
+ call gseti (gp, G_XDRAWAXES, 1)
+ call gt_labax (gp, ID_GT(id))
+ call gt_sets (ID_GT(id), GTSUBTITLE, "")
+
+ call ggwind (gp, xmin, xmax, ymin, ymax)
+ xminz = id_zshiftr (id, xmin, 0)
+ xmaxz = id_zshiftr (id, xmax, 0)
+ call gswind (gp, xminz, xmaxz, ymin, ymax)
+ call gseti (gp, G_XDRAWAXES, 2)
+ call gseti (gp, G_YDRAWAXES, 0)
+ call glabax (gp, "", "", "")
+
+ call gswind (gp, xmin, xmax, ymin, ymax)
+ call gctran (gp, xmin, ymin, xmax, ymax, 1, 0)
+ call gctran (gp, xmax, ymax, xmin, ymin, 0, 1)
+ } else
+ call gt_labax (gp, ID_GT(id))
+
+ call gt_plot (gp, ID_GT(id), Memr[x], Memr[y], n)
+
+ do i = 1, ID_NFEATURES(id)
+ call id_mark (id, i)
+
+ call sfree (sp)
+end
+
+
+# ID_GRAPH2 -- Make review graph for current feature.
+
+procedure id_graph2 (id)
+
+pointer id # ID pointer
+
+int i, j, k, n
+real xmin, xmax, ymin, ymax, dy, xminz, xmaxz, id_zshiftr()
+pointer gp, sh, sp, x, y, str
+
+begin
+ gp = ID_GP(id)
+ sh = ID_SH(id)
+
+ call smark (sp)
+ call salloc (x, SN(sh), TY_REAL)
+ y = SY(sh)
+ n = SN(sh)
+
+ call achtdr (FITDATA(id,1), Memr[x], n)
+
+ xmin = real (FIT(id,ID_CURRENT(id))) - ID_ZWIDTH(id) / 2.
+ xmax = real (FIT(id,ID_CURRENT(id))) + ID_ZWIDTH(id) / 2.
+
+ i = 0
+ do k = 1, n {
+ if ((Memr[x+k-1] < xmin) || (Memr[x+k-1] > xmax))
+ next
+ if (i == 0)
+ i = k
+ j = k
+ }
+ k = j - i + 1
+
+ call alimr (Memr[y+i-1], k, ymin, ymax)
+ dy = ymax - ymin
+
+ call gclear (gp)
+ call gswind (gp, xmin, xmax, ymin - .2 * dy, ymax + .2 * dy)
+# if (ID_GT(id) != NULL) {
+# call gseti (gp, G_XTRAN, GT_XTRAN(ID_GT(id)))
+# call gseti (gp, G_YTRAN, GT_YTRAN(ID_GT(id)))
+# }
+ if (ID_TASK(id) == RVIDLINES && ID_REDSHIFT(id) != 0D0) {
+ call salloc (str, SZ_LINE, TY_CHAR)
+ if (ID_ZHELIO(id) == 0D0)
+ call sprintf (Memc[str], SZ_LINE,
+ "Vobs = %.5g, Zobs = %.5g\n\n\n")
+ else
+ call sprintf (Memc[str], SZ_LINE,
+ "Vhelio = %.5g, Zhelio = %.5g\n\n\n")
+ call pargd ((ID_REDSHIFT(id)+ID_ZHELIO(id)) * VLIGHT)
+ call pargd (ID_REDSHIFT(id)+ID_ZHELIO(id))
+ call gt_sets (ID_GT(id), GTSUBTITLE, Memc[str])
+ call gseti (gp, G_XDRAWAXES, 1)
+ call gt_labax (gp, ID_GT(id))
+ call gt_sets (ID_GT(id), GTSUBTITLE, "")
+
+ call ggwind (gp, xmin, xmax, ymin, ymax)
+ xminz = id_zshiftr (id, xmin, 0)
+ xmaxz = id_zshiftr (id, xmax, 0)
+ call gswind (gp, xminz, xmaxz, ymin, ymax)
+ call gseti (gp, G_XDRAWAXES, 2)
+ call gseti (gp, G_YDRAWAXES, 0)
+ call glabax (gp, "", "", "")
+
+ call gswind (gp, xmin, xmax, ymin, ymax)
+ call gctran (gp, xmin, ymin, xmax, ymax, 1, 0)
+ call gctran (gp, xmax, ymax, xmin, ymin, 0, 1)
+ } else
+ call gt_labax (gp, ID_GT(id))
+
+ call gt_plot (gp, ID_GT(id), Memr[x], Memr[y], n)
+
+ do i = 1, ID_NFEATURES(id)
+ call id_mark (id, i)
+
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/ididentify.x b/noao/rv/rvidlines/ididentify.x
new file mode 100644
index 00000000..9874865e
--- /dev/null
+++ b/noao/rv/rvidlines/ididentify.x
@@ -0,0 +1,795 @@
+include <error.h>
+include <imhdr.h>
+include <gset.h>
+include <smw.h>
+include "identify.h"
+
+define HELP "noao$onedspec/identify/identify.key"
+define RVHELP "noao$rv/rvidlines/rvidlines.key"
+define ICFITHELP "noao$lib/scr/idicgfit.key"
+define PROMPT "identify options"
+
+define PAN 1 # Pan graph
+define ZOOM 2 # Zoom graph
+
+# ID_IDENTIFY -- Identify features in an image.
+# This is the main interactive loop.
+
+procedure id_identify (id)
+
+pointer id # ID pointer
+
+real wx, wy, wx2, wy2
+int wcs, key
+char cmd[SZ_LINE]
+
+char newimage[SZ_FNAME]
+int i, j, last, all, prfeature, nfeatures1, npeaks, ftype, newline[2]
+bool answer
+double pix, fit, user, shift, pix_shift, z_shift, xg[10]
+pointer peaks, label
+
+bool clgetb()
+pointer gopen()
+int clgcur(), scan(), nscan(), find_peaks(), errcode(), id_gid(), nowhite()
+double id_center(), fit_to_pix(), id_fitpt()
+double id_shift(), id_rms(), id_zshiftd(), id_zval()
+errchk id_gdata(), id_graph(), id_dbread(), xt_mk1d(), id_log()
+
+define newim_ 10
+define newkey_ 20
+define beep_ 99
+
+begin
+newim_
+ # Open the image and return if there is an error.
+ iferr (call id_map (id)) {
+ call erract (EA_WARN)
+ return
+ }
+
+ # Get the image data and return if there is an error.
+ iferr (call id_gdata (id)) {
+ call erract (EA_WARN)
+ return
+ }
+
+ # Get the database entry for the image if it exists.
+ iferr {
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO, YES)
+ ID_NEWDBENTRY(id) = NO
+ } then
+ if ((ID_NFEATURES(id) > 0) || (ID_CV(id) != NULL))
+ ID_NEWDBENTRY(id) = YES
+
+ # Set the coordinate information.
+ iferr (call id_fitdata (id))
+ ;
+
+ # Set fitting limits.
+ call ic_putr (ID_IC(id), "xmin", real (PIXDATA(id,1)))
+ call ic_putr (ID_IC(id), "xmax", real (PIXDATA(id,ID_NPTS(id))))
+ call ic_pstr (ID_IC(id), "help", ICFITHELP)
+
+ # Open graphics.
+ call clgstr ("graphics", newimage, SZ_FNAME)
+ ID_GP(id) = gopen (newimage, NEW_FILE, STDGRAPH)
+
+ # Initialize.
+ ID_GTYPE(id) = PAN
+ all = 0
+ last = ID_CURRENT(id)
+ newimage[1] = EOS
+ newline[1] = ID_LINE(id,1)
+ newline[2] = ID_LINE(id,2)
+ ID_REFIT(id) = NO
+ ID_NEWFEATURES(id) = NO
+ ID_NEWCV(id) = NO
+ wy = INDEF
+ key = 'r'
+
+ repeat {
+ prfeature = YES
+ if (all != 0)
+ all = mod (all + 1, 3)
+
+ switch (key) {
+ case '?': # Print help
+ if (ID_TASK(id) == IDENTIFY)
+ call gpagefile (ID_GP(id), HELP, PROMPT)
+ else
+ call gpagefile (ID_GP(id), RVHELP, PROMPT)
+ case ':': # Process colon commands
+ if (cmd[1] == '/')
+ call gt_colon (cmd, ID_GP(id), ID_GT(id), ID_NEWGRAPH(id))
+ else
+ call id_colon (id, cmd, newimage, prfeature)
+ case ' ': # Go to current feature
+ case '.': # Go to nearest feature
+ if (ID_NFEATURES(id) == 0)
+ goto beep_
+ call id_nearest (id, double (wx))
+ case '-': # Go to previous feature
+ if (ID_CURRENT(id) == 1)
+ goto beep_
+ ID_CURRENT(id) = ID_CURRENT(id) - 1
+ case '+', 'n': # Go to next feature
+ if (ID_CURRENT(id) == ID_NFEATURES(id))
+ goto beep_
+ ID_CURRENT(id) = ID_CURRENT(id) + 1
+ case 'a': # Set all flag for next key
+ all = 1
+ case 'b': # Mark blended features
+ i = 1
+ fit = wx
+ xg[i] = fit_to_pix (id, fit)
+ call printf ("mark other components (exit with 'q'):")
+ while (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)!=EOF) {
+ if (key == 'q')
+ break
+ i = i + 1
+ fit = wx
+ xg[i] = fit_to_pix (id, fit)
+ if (i == 10)
+ break
+ }
+
+ call printf ("mark two background points: ")
+ j = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+ j = clgcur ("cursor", wx2, wy2, wcs, key, cmd, SZ_LINE)
+ wx = fit_to_pix (id, double (wx))
+ wx2 = fit_to_pix (id, double (wx2))
+
+ switch (ID_FTYPE(id)) {
+ case EMISSION:
+ ftype = GEMISSION
+ case ABSORPTION:
+ ftype = GABSORPTION
+ default:
+ ftype = ID_FTYPE(id)
+ }
+ iferr (call id_gcenter (id, xg, i, wx, wy, wx2, wy2,
+ ID_FWIDTH(id), ftype, YES)) {
+ call erract (EA_WARN)
+ prfeature = NO
+ goto newkey_
+ }
+
+ do j = 1, i {
+ pix = xg[j]
+ fit = id_fitpt (id, pix)
+ user = fit
+ call id_newfeature (id, pix, fit, user, 1.0D0,
+ ID_FWIDTH(id), ftype, NULL)
+ USER(id,ID_CURRENT(id)) = INDEFD
+ call id_match (id, FIT(id,ID_CURRENT(id)),
+ USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ call id_mark (id, ID_CURRENT(id))
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("%10.8g ")
+ call pargd (
+ id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ fit = id_zshiftd (id, user, 1)
+ call id_match (id, fit, USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd,
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ }
+ case 'c': # Recenter features
+ if (all != 0) {
+ for (i = 1; i <= ID_NFEATURES(id); i = i + 1) {
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, i)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ FWIDTH(id,i) = ID_FWIDTH(id)
+ PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i),
+ FTYPE(id,i), NO)
+ if (!IS_INDEFD (PIX(id,i))) {
+ FIT(id,i) = id_fitpt (id, PIX(id,i))
+ call id_mark (id, i)
+ } else {
+ call id_delete (id, i)
+ i = i - 1
+ }
+ }
+ ID_NEWFEATURES(id) = YES
+ } else {
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call id_nearest (id, double (wx))
+ pix = PIX(id,ID_CURRENT(id))
+ pix = id_center (id, pix, 1, ID_FWIDTH(id),
+ FTYPE(id,ID_CURRENT(id)), NO)
+ if (!IS_INDEFD (pix)) {
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ PIX(id,ID_CURRENT(id)) = pix
+ FWIDTH(id,ID_CURRENT(id)) = ID_FWIDTH(id)
+ FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_mark (id, ID_CURRENT(id))
+ ID_NEWFEATURES(id) = YES
+ } else {
+ call printf ("Centering failed\n")
+ prfeature = NO
+ }
+ }
+ case 'd': # Delete features
+ if (all != 0) {
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ ID_NEWFEATURES(id) = YES
+ ID_NEWGRAPH(id) = YES
+ } else {
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call id_nearest (id, double (wx))
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_delete (id, ID_CURRENT(id))
+ ID_CURRENT(id) = min (ID_NFEATURES(id), ID_CURRENT(id))
+ last = 0
+ }
+ case 'f': # Fit dispersion function
+ if (ID_TASK(id) == IDENTIFY) {
+ call id_dofit (id, YES)
+ } else {
+ call id_velocity (id, YES)
+ prfeature = NO
+ }
+ case 'g': # Fit shift
+ if (ID_TASK(id) == RVIDLINES)
+ goto beep_
+
+ call id_doshift (id, YES)
+ prfeature = NO
+ case 'i': # Initialize
+ call dcvfree (ID_CV(id))
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+ ID_NEWCV(id) = YES
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ ID_NEWFEATURES(id) = YES
+ ID_NEWGRAPH(id) = YES
+ case 'j': # Go to previous line
+ newline[1] = ID_LINE(id,1) - ID_NSUM(id,1)
+ if (newline[1] < 1) {
+ newline[1] = newline[1] + ID_MAXLINE(id,1)
+ newline[2] = ID_LINE(id,2) - ID_NSUM(id,2)
+ if (newline[2] < 1)
+ newline[2] = newline[2] + ID_MAXLINE(id,2)
+ }
+ case 'k': # Go to next line
+ newline[1] = ID_LINE(id,1) + ID_NSUM(id,1)
+ if (newline[1] > ID_MAXLINE(id,1)) {
+ newline[1] = newline[1] - ID_MAXLINE(id,1)
+ newline[2] = ID_LINE(id,2) + ID_NSUM(id,2)
+ if (newline[2] > ID_MAXLINE(id,2))
+ newline[2] = newline[2] - ID_MAXLINE(id,2)
+ }
+ case 'l': # Find features from line list
+ if (ID_TASK(id) == IDENTIFY) {
+ if (ID_NFEATURES(id) >= 2)
+ call id_dofit (id, NO)
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata(id))
+ ;
+ call id_fitfeatures(id)
+ ID_NEWCV(id) = NO
+ }
+ call id_linelist (id)
+ if (ID_NEWFEATURES(id) == YES)
+ ID_REFIT(id) = YES
+ } else {
+ call id_velocity (id, NO)
+ call id_linelist (id)
+ if (ID_NEWFEATURES(id) == YES) {
+ call id_velocity (id, NO)
+ ID_NEWGRAPH(id) = YES
+ }
+ }
+ case 'm': # Mark new feature
+ fit = wx
+ pix = fit_to_pix (id, fit)
+ pix = id_center (id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), YES)
+ if (IS_INDEFD (pix))
+ goto beep_
+ fit = id_fitpt (id, pix)
+ user = fit
+ call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id),
+ ID_FTYPE(id), NULL)
+ USER(id,ID_CURRENT(id)) = INDEFD
+ call id_match (id, FIT(id,ID_CURRENT(id)),
+ USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ call id_mark (id, ID_CURRENT(id))
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("%10.8g ")
+ call pargd (
+ id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ fit = id_zshiftd (id, user, 1)
+ call id_match (id, fit, USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ case 'o': # Go to a specified line
+ call printf ("Line/Column/Band (%d %d): ")
+ call pargi (ID_LINE(id,1))
+ call pargi (ID_LINE(id,2))
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargi (j)
+ if (nscan() == 1) {
+ if (j < 1 || j > ID_MAXLINE(id,1))
+ goto beep_
+ newline[1] = j
+ call gargi (j)
+ if (nscan() == 2) {
+ if (j < 1 || j > ID_MAXLINE(id,2))
+ goto beep_
+ newline[2] = j
+ }
+ }
+ }
+ case 'p': # Switch to pan mode
+ if (ID_GTYPE(id) != PAN) {
+ ID_GTYPE(id) = PAN
+ ID_NEWGRAPH(id) = YES
+ }
+ case 'q': # Exit loop
+ break
+ case 'r': # Redraw the graph
+ ID_NEWGRAPH(id) = YES
+ case 's', 'x': # Shift or correlate features
+ if (ID_TASK(id) == RVIDLINES)
+ goto beep_
+
+ # Get coordinate shift.
+ switch (key) {
+ case 's':
+ call printf ("User coordinate (%10.8g): ")
+ call pargr (wx)
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ if (nscan() == 1)
+ shift = wx - user
+ } else
+ shift = 0.
+ case 'x':
+ if (ID_NFEATURES(id) > 5)
+ shift = id_shift (id)
+ else
+ goto beep_
+ }
+
+ ID_NEWFEATURES(id) = YES
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ prfeature = NO
+
+ if (ID_NFEATURES(id) < 1) {
+ call printf ("User coordinate shift=%5f\n")
+ call pargd (shift)
+ ID_SHIFT(id) = ID_SHIFT(id) + shift
+ goto newkey_
+ }
+
+ # Recenter features.
+ pix_shift = 0.
+ z_shift = 0.
+ nfeatures1 = ID_NFEATURES(id)
+
+ j = 0.
+ do i = 1, ID_NFEATURES(id) {
+ pix = fit_to_pix (id, FIT(id,i) + shift)
+ pix = id_center (id, pix, 1, FWIDTH(id,i), FTYPE(id,i), NO)
+ if (IS_INDEFD (pix)) {
+ if (ID_CURRENT(id) == i)
+ ID_CURRENT(id) = i + 1
+ next
+ }
+ fit = id_fitpt (id, pix)
+
+ pix_shift = pix_shift + pix - PIX(id,i)
+ if (FIT(id,i) != 0.)
+ z_shift = z_shift + id_zval (id, fit, FIT(id,i))
+
+ j = j + 1
+ PIX(id,j) = pix
+ 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)
+ if (ID_CURRENT(id) == i)
+ ID_CURRENT(id) = j
+ }
+ if (j != ID_NFEATURES(id)) {
+ ID_NFEATURES(id) = j
+ ID_CURRENT(id) = min (ID_CURRENT(id), ID_NFEATURES(id))
+ }
+
+ if (ID_NFEATURES(id) < 1) {
+ call printf ("User coordinate shift=%5f")
+ call pargd (shift)
+ call printf (", No features found during recentering\n")
+ ID_SHIFT(id) = ID_SHIFT(id) + shift
+ goto newkey_
+ }
+
+ # Adjust shift.
+ pix = ID_SHIFT(id)
+ call id_doshift (id, NO)
+ call id_fitfeatures (id)
+
+ # Print results.
+ call printf ("Recentered=%d/%d")
+ call pargi (ID_NFEATURES(id))
+ call pargi (nfeatures1)
+ call printf (
+ ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g\n")
+ call pargd (pix_shift / ID_NFEATURES(id))
+ call pargd (pix - ID_SHIFT(id))
+ call pargd (z_shift / ID_NFEATURES(id))
+ call pargd (id_rms(id))
+ case 't': # Move the current feature
+ if (ID_CURRENT(id) < 1)
+ goto beep_
+ pix = fit_to_pix (id, double (wx))
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ PIX(id,ID_CURRENT(id)) = pix
+ FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_mark (id, ID_CURRENT(id))
+ ID_NEWFEATURES(id) = YES
+ case 'u': # Set user coordinate
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("%10.8g ")
+ call pargd (id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ ID_NEWFEATURES(id) = YES
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ call printf ("Weight (%g): ")
+ call pargd (WTS(id,ID_CURRENT(id)))
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ if (nscan() > 0)
+ WTS(id,ID_CURRENT(id)) = user
+ }
+ case 'w': # Window graph
+ call gt_window (ID_GT(id), ID_GP(id), "cursor", ID_NEWGRAPH(id))
+ case 'y': # Find peaks
+ call alimr (IMDATA(id,1), ID_NPTS(id), wx, wy)
+ call malloc (peaks, ID_NPTS(id), TY_REAL)
+ if (ID_FTYPE(id) == ABSORPTION || ID_FTYPE(id) == GABSORPTION)
+ npeaks = find_peaks (IMDATA(id,1), Memr[peaks],
+ ID_NPTS(id), -1., int (ID_MINSEP(id)), 0,
+ ID_MAXFEATURES(id), wy, false)
+ else
+ npeaks = find_peaks (IMDATA(id,1), Memr[peaks],
+ ID_NPTS(id), 0., int (ID_MINSEP(id)), 0,
+ ID_MAXFEATURES(id), wx, false)
+ for (j = 1; j <= ID_NFEATURES(id); j = j + 1) {
+ for (i = 1; i <= npeaks; i = i + 1) {
+ if (!IS_INDEF (Memr[peaks+i-1])) {
+ pix = Memr[peaks+i-1]
+ if (abs (pix - PIX(id,j)) < ID_MINSEP(id))
+ Memr[peaks+i-1] = INDEF
+ }
+ }
+ }
+ for (i = 1; i <= npeaks; i = i + 1) {
+ if (IS_INDEF(Memr[peaks+i-1]))
+ next
+ pix = Memr[peaks+i-1]
+ pix = id_center(id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), NO)
+ if (IS_INDEFD (pix))
+ next
+ fit = id_fitpt (id, pix)
+ user = INDEFD
+ call id_match (id, fit, user, label, ID_MATCH(id))
+ call id_newfeature (id, pix, fit, user, 1.0D0,
+ ID_FWIDTH(id), ID_FTYPE(id), label)
+ call id_mark (id, ID_CURRENT(id))
+ }
+ call mfree (peaks, TY_REAL)
+ case 'z': # Go to zoom mode
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ if (ID_GTYPE(id) == PAN)
+ ID_NEWGRAPH(id) = YES
+ ID_GTYPE(id) = ZOOM
+ call id_nearest (id, double (wx))
+ case 'I':
+ call fatal (0, "Interrupt")
+ default:
+beep_ call printf ("\007")
+ }
+
+newkey_
+ # Set update flag if anything has changed.
+ if ((ID_NEWFEATURES(id) == YES) || (ID_NEWCV(id) == YES))
+ ID_NEWDBENTRY(id) = YES
+
+ # If a new image exit loop, update database, and start over.
+ if (newimage[1] != EOS)
+ break
+
+ # If a new line, save features and set new line.
+ if (newline[1] != ID_LINE(id,1) || newline[2] != ID_LINE(id,2)) {
+ call id_saveid (id, ID_LINE(id,1))
+ ID_LINE(id,1) = newline[1]
+ ID_LINE(id,2) = newline[2]
+ call id_gdata (id)
+ if (id_gid (id, newline) == EOF) {
+ iferr {
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ NO, NO)
+ ID_NEWDBENTRY(id) = NO
+ ID_NEWFEATURES(id) = NO
+ } then
+ if ((ID_NFEATURES(id) > 0) || (ID_CV(id) != NULL))
+ ID_NEWDBENTRY(id) = YES
+ }
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ wy = INDEF
+ }
+
+ # Refit dispersion function
+ if (ID_REFIT(id) == YES) {
+ call id_dofit (id, NO)
+ ID_REFIT(id) = NO
+ }
+
+ # If there is a new dispersion solution evaluate the coordinates
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata (id))
+ ;
+ call id_fitfeatures (id)
+ ID_NEWCV(id) = NO
+ }
+
+ # Draw new graph in zoom mode if current feature has changed.
+ if ((ID_GTYPE(id) == ZOOM) && (last != ID_CURRENT(id)))
+ ID_NEWGRAPH(id) = YES
+
+ # Draw new graph.
+ if (ID_NEWGRAPH(id) == YES) {
+ call id_graph (id, ID_GTYPE(id))
+ ID_NEWGRAPH(id) = NO
+ }
+
+ # Set cursor and print status of current feature (unless canceled).
+ if (ID_CURRENT(id) > 0) {
+ if (IS_INDEF (wy)) {
+ i = max (1, min (ID_NPTS(id), int (PIX(id,ID_CURRENT(id)))))
+ wy = IMDATA(id,i)
+ }
+
+ call gscur (ID_GP(id), real (FIT(id,ID_CURRENT(id))), wy)
+ if (errcode() == OK && prfeature == YES) {
+ i = ID_CURRENT(id)
+ pix = PIX(id,i)
+ fit = FIT(id,i)
+ user = USER(id,i)
+ if (ID_TASK(id) == IDENTIFY) {
+ if (IS_INDEFD(user))
+ shift = INDEF
+ else
+ shift = fit - user
+ call printf ("%10.2f %10.8g %10.8g %10.3g %s")
+ call pargd (pix)
+ call pargd (fit)
+ call pargd (user)
+ call pargd (shift)
+ label = Memi[ID_LABEL(id)+i-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ } else {
+ if (IS_INDEFD(user))
+ shift = INDEF
+ else {
+ shift = id_zval (id, fit, user) - ID_REDSHIFT(id)
+ if (abs (shift) < 0.01)
+ shift = shift * VLIGHT
+ }
+ call printf ("%10.2f %10.8g %10.8g %10.8g %10.3g %s")
+ call pargd (pix)
+ call pargd (fit)
+ call pargd (id_zshiftd (id, fit, 0))
+ call pargd (user)
+ if (IS_INDEFD(user))
+ call pargd (INDEFD)
+ else
+ call pargd (shift)
+ label = Memi[ID_LABEL(id)+i-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ }
+ }
+ }
+
+ # Print delayed error message
+ if (errcode() != OK)
+ call erract (EA_WARN)
+
+ last = ID_CURRENT(id)
+ } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+ call gclose (ID_GP(id))
+
+ # Write velocity logfile.
+ if (ID_TASK(id) == RVIDLINES) {
+ call clgstr ("logfile", cmd, SZ_LINE)
+ if (nowhite (cmd, cmd, SZ_LINE) == 0)
+ answer = false
+ else
+ answer = true
+ if (answer && !clgetb ("autowrite")) {
+ call printf ("Write velocity data to the logfile (yes)? ")
+ call flush (STDOUT)
+ if (scan() != EOF)
+ call gargb (answer)
+ }
+ if (answer) {
+ iferr {
+ if (ID_NID(id) == 0)
+ call id_log (id, cmd, NULL)
+ else {
+ call id_saveid (id, ID_LINE(id,1))
+ for (i = 0; i < ID_NID(id); i = i + 1) {
+ j = Memi[ID_ID(id)+i]
+ j =id_gid (id, ID_LINE(j,1))
+ call id_dbsection (id, Memc[ID_IMAGE(id)],
+ ID_AP(id,1), Memc[ID_SECTION(id)],
+ SZ_FNAME)
+ call id_log (id, cmd, NULL)
+ }
+ }
+ } then
+ call erract (EA_WARN)
+ }
+ }
+
+ # Warn user that feature data is newer than database entry.
+ if (ID_NEWDBENTRY(id) == YES)
+ answer = true
+ else {
+ answer = false
+ for (i = 0; i < ID_NID(id); i = i + 1) {
+ if (ID_NEWDBENTRY(Memi[ID_ID(id)+i]) == YES) {
+ answer = true
+ break
+ }
+ }
+ }
+ if (answer) {
+ if (!clgetb ("autowrite")) {
+ call printf ("Write feature data to the database (yes)? ")
+ call flush (STDOUT)
+ if (scan() != EOF)
+ call gargb (answer)
+ }
+ if (answer) {
+ newline[1] = ID_LINE(id,1)
+ newline[2] = ID_LINE(id,2)
+ if (ID_NEWDBENTRY(id) == YES)
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ for (i = 0; i < ID_NID(id); i = i + 1) {
+ j = Memi[ID_ID(id)+i]
+ if (ID_NEWDBENTRY(j) == YES &&
+ (ID_LINE(j,1)!=newline[1]||ID_LINE(j,2)!=newline[2])) {
+ j =id_gid (id, ID_LINE(j,1))
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ NO)
+ }
+ }
+ }
+ }
+
+ call flush (STDOUT)
+
+ # Free image data.
+ call mfree (ID_PIXDATA(id), TY_DOUBLE)
+ call mfree (ID_FITDATA(id), TY_DOUBLE)
+ call id_free1 (id)
+
+ call smw_close (MW(ID_SH(id)))
+ call imunmap (IM(ID_SH(id)))
+ call shdr_close (ID_SH(id))
+
+ # If a new image was requested with colon command start over.
+ if (newimage[1] != EOS) {
+ call strcpy (newimage, Memc[ID_IMAGE(id)], SZ_FNAME)
+ goto newim_
+ }
+end
diff --git a/noao/rv/rvidlines/idinit.x b/noao/rv/rvidlines/idinit.x
new file mode 100644
index 00000000..b17f94a1
--- /dev/null
+++ b/noao/rv/rvidlines/idinit.x
@@ -0,0 +1,352 @@
+include <gset.h>
+include <math/curfit.h>
+include "identify.h"
+
+# ID_INIT -- Allocate identify structure
+
+procedure id_init (id, taskid)
+
+int taskid #I Task ID
+pointer id #O ID pointer
+
+begin
+ call calloc (id, LEN_IDSTRUCT, TY_STRUCT)
+
+ ID_TASK(id) = taskid
+ ID_NALLOC(id) = 20
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ ID_NID(id) = 0
+ ID_DT(id) = NULL
+
+ call malloc (ID_IMAGE(id), SZ_FNAME, TY_CHAR)
+ call malloc (ID_SECTION(id), SZ_FNAME, TY_CHAR)
+ call malloc (ID_DATABASE(id), SZ_FNAME, TY_CHAR)
+ call malloc (ID_COORDLIST(id), SZ_FNAME, TY_CHAR)
+
+ call malloc (ID_PIX(id), ID_NALLOC(id), TY_DOUBLE)
+ call malloc (ID_FIT(id), ID_NALLOC(id), TY_DOUBLE)
+ call malloc (ID_USER(id), ID_NALLOC(id), TY_DOUBLE)
+ call malloc (ID_WTS(id), ID_NALLOC(id), TY_DOUBLE)
+ call malloc (ID_FWIDTHS(id), ID_NALLOC(id), TY_REAL)
+ call malloc (ID_FTYPES(id), ID_NALLOC(id), TY_INT)
+ call calloc (ID_LABEL(id), ID_NALLOC(id), TY_POINTER)
+end
+
+
+# ID_FREE -- Free identify structure.
+
+procedure id_free (id)
+
+pointer id # ID pointer
+
+int i
+pointer ptr
+
+begin
+ call mfree (ID_IMAGE(id), TY_CHAR)
+ call mfree (ID_SECTION(id), TY_CHAR)
+ call mfree (ID_DATABASE(id), TY_CHAR)
+ call mfree (ID_COORDLIST(id), TY_CHAR)
+ call mfree (ID_APS(id), TY_INT)
+
+ ptr = ID_LABEL(id)
+ do i = 1, ID_NALLOC(id) {
+ call mfree (Memi[ptr], TY_CHAR)
+ ptr = ptr + 1
+ }
+
+ call mfree (ID_PIX(id), TY_DOUBLE)
+ call mfree (ID_FIT(id), TY_DOUBLE)
+ call mfree (ID_USER(id), TY_DOUBLE)
+ call mfree (ID_WTS(id), TY_DOUBLE)
+ call mfree (ID_FWIDTHS(id), TY_REAL)
+ call mfree (ID_FTYPES(id), TY_INT)
+ call mfree (ID_LABEL(id), TY_POINTER)
+
+ if (ID_DT(id) != NULL)
+ call dtunmap (ID_DT(id))
+ call id_free1 (id)
+ call id_unmapll (id)
+ call gt_free (ID_GT(id))
+ call dcvfree (ID_CV(id))
+ call ic_closed (ID_IC(id))
+
+ call mfree (id, TY_STRUCT)
+end
+
+
+# ID_FREE1 -- Free saved identify structures.
+
+procedure id_free1 (id)
+
+pointer id # ID pointer
+
+int i, j
+pointer id1, ptr
+
+begin
+ for (i = 0; i < ID_NID(id); i = i + 1) {
+ id1 = Memi[ID_ID(id)+i]
+
+ ptr = ID_LABEL(id1)
+ do j = 1, ID_NALLOC(id1) {
+ call mfree (Memi[ptr], TY_CHAR)
+ ptr = ptr + 1
+ }
+
+ call mfree (ID_PIX(id1), TY_DOUBLE)
+ call mfree (ID_FIT(id1), TY_DOUBLE)
+ call mfree (ID_USER(id1), TY_DOUBLE)
+ call mfree (ID_WTS(id1), TY_DOUBLE)
+ call mfree (ID_FWIDTHS(id1), TY_REAL)
+ call mfree (ID_FTYPES(id1), TY_INT)
+ call mfree (ID_LABEL(id1), TY_POINTER)
+ if (ID_CV(id1) != NULL)
+ call dcvfree (ID_CV(id1))
+ if (ID_IC(id1) != NULL)
+ call ic_closed (ID_IC(id1))
+ call shdr_close (ID_SH(id1))
+ call mfree (id1, TY_STRUCT)
+ }
+ call mfree (ID_ID(id), TY_POINTER)
+
+ ID_NID(id) = 0
+end
+
+
+# ID_SAVEID -- Save identify information by line.
+
+procedure id_saveid (id, line)
+
+pointer id # IDENTIFY structure
+int line[2] # Save as line
+
+int i
+pointer id1
+
+begin
+ # Check if already saved. If not saved allocate memory.
+ for (i = 1; i <= ID_NID(id); i = i + 1) {
+ id1 = Memi[ID_ID(id)+i-1]
+ if (ID_LINE(id1,1) == line[1] && ID_LINE(id1,2) == line[2])
+ break
+ }
+
+ call id_sid (id, i)
+end
+
+
+
+# ID_SID -- Save identify information by index.
+
+procedure id_sid (id, n)
+
+pointer id # IDENTIFY structure
+int n # Save as index
+
+int i, j, dcvstati(), strlen()
+pointer sp, id1, coeffs, ptr1, ptr2
+
+begin
+ if (n > ID_NID(id)) {
+ if (n == 1)
+ call malloc (ID_ID(id), 1, TY_POINTER)
+ else
+ call realloc (ID_ID(id), n, TY_POINTER)
+ call calloc (id1, LEN_IDSTRUCT, TY_STRUCT)
+ Memi[ID_ID(id)+n-1] = id1
+ ID_NID(id) = n
+ } else
+ id1 = Memi[ID_ID(id)+n-1]
+
+ # Allocate or reallocate memory for features and copy them.
+ if (ID_NFEATURES(id) > 0) {
+ if (ID_NALLOC(id1) == 0) {
+ call malloc (ID_PIX(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call malloc (ID_FIT(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call malloc (ID_USER(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call malloc (ID_WTS(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call malloc (ID_FWIDTHS(id1), ID_NFEATURES(id), TY_REAL)
+ call malloc (ID_FTYPES(id1), ID_NFEATURES(id), TY_INT)
+ call calloc (ID_LABEL(id1), ID_NFEATURES(id), TY_POINTER)
+ } else if (ID_NALLOC(id1) != ID_NFEATURES(id)) {
+ call realloc (ID_PIX(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call realloc (ID_FIT(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call realloc (ID_USER(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call realloc (ID_WTS(id1), ID_NFEATURES(id), TY_DOUBLE)
+ call realloc (ID_FWIDTHS(id1), ID_NFEATURES(id), TY_REAL)
+ call realloc (ID_FTYPES(id1), ID_NFEATURES(id), TY_INT)
+ call realloc (ID_LABEL(id1), ID_NFEATURES(id), TY_POINTER)
+
+ j = ID_NALLOC(id1)
+ i = ID_NFEATURES(id) - j
+ if (i > 0)
+ call aclri (Memi[ID_LABEL(id1)+j], i)
+ }
+ call amovd (PIX(id,1), PIX(id1,1), ID_NFEATURES(id))
+ call amovd (FIT(id,1), FIT(id1,1), ID_NFEATURES(id))
+ call amovd (USER(id,1), USER(id1,1), ID_NFEATURES(id))
+ call amovd (WTS(id,1), WTS(id1,1), ID_NFEATURES(id))
+ call amovr (FWIDTH(id,1), FWIDTH(id1,1), ID_NFEATURES(id))
+ call amovi (FTYPE(id,1), FTYPE(id1,1), ID_NFEATURES(id))
+
+ ptr1 = ID_LABEL(id)
+ ptr2 = ID_LABEL(id1)
+ do i = 1, ID_NFEATURES(id) {
+ call mfree (Memi[ptr2], TY_CHAR)
+ if (Memi[ptr1] != NULL) {
+ j = strlen (Memc[Memi[ptr1]])
+ call malloc (Memi[ptr2], j, TY_CHAR)
+ call strcpy (Memc[Memi[ptr1]], Memc[Memi[ptr2]], j)
+ }
+ ptr1 = ptr1 + 1
+ ptr2 = ptr2 + 1
+ }
+
+ ID_NALLOC(id1) = ID_NFEATURES(id)
+ }
+
+ # Use a SAVE and RESTORE to copy the CURFIT data.
+ if (ID_CV(id1) != NULL)
+ call dcvfree (ID_CV(id1))
+ if (ID_CV(id) != NULL) {
+ call smark (sp)
+ i = dcvstati (ID_CV(id), CVNSAVE)
+ call salloc (coeffs, i, TY_DOUBLE)
+ call dcvsave (ID_CV(id), Memd[coeffs])
+ call dcvrestore (ID_CV(id1), Memd[coeffs])
+ call sfree (sp)
+
+ if (ID_IC(id1) == NULL)
+ call ic_open (ID_IC(id1))
+ call ic_copy (ID_IC(id), ID_IC(id1))
+ }
+
+ #ID_LINE(id1,1) = line[1]
+ #ID_LINE(id1,2) = line[2]
+ ID_LINE(id1,1) = ID_LINE(id,1)
+ ID_LINE(id1,2) = ID_LINE(id,2)
+ ID_AP(id1,1) = ID_AP(id,1)
+ ID_AP(id1,2) = ID_AP(id,2)
+ ID_NFEATURES(id1) = ID_NFEATURES(id)
+ ID_SHIFT(id1) = ID_SHIFT(id)
+ ID_REDSHIFT(id1) = ID_REDSHIFT(id)
+ ID_RMSRED(id1) = ID_RMSRED(id)
+ ID_ZHELIO(id1) = ID_ZHELIO(id)
+ ID_CURRENT(id1) = ID_CURRENT(id)
+
+ ID_NEWFEATURES(id1) = ID_NEWFEATURES(id)
+ ID_NEWCV(id1) = ID_NEWCV(id)
+ ID_NEWDBENTRY(id1) = ID_NEWDBENTRY(id)
+end
+
+
+# ID_GID -- Get saved identify information for specified line.
+
+int procedure id_gid (id, line)
+
+pointer id # IDENTIFY structure
+int line[2] # Line number to get
+
+int i, id_getid()
+
+begin
+ # Check if saved.
+ for (i = 1; i <= ID_NID(id); i = i + 1) {
+ if (ID_LINE(Memi[ID_ID(id)+i-1],1) == line[1] &&
+ ID_LINE(Memi[ID_ID(id)+i-1],2) == line[2])
+ break
+ }
+
+ return (id_getid (id, i))
+end
+
+
+# ID_GETID -- Get saved identify information for specified index.
+
+int procedure id_getid (id, n)
+
+pointer id # IDENTIFY structure
+int n # Index of saved features to be returned
+
+int i, j, dcvstati(), strlen()
+pointer sp, id1, coeffs, ptr1, ptr2
+
+begin
+ if (n < 1 || n > ID_NID(id))
+ return (EOF)
+
+ id1 = Memi[ID_ID(id)+n-1]
+
+ # Reallocate memory for features and copy them.
+ if (ID_NFEATURES(id1) > 0) {
+ if (ID_NALLOC(id1) != ID_NALLOC(id)) {
+ call realloc (ID_PIX(id), ID_NALLOC(id1), TY_DOUBLE)
+ call realloc (ID_FIT(id), ID_NALLOC(id1), TY_DOUBLE)
+ call realloc (ID_USER(id), ID_NALLOC(id1), TY_DOUBLE)
+ call realloc (ID_WTS(id), ID_NALLOC(id1), TY_DOUBLE)
+ call realloc (ID_FWIDTHS(id), ID_NALLOC(id1), TY_REAL)
+ call realloc (ID_FTYPES(id), ID_NALLOC(id1), TY_INT)
+ call realloc (ID_LABEL(id), ID_NALLOC(id1), TY_POINTER)
+
+ j = ID_NALLOC(id)
+ i = ID_NALLOC(id1) - j
+ if (i > 0)
+ call aclri (Memi[ID_LABEL(id)+j], i)
+ }
+ call amovd (PIX(id1,1), PIX(id,1), ID_NFEATURES(id1))
+ call amovd (FIT(id1,1), FIT(id,1), ID_NFEATURES(id1))
+ call amovd (USER(id1,1), USER(id,1), ID_NFEATURES(id1))
+ call amovd (WTS(id1,1), WTS(id,1), ID_NFEATURES(id1))
+ call amovr (FWIDTH(id1,1), FWIDTH(id,1), ID_NFEATURES(id1))
+ call amovi (FTYPE(id1,1), FTYPE(id,1), ID_NFEATURES(id1))
+
+ ptr1 = ID_LABEL(id1)
+ ptr2 = ID_LABEL(id)
+ do i = 1, ID_NFEATURES(id1) {
+ call mfree (Memi[ptr2], TY_CHAR)
+ if (Memi[ptr1] != NULL) {
+ j = strlen (Memc[Memi[ptr1]])
+ call malloc (Memi[ptr2], j, TY_CHAR)
+ call strcpy (Memc[Memi[ptr1]], Memc[Memi[ptr2]], j)
+ }
+ ptr1 = ptr1 + 1
+ ptr2 = ptr2 + 1
+ }
+
+ ID_NALLOC(id) = ID_NALLOC(id1)
+ ID_NFEATURES(id) = ID_NFEATURES(id1)
+ ID_NEWFEATURES(id) = ID_NEWFEATURES(id1)
+ ID_CURRENT(id) = ID_CURRENT(id1)
+ ID_NEWDBENTRY(id) = ID_NEWDBENTRY(id1)
+ }
+
+ # Use a SAVE and RESTORE to copy the CURFIT data.
+ if (ID_CV(id1) != NULL) {
+ if (ID_CV(id) != NULL)
+ call dcvfree (ID_CV(id))
+ call smark (sp)
+ i = dcvstati (ID_CV(id1), CVNSAVE)
+ call salloc (coeffs, i, TY_DOUBLE)
+ call dcvsave (ID_CV(id1), Memd[coeffs])
+ call dcvrestore (ID_CV(id), Memd[coeffs])
+ call sfree (sp)
+
+ call ic_copy (ID_IC(id1), ID_IC(id))
+
+ ID_SHIFT(id) = ID_SHIFT(id1)
+ ID_REDSHIFT(id) = ID_REDSHIFT(id1)
+ ID_RMSRED(id) = ID_RMSRED(id1)
+ ID_ZHELIO(id) = ID_ZHELIO(id1)
+ ID_NEWCV(id) = ID_NEWCV(id1)
+ ID_NEWDBENTRY(id) = ID_NEWDBENTRY(id1)
+ }
+
+ ID_LINE(id,1) = ID_LINE(id1,1)
+ ID_LINE(id,2) = ID_LINE(id1,2)
+ ID_AP(id,1) = ID_AP(id1,1)
+ ID_AP(id,2) = ID_AP(id1,2)
+
+ return (OK)
+end
diff --git a/noao/rv/rvidlines/idlabel.x b/noao/rv/rvidlines/idlabel.x
new file mode 100644
index 00000000..cb5fa439
--- /dev/null
+++ b/noao/rv/rvidlines/idlabel.x
@@ -0,0 +1,30 @@
+define SKIP ($1==' '||$1=='\t'||$1=='"'||$1=='\'')
+
+# ID_LABEL -- Set label
+
+procedure id_label (str, label)
+
+char str[ARB] # String to be set
+pointer label # Label pointer to be set
+
+int i, j, strlen()
+pointer cp
+
+begin
+ call mfree (label, TY_CHAR)
+
+ for (i=1; str[i]!=EOS && SKIP(str[i]); i=i+1)
+ ;
+ for (j=strlen(str); j>=i && SKIP(str[j]); j=j-1)
+ ;
+
+ if (i <= j) {
+ call malloc (label, j-i+1, TY_CHAR)
+ cp = label
+ for (; i<=j; i=i+1) {
+ Memc[cp] = str[i]
+ cp = cp + 1
+ }
+ Memc[cp] = EOS
+ }
+end
diff --git a/noao/rv/rvidlines/idlinelist.x b/noao/rv/rvidlines/idlinelist.x
new file mode 100644
index 00000000..d427aad7
--- /dev/null
+++ b/noao/rv/rvidlines/idlinelist.x
@@ -0,0 +1,250 @@
+include <error.h>
+include <mach.h>
+include "identify.h"
+
+# ID_MAPLL -- Read the line list into memory.
+
+procedure id_mapll (id)
+
+pointer id # Identify structure
+
+int fd, nalloc, nlines
+pointer ll1, ll2, str
+
+int open(), fscan(), nscan(), nowhite()
+double value
+
+errchk open, fscan, malloc, realloc
+
+begin
+ ID_LL(id) = NULL
+
+ if (nowhite (Memc[ID_COORDLIST(id)],
+ Memc[ID_COORDLIST(id)], SZ_FNAME) == 0)
+ return
+ iferr (fd = open (Memc[ID_COORDLIST(id)], READ_ONLY, TEXT_FILE)) {
+ call erract (EA_WARN)
+ return
+ }
+
+ call malloc (str, SZ_LINE, TY_CHAR)
+ nalloc = 0
+ nlines = 0
+ while (fscan (fd) != EOF) {
+ call gargd (value)
+ if (nscan() != 1)
+ next
+
+ if (nalloc == 0) {
+ nalloc = 100
+ call malloc (ll1, nalloc, TY_DOUBLE)
+ call calloc (ll2, nalloc, TY_POINTER)
+ } else if (nlines == nalloc) {
+ nalloc = nalloc + 100
+ call realloc (ll1, nalloc, TY_DOUBLE)
+ call realloc (ll2, nalloc, TY_POINTER)
+ call aclri (Memi[ll2+nalloc-100], 100)
+ }
+
+ Memd[ll1+nlines] = value
+ call gargstr (Memc[str], SZ_LINE)
+ call id_label (Memc[str], Memi[ll2+nlines])
+
+ nlines = nlines + 1
+ }
+ call mfree (str, TY_CHAR)
+ call close (fd)
+
+ if (nlines == 0)
+ return
+
+ call realloc (ll1, nlines + 1, TY_DOUBLE)
+ call realloc (ll2, nlines + 1, TY_POINTER)
+ Memd[ll1+nlines] = INDEFD
+ call malloc (ID_LL(id), 2, TY_POINTER)
+ Memi[ID_LL(id)] = ll1
+ Memi[ID_LL(id)+1] = ll2
+end
+
+
+# ID_UNMAPLL -- Unmap the linelist.
+
+procedure id_unmapll (id)
+
+pointer id # Identify structure
+
+pointer ll1, ll2
+
+begin
+ if (ID_LL(id) == NULL)
+ return
+
+ ll1 = Memi[ID_LL(id)]
+ ll2 = Memi[ID_LL(id)+1]
+ while (!IS_INDEFD(Memd[ll1])) {
+ call mfree (Memi[ll2], TY_CHAR)
+ ll1 = ll1 + 1
+ ll2 = ll2 + 1
+ }
+
+ call mfree (Memi[ID_LL(id)], TY_DOUBLE)
+ call mfree (Memi[ID_LL(id)+1], TY_POINTER)
+ call mfree (ID_LL(id), TY_POINTER)
+end
+
+
+
+# ID_MATCH -- Match current feature against a line list.
+#
+# This is extremely inefficient. It can be greatly improved.
+
+procedure id_match (id, in, out, label, diff)
+
+pointer id # Identify structure
+double in # Coordinate to be matched
+double out # Matched coordinate
+pointer label # Pointer to label
+real diff # Maximum difference
+
+double zin, delta, deltamin, id_zshiftd()
+pointer ll1, ll2, tmp
+int strlen()
+
+begin
+ if (ID_LL(id) == NULL) {
+ out = id_zshiftd (id, in, 0)
+ label = NULL
+ return
+ }
+
+ zin = id_zshiftd (id, in, 0)
+ deltamin = MAX_REAL
+
+ ll1 = Memi[ID_LL(id)]
+ ll2 = Memi[ID_LL(id)+1]
+ while (!IS_INDEFD(Memd[ll1])) {
+ delta = abs (zin - Memd[ll1])
+ if (delta < deltamin) {
+ deltamin = delta
+ if (deltamin <= diff) {
+ out = Memd[ll1]
+ label = Memi[ll2]
+ }
+ }
+ ll1 = ll1 + 1
+ ll2 = ll2 + 1
+ }
+
+ if (label != NULL) {
+ call malloc (tmp, strlen (Memc[label]), TY_CHAR)
+ call strcpy (Memc[label], Memc[tmp], ARB)
+ label = tmp
+ }
+end
+
+# ID_LINELIST -- Add features from a line list.
+
+procedure id_linelist (id)
+
+pointer id # Identify structure
+
+int i, nfound, nextpix, lastpix, cursave
+double pix, fit, fit1, fit2, user, peak, minval, diff, diff1
+pointer sp, pixes, fits, users, labels, ll1, ll2, label
+
+double id_center(), fit_to_pix(), id_fitpt(), id_peak(), id_zshiftd()
+
+begin
+ if (ID_LL(id) == NULL)
+ return
+
+ call smark (sp)
+ call salloc (pixes, ID_MAXFEATURES(id), TY_DOUBLE)
+ call salloc (fits, ID_MAXFEATURES(id), TY_DOUBLE)
+ call salloc (users, ID_MAXFEATURES(id), TY_DOUBLE)
+ call salloc (labels, ID_MAXFEATURES(id), TY_POINTER)
+
+ nfound = 0
+ lastpix = 0
+ minval = MAX_REAL
+
+ fit1 = min (FITDATA(id,1), FITDATA(id,ID_NPTS(id)))
+ fit2 = max (FITDATA(id,1), FITDATA(id,ID_NPTS(id)))
+ ll1 = Memi[ID_LL(id)]
+ ll2 = Memi[ID_LL(id)+1]
+ while (!IS_INDEFD(Memd[ll1])) {
+ user = id_zshiftd (id, Memd[ll1], 1)
+ label = Memi[ll2]
+ ll1 = ll1 + 1
+ ll2 = ll2 + 1
+ if (user < fit1)
+ next
+ if (user > fit2)
+ break
+
+ pix = id_center (id, fit_to_pix (id, user), 1, ID_FWIDTH(id),
+ ID_FTYPE(id), NO)
+ if (!IS_INDEFD(pix)) {
+ fit = id_fitpt (id, pix)
+ diff = abs (fit - user)
+ if (diff > ID_MATCH(id))
+ next
+ if (lastpix > 0) {
+ if (abs (pix - Memd[pixes+lastpix-1]) < 0.01) {
+ diff1 = abs (Memd[fits+lastpix-1]-Memd[users+lastpix-1])
+ if (diff < diff1) {
+ Memd[pixes+lastpix-1] = pix
+ Memd[fits+lastpix-1] = fit
+ Memd[users+lastpix-1] = id_zshiftd (id, user, 0)
+ Memi[labels+lastpix-1] = label
+ }
+ next
+ }
+ }
+
+ peak = abs (id_peak (id, pix))
+ if (nfound < ID_MAXFEATURES(id)) {
+ nfound = nfound + 1
+ if (peak < minval) {
+ nextpix = nfound
+ minval = peak
+ }
+ Memd[pixes+nfound-1] = pix
+ Memd[fits+nfound-1] = id_fitpt (id, pix)
+ Memd[users+nfound-1] = id_zshiftd (id, user, 0)
+ Memi[labels+nfound-1] = label
+ lastpix = nfound
+ } else if (peak > minval) {
+ Memd[pixes+nextpix-1] = pix
+ Memd[fits+nextpix-1] = id_fitpt (id, pix)
+ Memd[users+nextpix-1] = id_zshiftd (id, user, 0)
+ Memi[labels+nextpix-1] = label
+ lastpix = nextpix
+
+ minval = MAX_REAL
+ do i = 1, nfound {
+ pix = Memd[pixes+i-1]
+ peak = abs (id_peak (id, pix))
+ if (peak < minval) {
+ nextpix = i
+ minval = peak
+ }
+ }
+ }
+ }
+ }
+
+ do i = 1, nfound {
+ pix = Memd[pixes+i-1]
+ fit = Memd[fits+i-1]
+ user = Memd[users+i-1]
+ label = Memi[labels+i-1]
+ call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id),
+ ID_FTYPE(id), label)
+ if (i == 1)
+ cursave = ID_CURRENT(id)
+ }
+ ID_CURRENT(id) = cursave
+
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/idlog.x b/noao/rv/rvidlines/idlog.x
new file mode 100644
index 00000000..58645fc1
--- /dev/null
+++ b/noao/rv/rvidlines/idlog.x
@@ -0,0 +1,190 @@
+include <smw.h>
+include <time.h>
+include "identify.h"
+
+
+# ID_LOG -- Write log
+
+procedure id_log (id, file, logfd)
+
+pointer id # ID pointer
+char file[ARB] # Log file
+int logfd # Log fd
+
+char str[SZ_TIME]
+int i, fd, nrms
+double z, zrms, zerr, zhelio, resid, rms, v, verr, vhelio, hjd
+
+int open()
+double id_zshiftd(), id_zval()
+long clktime()
+errchk open, id_velocity, id_vhelio
+
+begin
+ if (ID_NFEATURES(id) == 0)
+ return
+
+ if (logfd == NULL)
+ fd = open (file, APPEND, TEXT_FILE)
+ else
+ fd = logfd
+
+ if (ID_TASK(id) == IDENTIFY) {
+
+ call cnvtime (clktime (0), str, SZ_TIME)
+ call fprintf (fd, "\n%s\n")
+ call pargstr (str)
+ call fprintf (fd, "Features identified in image %s%s: %s\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargstr (TITLE(ID_SH(id)))
+
+ call fprintf (fd, " %8s %10s %10s %10s %5s %s\n")
+ call pargstr ("Pixel")
+ call pargstr ("Fit")
+ call pargstr ("User")
+ call pargstr ("Residual")
+ call pargstr ("Fwidth")
+ call pargstr ("Wt")
+ call pargstr ("Label")
+
+ rms = 0.
+ nrms = 0
+ do i = 1, ID_NFEATURES(id) {
+ call fprintf (fd,
+ "%2d %8.2f %10.8g %10.8g %10.4g %4f %s\n")
+ call pargi (i)
+ call pargd (PIX(id,i))
+ call pargd (FIT(id,i))
+ call pargd (USER(id,i))
+ if (IS_INDEFD (USER(id,i)))
+ call pargd (USER(id,i))
+ else {
+ resid = FIT(id,i) - USER(id,i)
+ call pargd (resid)
+ if (WTS(id,i) > 0.) {
+ rms = rms + resid ** 2
+ nrms = nrms + 1
+ }
+ }
+ call pargd (WTS(id,i))
+ if (Memi[ID_LABEL(id)+i-1] != NULL)
+ call pargstr (Memc[Memi[ID_LABEL(id)+i-1]])
+ else
+ call pargstr ("")
+ }
+
+ if (nrms > 1) {
+ call fprintf (fd, "RMS = %0.6g\n")
+ call pargd (sqrt (rms / nrms))
+ }
+
+ } else {
+ call id_velocity (id, NO)
+ z = ID_REDSHIFT(id)
+ zrms = ID_RMSRED(id)
+ zhelio = ID_ZHELIO(id)
+ v = ID_REDSHIFT(id) * VLIGHT
+ call id_vhelio (IM(ID_SH(id)), vhelio, hjd, fd)
+
+ call cnvtime (clktime (0), str, SZ_TIME)
+ call fprintf (fd, "\n%s\n")
+ call pargstr (str)
+ call fprintf (fd, "Features identified in image %s%s: %s\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargstr (TITLE(ID_SH(id)))
+
+ call fprintf (fd, "%10s %10s %10s %10s %10s %5s %s\n")
+ call pargstr ("Measured")
+ call pargstr ("User")
+ call pargstr ("Residual")
+ call pargstr ("Velocity")
+ call pargstr ("Residual")
+ call pargstr ("Wt")
+ call pargstr ("Label")
+
+ rms = 0.
+ nrms = 0
+ do i = 1, ID_NFEATURES(id) {
+ call fprintf (fd,
+ "%10.8g %10.8g %10.4g %10.8g %10.4g %4f %s\n")
+ call pargd (id_zshiftd (id, FIT(id,i), 0))
+ if (IS_INDEFD (USER(id,i))) {
+ call pargd (INDEFD)
+ call pargd (INDEFD)
+ call pargd (INDEFD)
+ call pargd (INDEFD)
+ } else {
+ call pargd (USER(id,i))
+ resid = id_zshiftd (id, FIT(id,i), 0) - USER(id,i)
+ call pargd (resid)
+ if (WTS(id,i) > 0.) {
+ rms = rms + resid ** 2
+ nrms = nrms + 1
+ }
+ verr = id_zval (id, FIT(id,i), USER(id,i)) * VLIGHT
+ call pargd (verr + vhelio)
+ call pargd (verr - v)
+ }
+ call pargd (WTS(id,i))
+ if (Memi[ID_LABEL(id)+i-1] != NULL)
+ call pargstr (Memc[Memi[ID_LABEL(id)+i-1]])
+ else
+ call pargstr ("")
+ }
+
+ if (nrms > 1) {
+ call fprintf (fd, "Wavelength RMS = %0.6g\n")
+ call pargd (sqrt (rms / nrms))
+ call fprintf (fd, "Velocity RMS = %8.5g\n")
+ call pargd (zrms * VLIGHT)
+ }
+
+ zerr = zrms
+ if (nrms > 1)
+ zerr = zerr / sqrt (nrms - 1.)
+ v = z * VLIGHT
+ verr = zerr * VLIGHT
+
+ call fprintf (fd, "\n")
+ call fprintf (fd, "%s %3d : Zobs = %10.5g, ")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_AP(id,1))
+ call pargd (z)
+ call fprintf (fd, "Mean err = %10.5g, Lines = %3d\n")
+ call pargd (zerr)
+ call pargi (nrms)
+ call fprintf (fd, "%s %3d : Vobs = %8.5g km/s, ")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_AP(id,1))
+ call pargd (v)
+ call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n")
+ call pargd (verr)
+ call pargi (nrms)
+ if (zhelio != 0D0) {
+ call fprintf (fd, "%s %3d : Zhelio = %10.5g, ")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_AP(id,1))
+ call pargd (z + zhelio)
+ call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n")
+ call pargd (zerr)
+ call pargi (nrms)
+ call fprintf (fd, "%s %3d : Vhelio = %8.5g km/s, ")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_AP(id,1))
+ call pargd (v + vhelio)
+ call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n")
+ call pargd (verr)
+ call pargi (nrms)
+ call fprintf (fd, "%s %3d : HJD = %g\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_AP(id,1))
+ call pargd (hjd)
+ }
+ call fprintf (fd, "\n")
+ }
+
+ if (logfd == NULL)
+ call close (fd)
+end
diff --git a/noao/rv/rvidlines/idmap.x b/noao/rv/rvidlines/idmap.x
new file mode 100644
index 00000000..5af57217
--- /dev/null
+++ b/noao/rv/rvidlines/idmap.x
@@ -0,0 +1,379 @@
+include <ctype.h>
+include <imhdr.h>
+include <pkg/gtools.h>
+include <smw.h>
+include "identify.h"
+
+# Sepcial section words.
+define SPECIAL "|first|middle|x|y|z|last|column|line|band|"
+define FIRST 1
+define MIDDLE 2
+define X 3
+define Y 4
+define Z 5
+define LAST 6
+define COLUMN 7
+define LINE 8
+define BAND 9
+
+# ID_MAP -- Map an image for IDENTIFY/REIDENTIFY
+# The image must 1, 2, or 3 dimensional. An image section may be given with
+# the image name or with the CL parameter "section". The CL parameter can
+# have one of the following formats:
+# 1. An IMIO image section
+# 2. [line|column|x|y|z] [#|middle|last] [#|middle|last]
+# 3. [#|middle|last] [#|middle|last] [line|column|x|y|z]
+# where # is a line or column number. The strings may be abbreviated.
+# The task returns and error if it cannot map the image or determine
+# the 1D line or column desired.
+
+procedure id_map (id)
+
+pointer id # IDENTIFY data structure pointer
+
+int i, j, k, l, a, b, c, x1[3], x2[3], xs[3]
+pointer sp, wrd1, wrd2, wrd3, im
+
+int imaccess(), strdic(), ctoi(), nscan()
+pointer immap()
+errchk immap, id_maphdr
+
+begin
+ # Separate the image name and image section and map the full image.
+ call imgsection (Memc[ID_IMAGE(id)], Memc[ID_SECTION(id)], SZ_FNAME)
+ call imgimage (Memc[ID_IMAGE(id)], Memc[ID_IMAGE(id)], SZ_FNAME)
+ call id_noextn (Memc[ID_IMAGE(id)])
+ im = immap (Memc[ID_IMAGE(id)], READ_ONLY, 0)
+
+ # If no image section is found use the "section" parameter.
+ if (Memc[ID_SECTION(id)] == EOS && IM_NDIM(im) > 1) {
+ call clgstr ("section", Memc[ID_SECTION(id)], SZ_FNAME)
+ call xt_stripwhite (Memc[ID_SECTION(id)])
+
+ # If not an image section construct one.
+ if (Memc[ID_SECTION(id)] != '[') {
+ call smark (sp)
+ call salloc (wrd1, SZ_FNAME, TY_CHAR)
+ call salloc (wrd2, SZ_FNAME, TY_CHAR)
+ call salloc (wrd3, SZ_FNAME, TY_CHAR)
+
+ call sscan (Memc[ID_SECTION(id)])
+
+ # Parse axis and elements.
+ call gargwrd (Memc[wrd1], SZ_FNAME)
+ call gargwrd (Memc[wrd2], SZ_FNAME)
+ call gargwrd (Memc[wrd3], SZ_FNAME)
+ switch (nscan()) {
+ case 0:
+ a = X
+ b = MIDDLE
+ c = MIDDLE
+ case 1:
+ a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL)
+ b = MIDDLE
+ c = MIDDLE
+ case 2:
+ a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL)
+ if (a >= X)
+ b = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL)
+ else {
+ b = a
+ a = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL)
+ call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME)
+ }
+ c = MIDDLE
+ call strcpy (Memc[wrd2], Memc[wrd3], SZ_FNAME)
+ case 3:
+ a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL)
+ if (a >= X) {
+ b = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL)
+ c = strdic (Memc[wrd3], Memc[wrd3], SZ_FNAME, SPECIAL)
+ } else {
+ b = a
+ a = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL)
+ if (a >= X) {
+ c = strdic (Memc[wrd3], Memc[wrd3],SZ_FNAME,SPECIAL)
+ call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME)
+ } else {
+ c = b
+ b = a
+ a = strdic (Memc[wrd3], Memc[wrd3],SZ_FNAME,SPECIAL)
+ call strcpy (Memc[wrd2], Memc[wrd3], SZ_FNAME)
+ call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME)
+ }
+ }
+ }
+
+ switch (a) {
+ case X, LINE:
+ i = 1
+ j = 2
+ k = 3
+ case Y, COLUMN:
+ i = 2
+ j = 1
+ k = 3
+ case Z, BAND:
+ i = 3
+ j = 1
+ k = 2
+ default:
+ call imunmap (im)
+ call error (1,
+ "Error in section specification or non-unique abbreviation")
+ }
+
+ switch (b) {
+ case FIRST:
+ ID_LINE(id,1) = 1
+ case MIDDLE:
+ ID_LINE(id,1) = (1 + IM_LEN(im,j)) / 2
+ case LAST:
+ ID_LINE(id,1) = IM_LEN(im,j)
+ default:
+ l = 1
+ if (ctoi (Memc[wrd2], l, ID_LINE(id,1)) == 0)
+ call error (1, "Error in section specification")
+ }
+
+ switch (c) {
+ case FIRST:
+ ID_LINE(id,2) = 1
+ case MIDDLE:
+ ID_LINE(id,2) = (1 + IM_LEN(im,k)) / 2
+ case LAST:
+ ID_LINE(id,2) = IM_LEN(im,k)
+ default:
+ l = 1
+ if (ctoi (Memc[wrd3], l, ID_LINE(id,2)) == 0)
+ call error (1, "Error in section specification")
+ }
+
+ # Format section.
+ switch (IM_NDIM(im)) {
+ case 2:
+ switch (i) {
+ case 1:
+ call sprintf (Memc[ID_SECTION(id)], SZ_FNAME, "[*,%d]")
+ case 2:
+ call sprintf (Memc[ID_SECTION(id)], SZ_FNAME, "[%d,*]")
+ default:
+ call error (1, "Error in section specification")
+ }
+ call pargi (ID_LINE(id,1))
+ case 3:
+ switch (i) {
+ case 1:
+ call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[*,%d,%d]")
+ case 2:
+ call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[%d,*,%d]")
+ case 3:
+ call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[%d,%d,*]")
+ }
+ call pargi (ID_LINE(id,1))
+ call pargi (ID_LINE(id,2))
+ case 4:
+ call error (1, "Image dimension greater than 3 not allowed")
+ }
+
+ call sfree (sp)
+ }
+ }
+
+ # Parse the image section.
+ x1[1] = 1; x2[1] = IM_LEN(im,1); xs[1] = 1
+ x1[2] = 1; x2[2] = IM_LEN(im,2); xs[2] = 1
+ x1[3] = 1; x2[3] = IM_LEN(im,3); xs[3] = 1
+ call id_section (Memc[ID_SECTION(id)], x1, x2, xs, 3)
+
+ # Set the axes. The axis to be identified is the longest one.
+ i = 1
+ if (IM_NDIM(im) > 1 && abs (x1[2]-x2[2]) >= abs (x1[i]-x2[i]))
+ i = 2
+ if (IM_NDIM(im) > 2 && abs (x1[3]-x2[3]) >= abs (x1[i]-x2[i]))
+ i = 3
+ if (IM_NDIM(im) > 3)
+ call error (1, "Image dimension greater than 3 not allowed")
+
+ switch (i) {
+ case 1:
+ j = 2
+ k = 3
+ case 2:
+ j = 1
+ k = 3
+ case 3:
+ j = 1
+ k = 2
+ }
+
+ ID_LINE(id,1) = (x1[j] + x2[j]) / 2
+ ID_LINE(id,2) = (x1[k] + x2[k]) / 2
+ ID_MAXLINE(id,1) = IM_LEN(im, j)
+ ID_MAXLINE(id,2) = IM_LEN(im, k)
+ ID_NSUM(id,1) = min (ID_MAXLINE(id,1), ID_NSUM(id,1))
+ ID_NSUM(id,2) = min (ID_MAXLINE(id,2), ID_NSUM(id,2))
+ call smw_daxis (NULL, NULL, i, ID_NSUM(id,1), ID_NSUM(id,2))
+
+ call id_maphdr (id, im)
+
+ # Open the image READ_WRITE if possible in order to add REFSPEC.
+ # This is not done earlier to avoid updating of the WCS.
+
+ call imunmap (im)
+ if (imaccess (Memc[ID_IMAGE(id)], READ_WRITE) == YES)
+ im = immap (Memc[ID_IMAGE(id)], READ_WRITE, 0)
+ else
+ im = immap (Memc[ID_IMAGE(id)], READ_ONLY, 0)
+ IM(ID_SH(id)) = im
+end
+
+
+# ID_MAPHDR -- Map image header.
+
+procedure id_maphdr (id, im)
+
+pointer id # ID pointer
+pointer im # IMIO pointer
+
+int i
+pointer mw, sh, smw_openim(), smw_sctran()
+errchk smw_openim(), shdr_open(), smw_sctran
+
+begin
+ mw = smw_openim (im)
+ if (SMW_TRANS(mw) == YES) {
+ if (SMW_PAXIS(mw,1) == 1)
+ call smw_daxis (mw, im, 2, INDEFI, INDEFI)
+ else
+ call smw_daxis (mw, im, 1, INDEFI, INDEFI)
+ call smw_saxes (mw, NULL, im)
+ }
+ call shdr_open (im, mw, ID_LINE(id,1), ID_LINE(id,2),
+ INDEFI, SHHDR, ID_SH(id))
+ sh = ID_SH(id)
+
+ if (SMW_FORMAT(mw) == SMW_MS || SMW_FORMAT(mw) == SMW_ES) {
+ ID_MAXLINE(id,1) = IM_LEN(im,2)
+ ID_MAXLINE(id,2) = IM_LEN(im,3)
+ ID_NSUM(id,1) = 1
+ ID_NSUM(id,2) = 1
+ ID_LINE(id,1) = max (1, min (ID_MAXLINE(id,1), ID_LINE(id,1)))
+ ID_LINE(id,2) = 1
+ call mfree (ID_APS(id), TY_INT)
+ call malloc (ID_APS(id), ID_MAXLINE(id,1), TY_INT)
+ do i = 1, ID_MAXLINE(id,1) {
+ call shdr_open (im, mw, i, 1, INDEFI, SHHDR, sh)
+ Memi[ID_APS(id)+i-1] = AP(sh)
+ }
+ ID_AP(id,1) = Memi[ID_APS(id)+ID_LINE(id,1)-1]
+ ID_AP(id,2) = 1
+ } else {
+ call mfree (ID_APS(id), TY_INT)
+ ID_AP(id,1) = ID_LINE(id,1)
+ ID_AP(id,2) = ID_LINE(id,2)
+ }
+ ID_NPTS(id) = IM_LEN(im, SMW_LAXIS(mw,1))
+
+ call gt_sets (ID_GT(id), GTXLABEL, LABEL(sh))
+ call ic_pstr (ID_IC(id), "ylabel", LABEL(sh))
+ call gt_sets (ID_GT(id), GTXUNITS, UNITS(sh))
+ call ic_pstr (ID_IC(id), "yunits", UNITS(sh))
+
+ # Set logical / physical transformations
+ i = 2 ** (SMW_PAXIS(mw,1) - 1)
+ ID_LP(id) = smw_sctran (mw, "logical", "physical", i)
+ ID_PL(id) = smw_sctran (mw, "physical", "logical", i)
+end
+
+
+# ID_SECTION -- Parse an image section into its elements.
+# 1. The default values must be set by the caller.
+# 2. A null image section is OK.
+# 3. The first nonwhitespace character must be '['.
+# 4. The last interpreted character must be ']'.
+#
+# This procedure should be replaced with an IMIO procedure at some
+# point.
+
+procedure id_section (section, x1, x2, xs, ndim)
+
+char section[ARB] # Image section
+int x1[ndim] # Starting pixel
+int x2[ndim] # Ending pixel
+int xs[ndim] # Step
+int ndim # Number of dimensions
+
+int i, ip, a, b, c, temp, ctoi()
+define error_ 99
+
+begin
+ # Decode the section string.
+ ip = 1
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == '[')
+ ip = ip + 1
+ else if (section[ip] == EOS)
+ return
+ else
+ goto error_
+
+ do i = 1, ndim {
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == ']')
+ break
+
+ # Default values
+ a = x1[i]
+ b = x2[i]
+ c = xs[i]
+
+ # Get a:b:c. Allow notation such as "-*:c"
+ # (or even "-:c") where the step is obviously negative.
+
+ if (ctoi (section, ip, temp) > 0) { # a
+ a = temp
+ if (section[ip] == ':') {
+ ip = ip + 1
+ if (ctoi (section, ip, b) == 0) # a:b
+ goto error_
+ } else
+ b = a
+ } else if (section[ip] == '-') { # -*
+ temp = a
+ a = b
+ b = temp
+ ip = ip + 1
+ if (section[ip] == '*')
+ ip = ip + 1
+ } else if (section[ip] == '*') # *
+ ip = ip + 1
+ if (section[ip] == ':') { # ..:step
+ ip = ip + 1
+ if (ctoi (section, ip, c) == 0)
+ goto error_
+ else if (c == 0)
+ goto error_
+ }
+ if (a > b && c > 0)
+ c = -c
+
+ x1[i] = a
+ x2[i] = b
+ xs[i] = c
+
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == ',')
+ ip = ip + 1
+ }
+
+ if (section[ip] != ']')
+ goto error_
+
+ return
+error_
+ call error (0, "Error in image section specification")
+end
diff --git a/noao/rv/rvidlines/idmark.x b/noao/rv/rvidlines/idmark.x
new file mode 100644
index 00000000..8f874560
--- /dev/null
+++ b/noao/rv/rvidlines/idmark.x
@@ -0,0 +1,97 @@
+include <gset.h>
+include <smw.h>
+include "identify.h"
+
+procedure id_mark (id, feature)
+
+pointer id # ID pointer
+int feature
+
+int pix, color, gstati()
+real x, y
+real mx, my, x1, x2, y1, y2, tick, gap
+pointer sp, format, label, ptr
+double smw_c1trand()
+
+define TICK .03 # Tick size in NDC
+define GAP .02 # Gap size in NDC
+
+begin
+ call ggwind (ID_GP(id), x1, x2, y1, y2)
+
+ x = FIT(id,feature)
+
+ if ((x < min (x1, x2)) || (x > max (x1, x2)))
+ return
+
+ pix = smw_c1trand (ID_PL(id), PIX(id,feature)) - NP1(ID_SH(id)) + 1
+ pix = max (1, min (pix, ID_NPTS(id)-1))
+
+ call smark (sp)
+ call salloc (format, SZ_LINE, TY_CHAR)
+ call salloc (label, SZ_LINE, TY_CHAR)
+ switch (FTYPE(id,feature)) {
+ case EMISSION, GEMISSION:
+ y = max (IMDATA(id,pix), IMDATA(id,pix+1))
+ tick = TICK
+ gap = GAP
+ call strcpy ("u=180;h=c;v=b;s=0.5", Memc[format], SZ_LINE)
+ case ABSORPTION, GABSORPTION:
+ y = min (IMDATA(id,pix), IMDATA(id,pix+1))
+ tick = -TICK
+ gap = -GAP
+ call strcpy ("u=0;h=c;v=t;s=0.5", Memc[format], SZ_LINE)
+ }
+
+ call gctran (ID_GP(id), x, y, mx, my, 1, 0)
+ call gctran (ID_GP(id), mx, my + gap, x1, y1, 0, 1)
+ call gctran (ID_GP(id), mx, my + gap + tick, x1, y2, 0, 1)
+ color = gstati (ID_GP(id), G_PLCOLOR)
+ call gseti (ID_GP(id), G_PLCOLOR, color+1)
+ call gline (ID_GP(id), x1, y1, x1, y2)
+ call gseti (ID_GP(id), G_PLCOLOR, color)
+
+ call gctran (ID_GP(id), mx, my + tick + 2 * gap, x1, y2, 0, 1)
+ color = gstati (ID_GP(id), G_TXCOLOR)
+ call gseti (ID_GP(id), G_TXCOLOR, color+1)
+ switch (ID_LABELS(id)) {
+ case 2:
+ call sprintf (Memc[label], SZ_LINE, "%d")
+ call pargi (feature)
+ call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format])
+ case 3:
+ call sprintf (Memc[label], SZ_LINE, "%0.2f")
+ call pargd (PIX(id,feature))
+ call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format])
+ case 4:
+ if (!IS_INDEFD (USER(id,feature))) {
+ call sprintf (Memc[label], SZ_LINE, "%0.4f")
+ call pargd (USER(id,feature))
+ call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format])
+ }
+ case 5:
+ label = Memi[ID_LABEL(id)+feature-1]
+ if (label != NULL)
+ call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format])
+ case 6:
+ Memc[label] = EOS
+ ptr = Memi[ID_LABEL(id)+feature-1]
+ if (!IS_INDEFD (USER(id,feature))) {
+ if (ptr != NULL) {
+ call sprintf (Memc[label], SZ_LINE, "%0.4f %s")
+ call pargd (USER(id,feature))
+ call pargstr (Memc[ptr])
+ } else {
+ call sprintf (Memc[label], SZ_LINE, "%0.4f")
+ call pargd (USER(id,feature))
+ }
+ } else if (ptr != NULL)
+ call strcpy (Memc[ptr], Memc[label], SZ_LINE)
+ if (Memc[label] != EOS)
+ call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format])
+ }
+ call gseti (ID_GP(id), G_TXCOLOR, color)
+
+ call sfree (sp)
+ call gflush (ID_GP(id))
+end
diff --git a/noao/rv/rvidlines/idnearest.x b/noao/rv/rvidlines/idnearest.x
new file mode 100644
index 00000000..41aa4c61
--- /dev/null
+++ b/noao/rv/rvidlines/idnearest.x
@@ -0,0 +1,29 @@
+include "identify.h"
+
+# ID_NEAREST -- Find the nearest feature to a given coordinate.
+
+procedure id_nearest (id, fitnear)
+
+pointer id # ID pointer
+double fitnear # Coordinate to find nearest feature
+
+int i
+double delta, delta1
+
+begin
+ if (ID_NFEATURES(id) < 1) {
+ ID_CURRENT(id) = 0
+ return
+ }
+
+ ID_CURRENT(id) = 1
+ delta = abs (FIT(id,1) - fitnear)
+
+ do i = 2, ID_NFEATURES(id) {
+ delta1 = abs (FIT(id,i) - fitnear)
+ if (delta1 < delta) {
+ ID_CURRENT(id) = i
+ delta = delta1
+ }
+ }
+end
diff --git a/noao/rv/rvidlines/idnewfeature.x b/noao/rv/rvidlines/idnewfeature.x
new file mode 100644
index 00000000..efa489b4
--- /dev/null
+++ b/noao/rv/rvidlines/idnewfeature.x
@@ -0,0 +1,87 @@
+include <mach.h>
+include "identify.h"
+
+# ID_NEWFEATURE -- Allocate and initialize memory for a new feature.
+
+procedure id_newfeature (id, pix, fit, user, wt, width, type, label)
+
+pointer id # ID pointer
+double pix # Pixel coordinate
+double fit # Fit coordinate
+double user # User coordinate
+double wt # Feature weight
+real width # Feature width
+int type # Feature type
+pointer label # Pointer to feature label
+
+int i, current, strlen()
+double delta
+
+define NALLOC 20 # Length of additional allocations
+
+begin
+ if (IS_INDEFD (pix))
+ return
+
+ delta = MAX_REAL
+ do i = 1, ID_NFEATURES(id) {
+ if (abs (pix - PIX(id,i)) < delta) {
+ delta = abs (pix - PIX(id,i))
+ current = i
+ }
+ }
+
+ if (delta >= ID_MINSEP(id)) {
+ ID_NFEATURES(id) = ID_NFEATURES(id) + 1
+ if (ID_NALLOC(id) < ID_NFEATURES(id)) {
+ ID_NALLOC(id) = ID_NALLOC(id) + NALLOC
+ call realloc (ID_PIX(id), ID_NALLOC(id), TY_DOUBLE)
+ call realloc (ID_FIT(id), ID_NALLOC(id), TY_DOUBLE)
+ call realloc (ID_USER(id), ID_NALLOC(id), TY_DOUBLE)
+ call realloc (ID_WTS(id), ID_NALLOC(id), TY_DOUBLE)
+ call realloc (ID_FWIDTHS(id), ID_NALLOC(id), TY_REAL)
+ call realloc (ID_FTYPES(id), ID_NALLOC(id), TY_INT)
+ call realloc (ID_LABEL(id), ID_NALLOC(id), TY_POINTER)
+ call aclri (Memi[ID_LABEL(id)+ID_NALLOC(id)-NALLOC], NALLOC)
+ }
+ for (current=ID_NFEATURES(id); (current>1)&&(pix<PIX(id,current-1));
+ current=current-1) {
+ PIX(id,current) = PIX(id,current-1)
+ FIT(id,current) = FIT(id,current-1)
+ USER(id,current) = USER(id,current-1)
+ WTS(id,current) = WTS(id,current-1)
+ FWIDTH(id,current) = FWIDTH(id,current-1)
+ FTYPE(id,current) = FTYPE(id,current-1)
+ Memi[ID_LABEL(id)+current-1] = Memi[ID_LABEL(id)+current-2]
+ }
+ PIX(id,current) = pix
+ FIT(id,current) = fit
+ USER(id,current) = user
+ WTS(id,current) = wt
+ FWIDTH(id,current) = width
+ FTYPE(id,current) = type
+ if (label != NULL) {
+ i = strlen (Memc[label])
+ call malloc (Memi[ID_LABEL(id)+current-1], i, TY_CHAR)
+ call strcpy (Memc[label], Memc[Memi[ID_LABEL(id)+current-1]], i)
+ } else
+ Memi[ID_LABEL(id)+current-1] = NULL
+ ID_NEWFEATURES(id) = YES
+ } else if (abs (fit-user) < abs (FIT(id,current)-USER(id,current))) {
+ PIX(id,current) = pix
+ FIT(id,current) = fit
+ USER(id,current) = user
+ WTS(id,current) = wt
+ FWIDTH(id,current) = width
+ FTYPE(id,current) = type
+ if (label != NULL) {
+ i = strlen (Memc[label])
+ call malloc (Memi[ID_LABEL(id)+current-1], i, TY_CHAR)
+ call strcpy (Memc[label], Memc[Memi[ID_LABEL(id)+current-1]], i)
+ } else
+ Memi[ID_LABEL(id)+current-1] = NULL
+ ID_NEWFEATURES(id) = YES
+ }
+
+ ID_CURRENT(id) = current
+end
diff --git a/noao/rv/rvidlines/idnoextn.x b/noao/rv/rvidlines/idnoextn.x
new file mode 100644
index 00000000..6c82d778
--- /dev/null
+++ b/noao/rv/rvidlines/idnoextn.x
@@ -0,0 +1,11 @@
+# ID_NOEXTN -- Remove standard image extensions.
+
+procedure id_noextn (image)
+
+char image[ARB] # Image name
+
+int strlen()
+
+begin
+ call xt_imroot (image, image, strlen (image))
+end
diff --git a/noao/rv/rvidlines/idpeak.x b/noao/rv/rvidlines/idpeak.x
new file mode 100644
index 00000000..9cba49c4
--- /dev/null
+++ b/noao/rv/rvidlines/idpeak.x
@@ -0,0 +1,23 @@
+include "identify.h"
+
+# ID_PEAK -- Find the peak value above continuum.
+
+double procedure id_peak (id, pix)
+
+pointer id # ID pointer
+double pix # Pixel position
+double peak # Peak value
+
+int c, l, u
+
+begin
+ if (IS_INDEFD(pix))
+ return (INDEFD)
+
+ c = nint (pix)
+ l = max (1, nint (pix - ID_FWIDTH(id)))
+ u = min (ID_NPTS(id), nint (pix + ID_FWIDTH(id)))
+ peak = IMDATA(id,c) - (IMDATA(id,l) + IMDATA(id,u)) / 2.
+
+ return (peak)
+end
diff --git a/noao/rv/rvidlines/idrms.x b/noao/rv/rvidlines/idrms.x
new file mode 100644
index 00000000..257675de
--- /dev/null
+++ b/noao/rv/rvidlines/idrms.x
@@ -0,0 +1,28 @@
+include "identify.h"
+
+# ID_RMS -- Compute RMS of fit about the user coordinates
+
+double procedure id_rms (id)
+
+pointer id # ID pointer
+
+int i, nrms
+double rms, id_zshiftd()
+
+begin
+ rms = 0.
+ nrms = 0
+ for (i=1; i<=ID_NFEATURES(id); i=i+1) {
+ if (!IS_INDEFD (USER(id,i)) && WTS(id,i) != 0.) {
+ rms = rms + (id_zshiftd (id, FIT(id,i), 0) - USER(id,i)) ** 2
+ nrms = nrms + 1
+ }
+ }
+
+ if (nrms > 0)
+ rms = sqrt (rms / nrms)
+ else
+ rms = INDEFD
+
+ return (rms)
+end
diff --git a/noao/rv/rvidlines/idshift.x b/noao/rv/rvidlines/idshift.x
new file mode 100644
index 00000000..3fa9f9bf
--- /dev/null
+++ b/noao/rv/rvidlines/idshift.x
@@ -0,0 +1,65 @@
+include "identify.h"
+
+define NBIN 10 # Bin parameter for mode determination
+
+# ID_SHIFT -- Determine a shift by correlating feature user positions
+# with peaks in the image data.
+
+double procedure id_shift (id)
+
+pointer id # ID pointer
+
+int i, j, npeaks, ndiff, find_peaks()
+real d, dmin
+double pix, id_center(), id_fitpt()
+pointer x, y, diff
+errchk malloc, find_peaks
+
+begin
+ # Find the peaks in the image data and center.
+ call malloc (x, ID_NPTS(id), TY_REAL)
+ npeaks = find_peaks (IMDATA(id,1), Memr[x], ID_NPTS(id), 0.,
+ int (ID_MINSEP(id)), 0, ID_MAXFEATURES(id), 0., false)
+
+ # Center the peaks and convert to user coordinates.
+ call malloc (y, npeaks, TY_DOUBLE)
+ j = 0
+ do i = 1, npeaks {
+ pix = id_center (id, double(Memr[x+i-1]), 1, ID_FWIDTH(id),
+ ID_FTYPE(id), NO)
+ if (!IS_INDEFD (pix)) {
+ Memd[y+j] = id_fitpt (id, pix)
+ j = j + 1
+ }
+ }
+ npeaks = j
+
+ # Compute differences with feature list.
+ ndiff = npeaks * ID_NFEATURES(id)
+ call malloc (diff, ndiff, TY_REAL)
+ ndiff = 0
+ do i = 1, ID_NFEATURES(id) {
+ do j = 1, npeaks {
+ Memr[diff+ndiff] = Memd[y+j-1] - FIT(id,i)
+ ndiff = ndiff + 1
+ }
+ }
+ call mfree (x, TY_REAL)
+ call mfree (y, TY_DOUBLE)
+
+ # Sort the differences and find the mode.
+ call asrtr (Memr[diff], Memr[diff], ndiff)
+
+ dmin = Memr[diff+ndiff-1] - Memr[diff]
+ do i = 0, ndiff-NBIN-1 {
+ j = i + NBIN
+ d = Memr[diff+j] - Memr[diff+i]
+ if (d < dmin) {
+ dmin = d
+ pix = Memr[diff+i] + d / 2.
+ }
+ }
+ call mfree (diff, TY_REAL)
+
+ return (pix)
+end
diff --git a/noao/rv/rvidlines/idshow.x b/noao/rv/rvidlines/idshow.x
new file mode 100644
index 00000000..a01fba6f
--- /dev/null
+++ b/noao/rv/rvidlines/idshow.x
@@ -0,0 +1,83 @@
+include "identify.h"
+
+# ID_SHOW -- Show parameter information.
+
+procedure id_show (id, file)
+
+pointer id # ID pointer
+char file[ARB] # File
+
+char line[SZ_LINE]
+int fd
+
+int open(), ic_geti()
+errchk open()
+
+begin
+ fd = open (file, APPEND, TEXT_FILE)
+
+ call sysid (line, SZ_LINE)
+ call fprintf (fd, "%s\n")
+ call pargstr (line)
+
+ call fprintf (fd, "image %s\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call fprintf (fd, "nsum %d\n")
+ call pargi (ID_NSUM(id,1))
+ switch (ID_FTYPE(id)) {
+ case EMISSION:
+ call fprintf (fd, "ftype emission\n")
+ case ABSORPTION:
+ call fprintf (fd, "ftype absorption\n")
+ case GEMISSION:
+ call fprintf (fd, "ftype gemission\n")
+ case GABSORPTION:
+ call fprintf (fd, "ftype gabsorption\n")
+ }
+ switch (ID_LABELS(id)) {
+ case 2:
+ call fprintf (fd, "labels index\n")
+ case 3:
+ call fprintf (fd, "labels pixel\n")
+ case 4:
+ call fprintf (fd, "labels coords\n")
+ case 5:
+ call fprintf (fd, "labels user\n")
+ case 6:
+ call fprintf (fd, "labels both\n")
+ default:
+ call fprintf (fd, "labels none\n")
+ }
+ call fprintf (fd, "maxfeatures %d\n")
+ call pargi (ID_MAXFEATURES(id))
+ call fprintf (fd, "match %g\n")
+ call pargr (ID_MATCH(id))
+ call fprintf (fd, "zwidth %g\n")
+ call pargr (ID_ZWIDTH(id))
+ call fprintf (fd, "fwidth %g\n")
+ call pargr (ID_FWIDTH(id))
+ call fprintf (fd, "database %s\n")
+ call pargstr (Memc[ID_DATABASE(id)])
+ call fprintf (fd, "coordlist %s\n")
+ call pargstr (Memc[ID_COORDLIST(id)])
+ call fprintf (fd, "cradius %g\n")
+ call pargr (ID_CRADIUS(id))
+ call fprintf (fd, "threshold %g\n")
+ call pargr (ID_THRESHOLD(id))
+ call fprintf (fd, "minsep %g\n")
+ call pargr (ID_MINSEP(id))
+ if (ID_CV(id) != NULL) {
+ call fprintf (fd, "function = %s\n")
+ call ic_gstr (ID_IC(id), "function", line, SZ_LINE)
+ call pargstr (line)
+ call fprintf (fd, "order = %d\n")
+ call pargi (ic_geti (ID_IC(id), "order"))
+ call fprintf (fd, "Fit at first pixel = %0.8g\n")
+ call pargd (FITDATA(id,1))
+ call fprintf (fd, "Average fit interval = %0.8g\n")
+ call pargd ((FITDATA(id,ID_NPTS(id))-FITDATA(id,1))/
+ (ID_NPTS(id)-1))
+ }
+
+ call close (fd)
+end
diff --git a/noao/rv/rvidlines/idvelocity.x b/noao/rv/rvidlines/idvelocity.x
new file mode 100644
index 00000000..c62e5c89
--- /dev/null
+++ b/noao/rv/rvidlines/idvelocity.x
@@ -0,0 +1,188 @@
+include <smw.h>
+include <units.h>
+include "identify.h"
+
+
+# ID_VELOCITY -- Compute velocity.
+
+procedure id_velocity (id, interactive)
+
+pointer id # ID pointer
+int interactive # Called interactively?
+
+int i, n
+double z, sumz, sumz2, sumw, zerr, zhelio, v, verr, id_zval()
+
+begin
+ sumz = 0
+ sumw = 0
+ n = 0
+ for (i=1; i <= ID_NFEATURES(id); i = i + 1) {
+ if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.)
+ next
+ z = id_zval (id, FIT(id,i), USER(id,i))
+ sumz = sumz + WTS(id,i) * z
+ sumw = sumw + WTS(id,i)
+ n = n + 1
+ }
+
+ if (sumw > 0.) {
+ zhelio = ID_ZHELIO(id)
+ sumz = sumz / sumw
+
+ sumz2 = 0.
+ for (i=1; i <= ID_NFEATURES(id); i = i + 1) {
+ if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.)
+ next
+ z = id_zval (id, FIT(id,i), USER(id,i))
+ sumz2 = sumz2 + WTS(id,i) * (z - sumz) ** 2
+ }
+ if (sumz2 > 0.)
+ sumz2 = sqrt (sumz2 / sumw)
+ else
+ sumz2 = 0.
+ zerr = sumz2
+ if (n > 1)
+ zerr = zerr / sqrt (n - 1.)
+
+ if (interactive == YES) {
+ v = (sumz + zhelio) * VLIGHT
+ verr = zerr * VLIGHT
+
+ if (zhelio == 0D0)
+ call printf (
+ "%s%s: Lines=%3d, Vobs=%.5g (%.5g), Zobs=%.5g (%.5g)\n")
+ else
+ call printf (
+ "%s%s: Lines=%3d, Vhelio=%.5g (%.5g), Zhelio=%.5g (%.5g)\n")
+
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargi (n)
+ call pargd (v)
+ call pargd (verr)
+ call pargd (sumz + zhelio)
+ call pargd (zerr)
+ }
+ ID_REDSHIFT(id) = sumz
+ ID_RMSRED(id) = sumz2
+ ID_ZHELIO(id) = zhelio
+ }
+end
+
+
+# ID_ZVAL -- Compute Z value.
+
+double procedure id_zval (id, x, xref)
+
+pointer id #I Identify pointer
+double x #I Coordinate
+double xref #I Reference coordinate
+double z #O Z value
+
+double y, yref
+pointer un
+
+begin
+ y = x
+ yref = xref
+
+ un = UN(ID_SH(id))
+ if (UN_LOG(un) == YES) {
+ y = 10D0 ** y
+ yref = 10D0 ** yref
+ }
+ if (UN_INV(un) == YES) {
+ y = 1D0 / y
+ yref = 1D0 / yref
+ }
+
+ switch (UN_CLASS(un)) {
+ case UN_WAVE:
+ z = (y - yref) / yref
+ case UN_FREQ, UN_ENERGY:
+ z = (yref - y) / y
+ case UN_VEL:
+ y = sqrt ((1 + y) / (1 - y))
+ yref = sqrt ((1 + yref) / (1 - yref))
+ z = (y - yref) / yref
+ case UN_DOP:
+ y = y + 1
+ yref = yref + 1
+ z = (y - yref) / yref
+ }
+
+ return (z)
+end
+
+
+
+# ID_ZSHIFTD -- Shift coordinate by redshift.
+
+double procedure id_zshiftd (id, x, dir)
+
+pointer id #I Identify pointer
+double x #I Coordinate
+int dir #I Direction (0=to rest, 1=from rest)
+double y #O Shifted coordinate
+
+pointer un
+
+begin
+ y = x
+
+ un = UN(ID_SH(id))
+ if (UN_LOG(un) == YES)
+ y = 10D0 ** y
+ if (UN_INV(un) == YES)
+ y = 1D0 / y
+
+ switch (UN_CLASS(un)) {
+ case UN_WAVE:
+ if (dir == 0)
+ y = y / (1 + ID_REDSHIFT(id))
+ else
+ y = y * (1 + ID_REDSHIFT(id))
+ case UN_FREQ, UN_ENERGY:
+ if (dir == 0)
+ y = y * (1 + ID_REDSHIFT(id))
+ else
+ y = y / (1 + ID_REDSHIFT(id))
+ case UN_VEL:
+ y = sqrt ((1 + y) / (1 - y))
+ if (dir == 0)
+ y = y / (1 + ID_REDSHIFT(id))
+ else
+ y = y * (1 + ID_REDSHIFT(id))
+ y = y ** 2
+ y = (y - 1) / (y + 1)
+ case UN_DOP:
+ y = (y + 1)
+ if (dir == 0)
+ y = (y + 1) / (1 + ID_REDSHIFT(id)) - 1
+ else
+ y = (y + 1) * (1 + ID_REDSHIFT(id)) - 1
+ }
+
+ if (UN_INV(un) == YES)
+ y = 1D0 / y
+ if (UN_LOG(un) == YES)
+ y = log10 (y)
+
+ return (y)
+end
+
+
+# ID_ZSHIFTR -- Shift coordinate by redshift.
+
+real procedure id_zshiftr (id, x, dir)
+
+pointer id #I Identify pointer
+real x #I Coordinate
+int dir #I Direction (0=to rest, 1=from rest)
+
+double id_zshiftd()
+
+begin
+ return (real (id_zshiftd (id, double(x), dir)))
+end
diff --git a/noao/rv/rvidlines/idvhelio.x b/noao/rv/rvidlines/idvhelio.x
new file mode 100644
index 00000000..803efc2d
--- /dev/null
+++ b/noao/rv/rvidlines/idvhelio.x
@@ -0,0 +1,102 @@
+include <error.h>
+
+# ID_VHELIO -- Compute helocentric velocity.
+
+procedure id_vhelio (im, vhelio, hjd, fd)
+
+pointer im #I IMIO pointer
+double vhelio #O Heliocentric velocity correction
+double hjd #O Heliocentric Julian Date
+int fd #I Log file descriptor
+
+bool newobs, obshead
+int year, month, day, flags, dummy
+double ra, dec, ep, ut, lt
+double epoch, vrot, vbary, vorb
+double latitude, longitude, altitude
+pointer sp, str1, str2, tmp, obs, kp
+
+int dtm_decode()
+double imgetd(), obsgetd()
+pointer clopset()
+
+int err
+data err/0/
+
+errchk imgetd, imgstr, obsopen
+
+
+begin
+ call smark (sp)
+ call salloc (str1, SZ_LINE, TY_CHAR)
+ call salloc (str2, SZ_LINE, TY_CHAR)
+
+ obs = NULL
+
+ iferr {
+ # Get the observatory data.
+ call clgstr ("observatory", Memc[str1], SZ_LINE)
+ tmp = NULL
+ call obsimopen (tmp, im, Memc[str1], NO, newobs, obshead)
+ obs = tmp
+
+ latitude = obsgetd (obs, "latitude")
+ longitude = obsgetd (obs, "longitude")
+ altitude = obsgetd (obs, "altitude")
+
+ # Get the image header data.
+ kp = clopset ("keywpars")
+
+ call clgpset (kp, "date_obs", Memc[str1], SZ_LINE)
+ call imgstr (im, Memc[str1], Memc[str2], SZ_LINE)
+ if (dtm_decode (Memc[str2],year,month,day,ut,flags) == ERR)
+ call error (1, "Error in date string")
+
+ call clgpset (kp, "ut", Memc[str1], SZ_LINE)
+ call imgstr (im, Memc[str1], Memc[str2], SZ_LINE)
+ if (dtm_decode (Memc[str2],dummy,dummy,dummy,ut,flags) == ERR) {
+ iferr (ut = imgetd (im, Memc[str1]))
+ call error (1, "Error in UT keyword")
+ }
+
+ call clgpset (kp, "ra", Memc[str1], SZ_LINE)
+ ra = imgetd (im, Memc[str1])
+ call clgpset (kp, "dec", Memc[str1], SZ_LINE)
+ dec = imgetd (im, Memc[str1])
+ call clgpset (kp, "epoch", Memc[str1], SZ_LINE)
+ ep = imgetd (im, Memc[str1])
+ call clcpset (kp)
+
+ # Determine epoch of observation and precess coordinates.
+ call ast_date_to_epoch (year, month, day, ut, epoch)
+ call ast_precess (ra, dec, ep, ra, dec, epoch)
+
+ # Determine velocity components.
+ call ast_vorbit (ra, dec, epoch, vorb)
+ call ast_vbary (ra, dec, epoch, vbary)
+ call ast_vrotate (ra, dec, epoch, latitude, longitude,
+ altitude, vrot)
+ call ast_hjd (ra, dec, epoch, lt, hjd)
+
+ vhelio = vrot + vbary + vorb
+
+ if (fd != NULL)
+ call obslog (obs,
+ "RVIDLINES", "latitude longitude altitude", fd)
+
+ } then {
+ vhelio = 0.
+ if (err == 0) {
+ call eprintf ("Warning: Can't compute heliocentric velocity\n")
+ call erract (EA_WARN)
+ err = err + 1
+ }
+ }
+
+ iferr { # This IFERR is to clear errcode.
+ if (obs != NULL)
+ call obsclose (obs)
+ } then
+ ;
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/mkpkg b/noao/rv/rvidlines/mkpkg
new file mode 100644
index 00000000..a18105f8
--- /dev/null
+++ b/noao/rv/rvidlines/mkpkg
@@ -0,0 +1,47 @@
+# RVIDLINES
+
+$checkout libpkg.a ..
+$update libpkg.a
+$checkin libpkg.a ..
+$exit
+
+libpkg.a:
+ $ifeq (USE_GENERIC, yes)
+ $ifolder (peaks.x, peaks.gx)
+ $generic -k peaks.gx -o peaks.x $endif $endif
+
+ # Version from V2.10.3 ONEDSPEC that conflicts with V2.10.2 RV
+ iddeblend.x <mach.h>
+
+ idcenter.x identify.h <gset.h> <smw.h>
+ idcolon.x identify.h <error.h> <gset.h> <smw.h>
+ iddb.x identify.h <imset.h> <math/curfit.h> <smw.h>
+ iddelete.x identify.h
+ iddofit.x identify.h
+ iddoshift.x identify.h
+ idfitdata.x identify.h <smw.h>
+ idfixx.x <smw.h>
+ idgdata.x identify.h <imhdr.h> <imio.h> <pkg/gtools.h> <smw.h>
+ idgraph.x identify.h <gset.h> <pkg/gtools.h> <smw.h>
+ ididentify.x identify.h <error.h> <gset.h> <imhdr.h> <smw.h>
+ idinit.x identify.h <gset.h> <math/curfit.h>
+ idlabel.x
+ idlinelist.x identify.h <error.h> <mach.h>
+ idlog.x identify.h <smw.h> <time.h>
+ idmap.x identify.h <ctype.h> <imhdr.h> <pkg/gtools.h> <smw.h>
+ idmark.x identify.h <gset.h> <smw.h>
+ idnearest.x identify.h
+ idnewfeature.x identify.h <mach.h>
+ idnoextn.x
+ idpeak.x identify.h
+ idrms.x identify.h
+ idshift.x identify.h
+ idshow.x identify.h
+ idvelocity.x identify.h <smw.h> <units.h>
+ idvhelio.x <error.h>
+ peaks.x
+ reidentify.x identify.h <error.h> <gset.h> <imhdr.h>
+ t_identify.x identify.h <mach.h> <pkg/gtools.h>
+ t_reidentify.x identify.h <error.h> <fset.h> <gset.h> <pkg/gtools.h>\
+ <smw.h>
+ ;
diff --git a/noao/rv/rvidlines/peaks.gx b/noao/rv/rvidlines/peaks.gx
new file mode 100644
index 00000000..32764505
--- /dev/null
+++ b/noao/rv/rvidlines/peaks.gx
@@ -0,0 +1,447 @@
+# PEAKS -- The following procedures are general numerical functions
+# dealing with finding peaks in a data array.
+#
+# FIND_PEAKS Find the peaks in the data array.
+# FIND_LOCAL_MAXIMA Find the local minima/maxima in the data array.
+# IS_LOCAL_MAX Test a point to determine if it is a local min/max.
+# FIND_THRESHOLD Find the peaks with positions satisfying threshold
+# and contrast constraints.
+# FIND_ISOLATED Flag peaks which are within separation of a peak
+# with a higher peak value.
+# FIND_NMAX Select up to the nmax highest ranked peaks.
+# COMPARE Compare procedure for sort used in FIND_PEAKS.
+
+# FIND_PEAKS -- Find the peaks in the data array.
+#
+# The peaks are found using the following algorithm:
+#
+# 1. Find the local minima/maxima.
+# 2. Reject peaks below the threshold.
+# 3. Determine the ranks of the remaining peaks.
+# 4. Flag weaker peaks within separation of a stronger peak.
+# 5. Accept at most the nmax strongest peaks.
+#
+# Indefinite points are ignored. The peak positions are returned in the
+# array x.
+
+$for (r)
+int procedure find_peaks (data, x, npoints, contrast, separation, edge, nmax,
+ threshold, debug)
+
+# Procedure parameters:
+PIXEL data[npoints] # Input data array
+PIXEL x[npoints] # Output peak position array
+int npoints # Number of data points
+real contrast # Maximum contrast between strongest and weakest
+int separation # Minimum separation between peaks
+int edge # Minimum distance from the edge
+int nmax # Maximum number of peaks to be returned
+real threshold # Minimum threshold level for peaks
+bool debug # Print diagnostic information?
+
+int i, j
+int nlmax, nthreshold, nisolated, npeaks
+pointer sp, rank
+
+int find_local_maxima(), find_threshold(), find_isolated(), find_nmax()
+int compare()
+
+extern compare()
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # Find the local minima/maxima in data and put column positions in x..
+ if (contrast < 0.)
+ maxima = NO
+ else
+ maxima = YES
+
+ nlmax = find_local_maxima (data, x, npoints, debug)
+
+ # Reject local minima/maxima near the edge.
+ if (edge > 0) {
+ j = 0
+ do i = 1, nlmax {
+ if ((x[i] > edge) && (x[i] <= npoints - edge)) {
+ j = j + 1
+ x[j] = x[i]
+ }
+ }
+ nlmax = j
+ }
+
+ # Allocate a working array y.
+ call smark (sp)
+ call salloc (y, npoints, TY_PIXEL)
+
+ # Reject the local minima/maxima which do not satisfy the thresholds.
+ # The array y is set to the peak values of the remaining peaks.
+ nthreshold = find_threshold (data, x, Mem$t[y], nlmax,
+ contrast, threshold, debug)
+
+ # Rank the peaks by peak value.
+ call salloc (rank, nthreshold, TY_INT)
+ do i = 1, nthreshold
+ Memi[rank + i - 1] = i
+ call qsort (Memi[rank], nthreshold, compare)
+
+ # Reject the weaker peaks within sep of a stronger peak.
+ nisolated = find_isolated (x, Memi[rank], nthreshold, separation,
+ debug)
+
+ # Select the strongest nmax peaks.
+ npeaks = find_nmax (data, x, Memi[rank], nthreshold, nmax, debug)
+
+ call sfree (sp)
+ return (npeaks)
+end
+
+
+# FIND_LOCAL_MAXIMA -- Find the local minima/maxima in the data array.
+#
+# A data array is input and the local minima/maxima positions array is output.
+# The number of local minima/maxima found is returned.
+
+int procedure find_local_maxima (data, x, npoints, debug)
+
+PIXEL data[npoints] # Input data array
+PIXEL x[npoints] # Output local min/max positions array
+int npoints # Number of input points
+bool debug # Print debugging information?
+
+int i, nlmax
+
+bool is_local_max()
+
+begin
+ nlmax = 0
+ do i = 1, npoints {
+ if (is_local_max (i, data, npoints)) {
+ nlmax = nlmax + 1
+ x[nlmax] = i
+ }
+ }
+
+ if (debug) {
+ call printf (" Number of local minima/maxima found = %d.\n")
+ call pargi (nlmax)
+ }
+
+ return (nlmax)
+end
+
+
+# IS_LOCAL_MAX -- Test a point to determine if it is a local minima/maximum.
+#
+# Indefinite points are ignored.
+
+bool procedure is_local_max (index, data, npoints)
+
+# Procedure parameters:
+int index # Index to test for local minimum/maximum
+PIXEL data[npoints] # Data values
+int npoints # Number of points in the data vector
+
+int i, j, nright, nleft
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # INDEF points cannot be local minima/axima.
+ if (IS_INDEF (data[index]))
+ return (FALSE)
+
+ # Find the left and right indices where data values change and the
+ # number of points with the same value. Ignore INDEF points.
+ nleft = 0
+ for (i = index - 1; i >= 1; i = i - 1) {
+ if (!IS_INDEF (data[i])) {
+ if (data[i] != data[index])
+ break
+ nleft = nleft + 1
+ }
+ }
+ nright = 0
+ for (j = index + 1; i <= npoints; j = j + 1) {
+ if (!IS_INDEF (data[j])) {
+ if (data[j] != data[index])
+ break
+ nright = nright + 1
+ }
+ }
+
+ # Test for failure to be a local minima/axima
+ if (maxima == YES) {
+ if ((i == 0) && (j == npoints)) {
+ return (FALSE) # Data is constant
+ } else if (i == 0) {
+ if (data[j] > data[index])
+ return (FALSE) # Data increases to right
+ } else if (j == npoints) {
+ if (data[i] > data[index]) # Data increase to left
+ return (FALSE)
+ } else if ((data[i] > data[index]) || (data[j] > data[index])) {
+ return (FALSE) # Not a local maximum
+ } else if (!((nleft - nright == 0) || (nleft - nright == 1))) {
+ return (FALSE) # Not center of plateau
+ }
+ } else {
+ if ((i == 0) && (j == npoints)) {
+ return (FALSE) # Data is constant
+ } else if (i == 0) {
+ if (data[j] < data[index])
+ return (FALSE) # Data decreases to right
+ } else if (j == npoints) {
+ if (data[i] < data[index]) # Data decrease to left
+ return (FALSE)
+ } else if ((data[i] < data[index]) || (data[j] < data[index])) {
+ return (FALSE) # Not a local maximum
+ } else if (!((nleft - nright == 0) || (nleft - nright == 1))) {
+ return (FALSE) # Not center of plateau
+ }
+ }
+
+ # Point is a local minima/maxima
+ return (TRUE)
+end
+
+
+# FIND_THRESHOLD -- Find the peaks with positions satisfying threshold
+# and contrast constraints.
+#
+# The input is the data array, data, and the peak positions array, x.
+# The x array is resorted to the nthreshold peaks satisfying the constraints.
+# The corresponding nthreshold data values are returned the y array.
+# The number of peaks satisfying the constraints (nthreshold) is returned.
+
+int procedure find_threshold (data, x, y, npoints, contrast, threshold, debug)
+
+PIXEL data[ARB] # Input data values
+PIXEL x[npoints] # Input/Output peak positions
+PIXEL y[npoints] # Output peak data values
+int npoints # Number of peaks input
+real contrast # Contrast constraint
+real threshold # Threshold constraint
+bool debug # Print debugging information?
+
+int i, j, nthreshold
+PIXEL minval, maxval, lcut
+
+pointer dummy
+int maxima
+common /sort/ dummy, maxima
+
+begin
+ # Set the y array to be the values at the peak positions.
+ do i = 1, npoints {
+ j = x[i]
+ y[i] = data[j]
+ }
+
+ # Determine the min and max values of the peaks.
+ call alim$t (y, npoints, minval, maxval)
+
+ # Set the threshold based on the max of the absolute threshold and the
+ # contrast. Use arlt to set peaks below threshold to INDEF.
+
+ if (maxima == YES) {
+ lcut = max (PIXEL (threshold), PIXEL (contrast * maxval))
+ call arlt$t (y, npoints, lcut, INDEF)
+ } else {
+ lcut = threshold
+ call argt$t (y, npoints, lcut, INDEF)
+ }
+
+ if (debug) {
+ call printf (" Highest peak value = %g.\n")
+ call parg$t (maxval)
+ call printf (" Peak cutoff threshold = %g.\n")
+ call parg$t (lcut)
+ do i = 1, npoints {
+ if (IS_INDEF (y[i])) {
+ j = x[i]
+ call printf (
+ " Peak at column %d with value %g below threshold.\n")
+ call pargi (j)
+ call parg$t (data[j])
+ }
+ }
+ }
+
+ # Determine the number of acceptable peaks & resort the x and y arrays.
+ nthreshold = 0
+ do i = 1, npoints {
+ if (IS_INDEF (y[i]))
+ next
+ nthreshold = nthreshold + 1
+ x[nthreshold] = x[i]
+ y[nthreshold] = y[i]
+ }
+
+ if (debug) {
+ call printf (" Number of peaks above the threshold = %d.\n")
+ call pargi (nthreshold)
+ }
+
+ return (nthreshold)
+end
+
+# FIND_ISOLATED -- Flag peaks which are within separation of a peak
+# with a higher peak value.
+#
+# The peak positions, x, and their ranks, rank, are input.
+# The rank array contains the indices of the peak positions in order from
+# the highest peak value to the lowest peak value. Starting with
+# highest rank (rank[1]) all peaks of lower rank within separation
+# are marked by setting their positions to INDEF. The number of
+# unflaged peaks is returned.
+
+int procedure find_isolated (x, rank, npoints, separation, debug)
+
+# Procedure parameters:
+PIXEL x[npoints] # Positions of points
+int rank[npoints] # Rank of peaks
+int npoints # Number of peaks
+int separation # Minimum allowed separation
+bool debug # Print diagnostic information
+
+int i, j
+int nisolated
+
+begin
+ # Eliminate close neighbors. The eliminated
+ # peaks are marked by setting their positions to INDEF.
+ nisolated = 0
+ do i = 1, npoints {
+ if (IS_INDEF (x[rank[i]]))
+ next
+ nisolated = nisolated + 1
+ do j = i + 1, npoints {
+ if (IS_INDEF (x[rank[j]]))
+ next
+ if (abs (x[rank[i]] - x[rank[j]]) < separation) {
+ if (debug) {
+ call printf (
+ " Peak at column %d too near peak at column %d.\n")
+ call pargi (int (x[rank[j]]))
+ call pargi (int (x[rank[i]]))
+ }
+ x[rank[j]] = INDEF
+ }
+ }
+ }
+
+ if (debug) {
+ call printf (" Number of peaks separated by %d pixels = %d.\n")
+ call pargi (separation)
+ call pargi (nisolated)
+ }
+
+ # Return number of isolated peaks.
+ return (nisolated)
+end
+
+
+# FIND_NMAX -- Select up to the nmax highest ranked peaks.
+#
+# The data values, data, peak positions, x, and their ranks, rank, are input.
+# The data values are used only in printing debugging information.
+# Peak positions previously eliminated are flaged by the value INDEF.
+# The rank array contains the indices to the peak positions in order from
+# the highest peak value to the lowest peak value.
+# First all but the nmax highest ranked peaks (which have not been previously
+# eliminated) are eliminated by marking their positions with the value INDEF.
+# Then the remaining peaks are resorted to contain only the unflaged
+# peaks and the number of such peaks is returned.
+
+int procedure find_nmax (data, x, rank, npoints, nmax, debug)
+
+PIXEL data[ARB] # Input data values
+PIXEL x[npoints] # Peak positions
+int rank[npoints] # Ranks of peaks
+int npoints # Number of input peaks
+int nmax # Max number of peaks to be selected
+bool debug # Print debugging information?
+
+int i, j, npeaks
+
+begin
+ # Only mark peaks to reject if the number peaks is greater than nmax.
+ if (nmax < npoints) {
+ npeaks = 0
+ do i = 1, npoints {
+ if (IS_INDEF (x[rank[i]]))
+ next
+ npeaks = npeaks + 1
+ if (npeaks > nmax) {
+ if (debug) {
+ j = x[rank[i]]
+ call printf (
+ " Reject peak at column %d with rank %d and value %g.\n")
+ call pargi (j)
+ call pargi (i)
+ call parg$t (data[j])
+ }
+ x[rank[i]] = INDEF
+ }
+ }
+ }
+
+ # Eliminate INDEF points and determine the number of spectra found.
+ npeaks = 0
+ do i = 1, npoints {
+ if (IS_INDEF (x[i]))
+ next
+ npeaks = npeaks + 1
+ x[npeaks] = x[i]
+ }
+
+ return (npeaks)
+end
+
+
+# COMPARE -- Compare procedure for sort used in FIND_PEAKS.
+# Larger values are indexed first. INDEF values are indexed last.
+
+int procedure compare (index1, index2)
+
+# Procedure parameters:
+int index1 # Comparison index
+int index2 # Comparison index
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # INDEF points are considered to be smallest/largest possible values.
+ if (maxima == YES) {
+ if (IS_INDEF (Mem$t[y - 1 + index1]))
+ return (1)
+ else if (IS_INDEF (Mem$t[y - 1 + index2]))
+ return (-1)
+ else if (Mem$t[y - 1 + index1] < Mem$t[y - 1 + index2])
+ return (1)
+ else if (Mem$t[y - 1 + index1] > Mem$t[y - 1 + index2])
+ return (-1)
+ else
+ return (0)
+ } else {
+ if (IS_INDEF (Mem$t[y - 1 + index1]))
+ return (-1)
+ else if (IS_INDEF (Mem$t[y - 1 + index2]))
+ return (1)
+ else if (Mem$t[y - 1 + index1] < Mem$t[y - 1 + index2])
+ return (-1)
+ else if (Mem$t[y - 1 + index1] > Mem$t[y - 1 + index2])
+ return (1)
+ else
+ return (0)
+ }
+end
+$endfor
diff --git a/noao/rv/rvidlines/peaks.x b/noao/rv/rvidlines/peaks.x
new file mode 100644
index 00000000..8ee27f51
--- /dev/null
+++ b/noao/rv/rvidlines/peaks.x
@@ -0,0 +1,446 @@
+# PEAKS -- The following procedures are general numerical functions
+# dealing with finding peaks in a data array.
+#
+# FIND_PEAKS Find the peaks in the data array.
+# FIND_LOCAL_MAXIMA Find the local minima/maxima in the data array.
+# IS_LOCAL_MAX Test a point to determine if it is a local min/max.
+# FIND_THRESHOLD Find the peaks with positions satisfying threshold
+# and contrast constraints.
+# FIND_ISOLATED Flag peaks which are within separation of a peak
+# with a higher peak value.
+# FIND_NMAX Select up to the nmax highest ranked peaks.
+# COMPARE Compare procedure for sort used in FIND_PEAKS.
+
+# FIND_PEAKS -- Find the peaks in the data array.
+#
+# The peaks are found using the following algorithm:
+#
+# 1. Find the local minima/maxima.
+# 2. Reject peaks below the threshold.
+# 3. Determine the ranks of the remaining peaks.
+# 4. Flag weaker peaks within separation of a stronger peak.
+# 5. Accept at most the nmax strongest peaks.
+#
+# Indefinite points are ignored. The peak positions are returned in the
+# array x.
+
+
+int procedure find_peaks (data, x, npoints, contrast, separation, edge, nmax,
+ threshold, debug)
+
+# Procedure parameters:
+real data[npoints] # Input data array
+real x[npoints] # Output peak position array
+int npoints # Number of data points
+real contrast # Maximum contrast between strongest and weakest
+int separation # Minimum separation between peaks
+int edge # Minimum distance from the edge
+int nmax # Maximum number of peaks to be returned
+real threshold # Minimum threshold level for peaks
+bool debug # Print diagnostic information?
+
+int i, j
+int nlmax, nthreshold, nisolated, npeaks
+pointer sp, rank
+
+int find_local_maxima(), find_threshold(), find_isolated(), find_nmax()
+int compare()
+
+extern compare()
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # Find the local minima/maxima in data and put column positions in x..
+ if (contrast < 0.)
+ maxima = NO
+ else
+ maxima = YES
+
+ nlmax = find_local_maxima (data, x, npoints, debug)
+
+ # Reject local minima/maxima near the edge.
+ if (edge > 0) {
+ j = 0
+ do i = 1, nlmax {
+ if ((x[i] > edge) && (x[i] <= npoints - edge)) {
+ j = j + 1
+ x[j] = x[i]
+ }
+ }
+ nlmax = j
+ }
+
+ # Allocate a working array y.
+ call smark (sp)
+ call salloc (y, npoints, TY_REAL)
+
+ # Reject the local minima/maxima which do not satisfy the thresholds.
+ # The array y is set to the peak values of the remaining peaks.
+ nthreshold = find_threshold (data, x, Memr[y], nlmax,
+ contrast, threshold, debug)
+
+ # Rank the peaks by peak value.
+ call salloc (rank, nthreshold, TY_INT)
+ do i = 1, nthreshold
+ Memi[rank + i - 1] = i
+ call qsort (Memi[rank], nthreshold, compare)
+
+ # Reject the weaker peaks within sep of a stronger peak.
+ nisolated = find_isolated (x, Memi[rank], nthreshold, separation,
+ debug)
+
+ # Select the strongest nmax peaks.
+ npeaks = find_nmax (data, x, Memi[rank], nthreshold, nmax, debug)
+
+ call sfree (sp)
+ return (npeaks)
+end
+
+
+# FIND_LOCAL_MAXIMA -- Find the local minima/maxima in the data array.
+#
+# A data array is input and the local minima/maxima positions array is output.
+# The number of local minima/maxima found is returned.
+
+int procedure find_local_maxima (data, x, npoints, debug)
+
+real data[npoints] # Input data array
+real x[npoints] # Output local min/max positions array
+int npoints # Number of input points
+bool debug # Print debugging information?
+
+int i, nlmax
+
+bool is_local_max()
+
+begin
+ nlmax = 0
+ do i = 1, npoints {
+ if (is_local_max (i, data, npoints)) {
+ nlmax = nlmax + 1
+ x[nlmax] = i
+ }
+ }
+
+ if (debug) {
+ call printf (" Number of local minima/maxima found = %d.\n")
+ call pargi (nlmax)
+ }
+
+ return (nlmax)
+end
+
+
+# IS_LOCAL_MAX -- Test a point to determine if it is a local minima/maximum.
+#
+# Indefinite points are ignored.
+
+bool procedure is_local_max (index, data, npoints)
+
+# Procedure parameters:
+int index # Index to test for local minimum/maximum
+real data[npoints] # Data values
+int npoints # Number of points in the data vector
+
+int i, j, nright, nleft
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # INDEF points cannot be local minima/axima.
+ if (IS_INDEFR (data[index]))
+ return (FALSE)
+
+ # Find the left and right indices where data values change and the
+ # number of points with the same value. Ignore INDEF points.
+ nleft = 0
+ for (i = index - 1; i >= 1; i = i - 1) {
+ if (!IS_INDEFR (data[i])) {
+ if (data[i] != data[index])
+ break
+ nleft = nleft + 1
+ }
+ }
+ nright = 0
+ for (j = index + 1; i <= npoints; j = j + 1) {
+ if (!IS_INDEFR (data[j])) {
+ if (data[j] != data[index])
+ break
+ nright = nright + 1
+ }
+ }
+
+ # Test for failure to be a local minima/axima
+ if (maxima == YES) {
+ if ((i == 0) && (j == npoints)) {
+ return (FALSE) # Data is constant
+ } else if (i == 0) {
+ if (data[j] > data[index])
+ return (FALSE) # Data increases to right
+ } else if (j == npoints) {
+ if (data[i] > data[index]) # Data increase to left
+ return (FALSE)
+ } else if ((data[i] > data[index]) || (data[j] > data[index])) {
+ return (FALSE) # Not a local maximum
+ } else if (!((nleft - nright == 0) || (nleft - nright == 1))) {
+ return (FALSE) # Not center of plateau
+ }
+ } else {
+ if ((i == 0) && (j == npoints)) {
+ return (FALSE) # Data is constant
+ } else if (i == 0) {
+ if (data[j] < data[index])
+ return (FALSE) # Data decreases to right
+ } else if (j == npoints) {
+ if (data[i] < data[index]) # Data decrease to left
+ return (FALSE)
+ } else if ((data[i] < data[index]) || (data[j] < data[index])) {
+ return (FALSE) # Not a local maximum
+ } else if (!((nleft - nright == 0) || (nleft - nright == 1))) {
+ return (FALSE) # Not center of plateau
+ }
+ }
+
+ # Point is a local minima/maxima
+ return (TRUE)
+end
+
+
+# FIND_THRESHOLD -- Find the peaks with positions satisfying threshold
+# and contrast constraints.
+#
+# The input is the data array, data, and the peak positions array, x.
+# The x array is resorted to the nthreshold peaks satisfying the constraints.
+# The corresponding nthreshold data values are returned the y array.
+# The number of peaks satisfying the constraints (nthreshold) is returned.
+
+int procedure find_threshold (data, x, y, npoints, contrast, threshold, debug)
+
+real data[ARB] # Input data values
+real x[npoints] # Input/Output peak positions
+real y[npoints] # Output peak data values
+int npoints # Number of peaks input
+real contrast # Contrast constraint
+real threshold # Threshold constraint
+bool debug # Print debugging information?
+
+int i, j, nthreshold
+real minval, maxval, lcut
+
+pointer dummy
+int maxima
+common /sort/ dummy, maxima
+
+begin
+ # Set the y array to be the values at the peak positions.
+ do i = 1, npoints {
+ j = x[i]
+ y[i] = data[j]
+ }
+
+ # Determine the min and max values of the peaks.
+ call alimr (y, npoints, minval, maxval)
+
+ # Set the threshold based on the max of the absolute threshold and the
+ # contrast. Use arlt to set peaks below threshold to INDEF.
+
+ if (maxima == YES) {
+ lcut = max (real (threshold), real (contrast * maxval))
+ call arltr (y, npoints, lcut, INDEFR)
+ } else {
+ lcut = threshold
+ call argtr (y, npoints, lcut, INDEFR)
+ }
+
+ if (debug) {
+ call printf (" Highest peak value = %g.\n")
+ call pargr (maxval)
+ call printf (" Peak cutoff threshold = %g.\n")
+ call pargr (lcut)
+ do i = 1, npoints {
+ if (IS_INDEFR (y[i])) {
+ j = x[i]
+ call printf (
+ " Peak at column %d with value %g below threshold.\n")
+ call pargi (j)
+ call pargr (data[j])
+ }
+ }
+ }
+
+ # Determine the number of acceptable peaks & resort the x and y arrays.
+ nthreshold = 0
+ do i = 1, npoints {
+ if (IS_INDEFR (y[i]))
+ next
+ nthreshold = nthreshold + 1
+ x[nthreshold] = x[i]
+ y[nthreshold] = y[i]
+ }
+
+ if (debug) {
+ call printf (" Number of peaks above the threshold = %d.\n")
+ call pargi (nthreshold)
+ }
+
+ return (nthreshold)
+end
+
+# FIND_ISOLATED -- Flag peaks which are within separation of a peak
+# with a higher peak value.
+#
+# The peak positions, x, and their ranks, rank, are input.
+# The rank array contains the indices of the peak positions in order from
+# the highest peak value to the lowest peak value. Starting with
+# highest rank (rank[1]) all peaks of lower rank within separation
+# are marked by setting their positions to INDEF. The number of
+# unflaged peaks is returned.
+
+int procedure find_isolated (x, rank, npoints, separation, debug)
+
+# Procedure parameters:
+real x[npoints] # Positions of points
+int rank[npoints] # Rank of peaks
+int npoints # Number of peaks
+int separation # Minimum allowed separation
+bool debug # Print diagnostic information
+
+int i, j
+int nisolated
+
+begin
+ # Eliminate close neighbors. The eliminated
+ # peaks are marked by setting their positions to INDEF.
+ nisolated = 0
+ do i = 1, npoints {
+ if (IS_INDEFR (x[rank[i]]))
+ next
+ nisolated = nisolated + 1
+ do j = i + 1, npoints {
+ if (IS_INDEFR (x[rank[j]]))
+ next
+ if (abs (x[rank[i]] - x[rank[j]]) < separation) {
+ if (debug) {
+ call printf (
+ " Peak at column %d too near peak at column %d.\n")
+ call pargi (int (x[rank[j]]))
+ call pargi (int (x[rank[i]]))
+ }
+ x[rank[j]] = INDEFR
+ }
+ }
+ }
+
+ if (debug) {
+ call printf (" Number of peaks separated by %d pixels = %d.\n")
+ call pargi (separation)
+ call pargi (nisolated)
+ }
+
+ # Return number of isolated peaks.
+ return (nisolated)
+end
+
+
+# FIND_NMAX -- Select up to the nmax highest ranked peaks.
+#
+# The data values, data, peak positions, x, and their ranks, rank, are input.
+# The data values are used only in printing debugging information.
+# Peak positions previously eliminated are flaged by the value INDEF.
+# The rank array contains the indices to the peak positions in order from
+# the highest peak value to the lowest peak value.
+# First all but the nmax highest ranked peaks (which have not been previously
+# eliminated) are eliminated by marking their positions with the value INDEF.
+# Then the remaining peaks are resorted to contain only the unflaged
+# peaks and the number of such peaks is returned.
+
+int procedure find_nmax (data, x, rank, npoints, nmax, debug)
+
+real data[ARB] # Input data values
+real x[npoints] # Peak positions
+int rank[npoints] # Ranks of peaks
+int npoints # Number of input peaks
+int nmax # Max number of peaks to be selected
+bool debug # Print debugging information?
+
+int i, j, npeaks
+
+begin
+ # Only mark peaks to reject if the number peaks is greater than nmax.
+ if (nmax < npoints) {
+ npeaks = 0
+ do i = 1, npoints {
+ if (IS_INDEFR (x[rank[i]]))
+ next
+ npeaks = npeaks + 1
+ if (npeaks > nmax) {
+ if (debug) {
+ j = x[rank[i]]
+ call printf (
+ " Reject peak at column %d with rank %d and value %g.\n")
+ call pargi (j)
+ call pargi (i)
+ call pargr (data[j])
+ }
+ x[rank[i]] = INDEFR
+ }
+ }
+ }
+
+ # Eliminate INDEF points and determine the number of spectra found.
+ npeaks = 0
+ do i = 1, npoints {
+ if (IS_INDEFR (x[i]))
+ next
+ npeaks = npeaks + 1
+ x[npeaks] = x[i]
+ }
+
+ return (npeaks)
+end
+
+
+# COMPARE -- Compare procedure for sort used in FIND_PEAKS.
+# Larger values are indexed first. INDEF values are indexed last.
+
+int procedure compare (index1, index2)
+
+# Procedure parameters:
+int index1 # Comparison index
+int index2 # Comparison index
+
+pointer y
+int maxima
+common /sort/ y, maxima
+
+begin
+ # INDEF points are considered to be smallest/largest possible values.
+ if (maxima == YES) {
+ if (IS_INDEFR (Memr[y - 1 + index1]))
+ return (1)
+ else if (IS_INDEFR (Memr[y - 1 + index2]))
+ return (-1)
+ else if (Memr[y - 1 + index1] < Memr[y - 1 + index2])
+ return (1)
+ else if (Memr[y - 1 + index1] > Memr[y - 1 + index2])
+ return (-1)
+ else
+ return (0)
+ } else {
+ if (IS_INDEFR (Memr[y - 1 + index1]))
+ return (-1)
+ else if (IS_INDEFR (Memr[y - 1 + index2]))
+ return (1)
+ else if (Memr[y - 1 + index1] < Memr[y - 1 + index2])
+ return (-1)
+ else if (Memr[y - 1 + index1] > Memr[y - 1 + index2])
+ return (1)
+ else
+ return (0)
+ }
+end
diff --git a/noao/rv/rvidlines/reidentify.x b/noao/rv/rvidlines/reidentify.x
new file mode 100644
index 00000000..ed8133b2
--- /dev/null
+++ b/noao/rv/rvidlines/reidentify.x
@@ -0,0 +1,609 @@
+include <error.h>
+include <imhdr.h>
+include <gset.h>
+include "identify.h"
+
+define HELP "noao$onedspec/identify/identify.key"
+define RVHELP "noao$rv/rvidlines/rvidlines.key"
+define ICFITHELP "noao$lib/scr/idicgfit.key"
+define PROMPT "identify options"
+
+define PAN 1 # Pan graph
+define ZOOM 2 # Zoom graph
+
+# REIDENTIFY -- Reidentify features in an image.
+
+procedure reidentify (id)
+
+pointer id # ID pointer
+
+real wx, wy, wx2, wy2
+int wcs, key
+char cmd[SZ_LINE]
+
+char newimage[SZ_FNAME]
+int i, j, last, all, prfeature, nfeatures1, npeaks, ftype
+double pix, fit, user, shift, pix_shift, z_shift, xg[10]
+pointer peaks, label
+
+int clgcur(), scan(), nscan(), find_peaks(), errcode()
+double id_center(), fit_to_pix(), id_fitpt()
+double id_shift(), id_rms(), id_zshiftd(), id_zval()
+errchk id_graph()
+
+define newim_ 10
+define newkey_ 20
+define beep_ 99
+
+begin
+ # Initialize.
+ ID_GTYPE(id) = PAN
+ all = 0
+ last = ID_CURRENT(id)
+ newimage[1] = EOS
+ ID_REFIT(id) = NO
+ ID_NEWFEATURES(id) = NO
+ ID_NEWCV(id) = NO
+ wy = INDEF
+ key = 'r'
+
+ repeat {
+ prfeature = YES
+ if (all != 0)
+ all = mod (all + 1, 3)
+
+ switch (key) {
+ case '?': # Print help
+ if (ID_TASK(id) == IDENTIFY)
+ call gpagefile (ID_GP(id), HELP, PROMPT)
+ else
+ call gpagefile (ID_GP(id), RVHELP, PROMPT)
+ case ':': # Process colon commands
+ if (cmd[1] == '/')
+ call gt_colon (cmd, ID_GP(id), ID_GT(id), ID_NEWGRAPH(id))
+ else
+ call id_colon (id, cmd, newimage, prfeature)
+ case ' ': # Go to current feature
+ case '.': # Go to nearest feature
+ if (ID_NFEATURES(id) == 0)
+ goto beep_
+ call id_nearest (id, double (wx))
+ case '-': # Go to previous feature
+ if (ID_CURRENT(id) == 1)
+ goto beep_
+ ID_CURRENT(id) = ID_CURRENT(id) - 1
+ case '+', 'n': # Go to next feature
+ if (ID_CURRENT(id) == ID_NFEATURES(id))
+ goto beep_
+ ID_CURRENT(id) = ID_CURRENT(id) + 1
+ case 'a': # Set all flag for next key
+ all = 1
+ case 'b': # Mark blended features
+ i = 1
+ fit = wx
+ xg[i] = fit_to_pix (id, fit)
+ call printf ("mark other components (exit with 'q'):")
+ while (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)!=EOF) {
+ if (key == 'q')
+ break
+ i = i + 1
+ fit = wx
+ xg[i] = fit_to_pix (id, fit)
+ if (i == 10)
+ break
+ }
+
+ call printf ("mark two background points: ")
+ j = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+ j = clgcur ("cursor", wx2, wy2, wcs, key, cmd, SZ_LINE)
+ wx = fit_to_pix (id, double (wx))
+ wx2 = fit_to_pix (id, double (wx2))
+
+ switch (ID_FTYPE(id)) {
+ case EMISSION:
+ ftype = GEMISSION
+ case ABSORPTION:
+ ftype = GABSORPTION
+ default:
+ ftype = ID_FTYPE(id)
+ }
+ iferr (call id_gcenter (id, xg, i, wx, wy, wx2, wy2,
+ ID_FWIDTH(id), ftype, YES)) {
+ call erract (EA_WARN)
+ prfeature = NO
+ goto newkey_
+ }
+
+ do j = 1, i {
+ pix = xg[j]
+ fit = id_fitpt (id, pix)
+ user = fit
+ call id_newfeature (id, pix, fit, user, 1.0D0,
+ ID_FWIDTH(id), ftype, NULL)
+ USER(id,ID_CURRENT(id)) = INDEFD
+ call id_match (id, FIT(id,ID_CURRENT(id)),
+ USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ call id_mark (id, ID_CURRENT(id))
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("%10.8g ")
+ call pargd (
+ id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ fit = id_zshiftd (id, user, 1)
+ call id_match (id, fit, USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd,
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ }
+ case 'c': # Recenter features
+ if (all != 0) {
+ for (i = 1; i <= ID_NFEATURES(id); i = i + 1) {
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, i)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ FWIDTH(id,i) = ID_FWIDTH(id)
+ PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i),
+ FTYPE(id,i), NO)
+ if (!IS_INDEFD (PIX(id,i))) {
+ FIT(id,i) = id_fitpt (id, PIX(id,i))
+ call id_mark (id, i)
+ } else {
+ call id_delete (id, i)
+ i = i - 1
+ }
+ }
+ ID_NEWFEATURES(id) = YES
+ } else {
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call id_nearest (id, double (wx))
+ pix = PIX(id,ID_CURRENT(id))
+ pix = id_center (id, pix, 1, ID_FWIDTH(id),
+ FTYPE(id,ID_CURRENT(id)), NO)
+ if (!IS_INDEFD (pix)) {
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ PIX(id,ID_CURRENT(id)) = pix
+ FWIDTH(id,ID_CURRENT(id)) = ID_FWIDTH(id)
+ FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_mark (id, ID_CURRENT(id))
+ ID_NEWFEATURES(id) = YES
+ } else {
+ call printf ("Centering failed\n")
+ prfeature = NO
+ }
+ }
+ case 'd': # Delete features
+ if (all != 0) {
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ ID_NEWFEATURES(id) = YES
+ ID_NEWGRAPH(id) = YES
+ } else {
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call id_nearest (id, double (wx))
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_delete (id, ID_CURRENT(id))
+ ID_CURRENT(id) = min (ID_NFEATURES(id), ID_CURRENT(id))
+ last = 0
+ }
+ case 'f': # Fit dispersion function
+ if (ID_TASK(id) == IDENTIFY) {
+ call id_dofit (id, YES)
+ } else {
+ call id_velocity (id, YES)
+ prfeature = NO
+ }
+ case 'g': # Fit shift
+ if (ID_TASK(id) == RVIDLINES)
+ goto beep_
+
+ call id_doshift (id, YES)
+ prfeature = NO
+ case 'i': # Initialize
+ call dcvfree (ID_CV(id))
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+ ID_NEWCV(id) = YES
+ ID_NFEATURES(id) = 0
+ ID_CURRENT(id) = 0
+ ID_NEWFEATURES(id) = YES
+ ID_NEWGRAPH(id) = YES
+ case 'j', 'k', 'o':
+ call printf ("Command not available in REIDENTIFY")
+ prfeature = NO
+ case 'l': # Find features from line list
+ if (ID_TASK(id) == IDENTIFY) {
+ if (ID_NFEATURES(id) >= 2)
+ call id_dofit (id, NO)
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata(id))
+ ;
+ call id_fitfeatures(id)
+ ID_NEWCV(id) = NO
+ }
+ call id_linelist (id)
+ if (ID_NEWFEATURES(id) == YES)
+ ID_REFIT(id) = YES
+ } else {
+ call id_velocity (id, NO)
+ call id_linelist(id)
+ if (ID_NEWFEATURES(id) == YES) {
+ call id_velocity (id, NO)
+ ID_NEWGRAPH(id) = YES
+ }
+ }
+ case 'm': # Mark new feature
+ fit = wx
+ pix = fit_to_pix (id, fit)
+ pix = id_center (id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), YES)
+ if (IS_INDEFD (pix))
+ goto beep_
+ fit = id_fitpt (id, pix)
+ user = fit
+ call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id),
+ ID_FTYPE(id), NULL)
+ USER(id,ID_CURRENT(id)) = INDEFD
+ call id_match (id, FIT(id,ID_CURRENT(id)),
+ USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ call id_mark (id, ID_CURRENT(id))
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("[%10.8g] ")
+ call pargd (
+ id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ call id_match (id, user, USER(id,ID_CURRENT(id)),
+ Memi[ID_LABEL(id)+ID_CURRENT(id)-1],
+ ID_MATCH(id))
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ case 'p': # Switch to pan mode
+ if (ID_GTYPE(id) != PAN) {
+ ID_GTYPE(id) = PAN
+ ID_NEWGRAPH(id) = YES
+ }
+ case 'q': # Exit loop
+ break
+ case 'r': # Redraw the graph
+ ID_NEWGRAPH(id) = YES
+ case 's', 'x': # Shift or correlate features
+ if (ID_TASK(id) == RVIDLINES)
+ goto beep_
+
+ # Get coordinate shift.
+ switch (key) {
+ case 's':
+ call printf ("User coordinate (%10.8g): ")
+ call pargr (wx)
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ if (nscan() == 1)
+ shift = wx - user
+ } else
+ shift = 0.
+ case 'x':
+ if (ID_NFEATURES(id) > 5)
+ shift = id_shift (id)
+ else
+ goto beep_
+ }
+
+ ID_NEWFEATURES(id) = YES
+ ID_NEWCV(id) = YES
+ ID_NEWGRAPH(id) = YES
+ prfeature = NO
+
+ if (ID_NFEATURES(id) < 1) {
+ call printf ("User coordinate shift=%5f\n")
+ call pargd (shift)
+ ID_SHIFT(id) = ID_SHIFT(id) + shift
+ goto newkey_
+ }
+
+ # Recenter features.
+ pix_shift = 0.
+ z_shift = 0.
+ nfeatures1 = ID_NFEATURES(id)
+
+ j = 0.
+ do i = 1, ID_NFEATURES(id) {
+ pix = fit_to_pix (id, FIT(id,i) + shift)
+ pix = id_center (id, pix, 1, FWIDTH(id,i), FTYPE(id,i), NO)
+ if (IS_INDEFD (pix)) {
+ if (ID_CURRENT(id) == i)
+ ID_CURRENT(id) = i + 1
+ next
+ }
+ fit = id_fitpt (id, pix)
+
+ pix_shift = pix_shift + pix - PIX(id,i)
+ if (FIT(id,i) != 0.)
+ z_shift = z_shift + id_zval (id, fit, FIT(id,i))
+
+ j = j + 1
+ PIX(id,j) = pix
+ 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)
+ if (ID_CURRENT(id) == i)
+ ID_CURRENT(id) = j
+ }
+ if (j != ID_NFEATURES(id)) {
+ ID_NFEATURES(id) = j
+ ID_CURRENT(id) = min (ID_CURRENT(id), ID_NFEATURES(id))
+ }
+
+ if (ID_NFEATURES(id) < 1) {
+ call printf ("User coordinate shift=%5f")
+ call pargd (shift)
+ call printf (", No features found during recentering\n")
+ ID_SHIFT(id) = ID_SHIFT(id) + shift
+ goto newkey_
+ }
+
+ # Adjust shift.
+ pix = ID_SHIFT(id)
+ call id_doshift (id, NO)
+ call id_fitfeatures (id)
+
+ # Print results.
+ call printf ("Recentered=%d/%d")
+ call pargi (ID_NFEATURES(id))
+ call pargi (nfeatures1)
+ call printf (
+ ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g\n")
+ call pargd (pix_shift / ID_NFEATURES(id))
+ call pargd (pix - ID_SHIFT(id))
+ call pargd (z_shift / ID_NFEATURES(id))
+ call pargd (id_rms(id))
+ case 't': # Move the current feature
+ if (ID_CURRENT(id) < 1)
+ goto beep_
+ pix = fit_to_pix (id, double (wx))
+ call gseti (ID_GP(id), G_PLTYPE, 0)
+ call id_mark (id, ID_CURRENT(id))
+ PIX(id,ID_CURRENT(id)) = pix
+ FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix)
+ call gseti (ID_GP(id), G_PLTYPE, 1)
+ call id_mark (id, ID_CURRENT(id))
+ ID_NEWFEATURES(id) = YES
+ case 'u': # Set user coordinate
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ call printf ("%10.2f %10.8g ")
+ call pargd (PIX(id,ID_CURRENT(id)))
+ call pargd (FIT(id,ID_CURRENT(id)))
+ if (ID_REDSHIFT(id) != 0.) {
+ call printf ("[%10.8g] ")
+ call pargd (id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0))
+ }
+ call printf ("(%10.8g %s): ")
+ call pargd (USER(id,ID_CURRENT(id)))
+ label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ call gargwrd (cmd, SZ_LINE)
+ i = nscan()
+ if (i > 0) {
+ USER(id,ID_CURRENT(id)) = user
+ ID_NEWFEATURES(id) = YES
+ }
+ if (i > 1) {
+ call reset_scan ()
+ call gargd (user)
+ call gargstr (cmd, SZ_LINE)
+ call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1])
+ }
+ }
+ call printf ("Weight (%g): ")
+ call pargd (WTS(id,ID_CURRENT(id)))
+ call flush (STDOUT)
+ if (scan() != EOF) {
+ call gargd (user)
+ if (nscan() > 0)
+ WTS(id,ID_CURRENT(id)) = user
+ }
+ case 'w': # Window graph
+ call gt_window (ID_GT(id), ID_GP(id), "cursor", ID_NEWGRAPH(id))
+ case 'y': # Find peaks
+ call malloc (peaks, ID_NPTS(id), TY_REAL)
+ npeaks = find_peaks (IMDATA(id,1), Memr[peaks], ID_NPTS(id), 0.,
+ int (ID_MINSEP(id)), 0, ID_MAXFEATURES(id), 0., false)
+ for (j = 1; j <= ID_NFEATURES(id); j = j + 1) {
+ for (i = 1; i <= npeaks; i = i + 1) {
+ if (!IS_INDEF (Memr[peaks+i-1])) {
+ pix = Memr[peaks+i-1]
+ if (abs (pix - PIX(id,j)) < ID_MINSEP(id))
+ Memr[peaks+i-1] = INDEF
+ }
+ }
+ }
+ for (i = 1; i <= npeaks; i = i + 1) {
+ if (IS_INDEF(Memr[peaks+i-1]))
+ next
+ pix = Memr[peaks+i-1]
+ pix = id_center(id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), NO)
+ if (IS_INDEFD (pix))
+ next
+ fit = id_fitpt (id, pix)
+ user = INDEFD
+ call id_match (id, fit, user, label, ID_MATCH(id))
+ call id_newfeature (id, pix, fit, user, 1.0D0,
+ ID_FWIDTH(id), ID_FTYPE(id), label)
+ call id_mark (id, ID_CURRENT(id))
+ }
+ call mfree (peaks, TY_REAL)
+ case 'z': # Go to zoom mode
+ if (ID_NFEATURES(id) < 1)
+ goto beep_
+ if (ID_GTYPE(id) == PAN)
+ ID_NEWGRAPH(id) = YES
+ ID_GTYPE(id) = ZOOM
+ call id_nearest (id, double (wx))
+ case 'I':
+ call fatal (0, "Interrupt")
+ default:
+beep_ call printf ("\007")
+ }
+
+newkey_
+ # Set update flag if anything has changed.
+ if ((ID_NEWFEATURES(id) == YES) || (ID_NEWCV(id) == YES))
+ ID_NEWDBENTRY(id) = YES
+
+ # If a new image exit loop, update database, and start over.
+ if (newimage[1] != EOS) {
+ call printf ("Can't change image in REIDENTIFY")
+ newimage[1] = EOS
+ prfeature = NO
+ }
+
+ # Refit dispersion function
+ if (ID_REFIT(id) == YES) {
+ call id_dofit (id, NO)
+ ID_REFIT(id) = NO
+ }
+
+ # If there is a new dispersion solution evaluate the coordinates
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata (id))
+ ;
+ call id_fitfeatures (id)
+ ID_NEWCV(id) = NO
+ }
+
+ # Draw new graph in zoom mode if current feature has changed.
+ if ((ID_GTYPE(id) == ZOOM) && (last != ID_CURRENT(id)))
+ ID_NEWGRAPH(id) = YES
+
+ # Draw new graph.
+ if (ID_NEWGRAPH(id) == YES) {
+ call id_graph (id, ID_GTYPE(id))
+ ID_NEWGRAPH(id) = NO
+ }
+
+ # Set cursor and print status of current feature (unless canceled).
+ if (ID_CURRENT(id) > 0) {
+ if (IS_INDEF (wy)) {
+ i = max (1, min (ID_NPTS(id), int (PIX(id,ID_CURRENT(id)))))
+ wy = IMDATA(id,i)
+ }
+
+ call gscur (ID_GP(id), real (FIT(id,ID_CURRENT(id))), wy)
+ if (errcode() == OK && prfeature == YES) {
+ i = ID_CURRENT(id)
+ pix = PIX(id,i)
+ fit = FIT(id,i)
+ user = USER(id,i)
+ if (ID_TASK(id) == IDENTIFY) {
+ if (IS_INDEFD(user))
+ shift = INDEF
+ else
+ shift = fit - user
+ call printf ("%10.2f %10.8g %10.8g %10.3g %s")
+ call pargd (pix)
+ call pargd (fit)
+ call pargd (user)
+ call pargd (shift)
+ label = Memi[ID_LABEL(id)+i-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ } else {
+ if (IS_INDEFD(user))
+ shift = INDEF
+ else {
+ shift = id_zval (id, fit, user) - ID_REDSHIFT(id)
+ if (abs (shift) < 0.01)
+ shift = shift * VLIGHT
+ }
+ call printf ("%10.2f %10.8g %10.8g %10.8g %10.4g %s")
+ call pargd (pix)
+ call pargd (fit)
+ call pargd (id_zshiftd (id, fit, 0))
+ call pargd (user)
+ if (IS_INDEFD(user))
+ call pargd (INDEFD)
+ else
+ call pargd (shift)
+ label = Memi[ID_LABEL(id)+i-1]
+ if (label != NULL)
+ call pargstr (Memc[label])
+ else
+ call pargstr ("")
+ }
+ }
+ }
+
+ # Print delayed error message
+ if (errcode() != OK)
+ call erract (EA_WARN)
+
+ last = ID_CURRENT(id)
+ } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+end
diff --git a/noao/rv/rvidlines/rvidlines.key b/noao/rv/rvidlines/rvidlines.key
new file mode 100644
index 00000000..cfccacde
--- /dev/null
+++ b/noao/rv/rvidlines/rvidlines.key
@@ -0,0 +1,100 @@
+ RVIDLINES HELP
+
+
+STATUS LINE
+
+The status line gives the pixel position, observed wavelength,
+rest wavelength, user wavelength, velocity residual, and optional
+line identification:
+
+ pixel observed rest user Vresidual [identification]
+
+
+CURSOR KEY SUMMARY
+
+? Help l Match list w Window graph
+a Affect all features m Mark feature y Find peaks
+b Deblend n Next feature z Zoom graph
+c Center feature(s) o Go to line + Next feature
+d Delete feature(s) p Pan graph - Previous feature
+f Fit velocity q Quit . Nearest feature
+i Initialize r Redraw graph I Interrupt
+j Preceding line t Reset position
+k Next line u Enter coordinate
+
+
+COLON COMMAND SUMMARY
+
+:show [file] :features [file] :coordlist [file]
+:cradius [value] :threshold [value] :database [file]
+:ftype [type] :fwidth [value] :image [image]
+:labels [type] :match [value] :maxfeatures [value]
+:minsep [value] :read [image [ap]] :write [image [ap]]
+:zwidth [value] :add [image [ap]]
+
+
+CURSOR KEYS
+
+? Clear the screen and print menu of options
+a Apply next (c)enter or (d)elete operation to (a)ll features
+b Mark and de(b)lend features by Gaussian fitting
+c (C)enter the feature nearest the cursor
+d (D)elete the feature nearest the cursor
+f (F)it a redshift and velocity from the fitted and user coordinates
+i (I)nitialize (delete features and coordinate fit)
+j Go to the preceding image line or column in a 2D or multispec image
+k Go to the next image line or column in a 2D or multispec image
+l Match coordinates in the coordinate (l)ist
+m (M)ark a new feature near the cursor and enter coordinate and label
+n Move the cursor or zoom to the (n)ext feature (same as +)
+o Go to the specified image line or column in a 2D or multispec image
+p (P)an to user defined window after (z)ooming on a feature
+q (Q)uit and continue with next image (also carriage return)
+r (R)edraw the graph
+t Reset the position of a feature without centering
+u Enter a new (u)ser coordinate and label for the current feature
+w (W)indow the graph. Use '?' to window prompt for more help.
+y Automatically find "maxfeatures" strongest peaks and identify them
+z (Z)oom on the feature nearest the cursor
+. Move the cursor or zoom to the feature nearest the cursor
++ Move the cursor or zoom to the next feature
+- Move the cursor or zoom to the previous feature
+I Interrupt task and exit immediately. Database information is not saved.
+
+
+COLON COMMANDS
+
+The parameters are listed or set with the following commands which may be
+abbreviated. To list the value of a parameter type the command alone.
+
+:show file Show the values of all the parameters
+:features file Write feature list to file (default is STDOUT)
+
+:coordlist file Coordinate list file
+:cradius value Centering radius in pixels
+:threshold value Detection threshold for feature centering
+:database name Database for recording feature records
+:ftype value Feature type
+ (emission|absorption|gemission|gabsorption)
+:fwidth value Feature width in pixels
+:image imagename Set a new image or show the current image
+:labels value Feature label type
+ (none|index|pixel|coords|user|both)
+:match value Coordinate list matching distance
+:maxfeatures value Maximum number of features automatically found
+:minsep value Minimum separation allowed between features
+:read name ap Read a record from the database
+ (name and ap default to the current spectrum)
+:write name ap Write a record to the database
+ (name and ap default to the current spectrum)
+:add name ap Add features from the database
+ (name and ap default to the current spectrum)
+:zwidth value Zoom width in user units
+
+Labels:
+ none - No labels
+ index - Sequential numbers in order of increasing pixel position
+ pixel - Pixel coordinates
+ coords - User coordinates such as wavelength
+ user - User labels
+ both - Combination of coords and user
diff --git a/noao/rv/rvidlines/t_identify.x b/noao/rv/rvidlines/t_identify.x
new file mode 100644
index 00000000..aac9596f
--- /dev/null
+++ b/noao/rv/rvidlines/t_identify.x
@@ -0,0 +1,108 @@
+include <mach.h>
+include <pkg/gtools.h>
+include "identify.h"
+
+# T_IDENTIFY -- Identify features and determine dispersion solutions
+
+procedure t_identify ()
+
+begin
+ call identify (IDENTIFY)
+end
+
+
+# T_RVIDLINES -- Identify features and determine radial velocities
+
+procedure t_rvidlines ()
+
+begin
+ call identify (RVIDLINES)
+end
+
+
+# IDENTIFY -- Initialize task
+
+procedure identify (taskid)
+
+int taskid #I Task ID
+
+int list, clscan(), clgeti(), clgwrd(), nscan(), imtopenp(), imtgetim()
+real clgetr()
+pointer sp, str, id, gt_init()
+
+begin
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Allocate the basic data structure.
+ call id_init (id, taskid)
+
+ # Get task parameters.
+ list = imtopenp ("images")
+ if (clscan ("nsum") != EOF) {
+ call gargi (ID_NSUM(id,1))
+ call gargi (ID_NSUM(id,2))
+ if (nscan() == 0)
+ call error (1, "Error in 'nsum' parameter")
+ if (nscan() == 1)
+ ID_NSUM(id,2) = ID_NSUM(id,1)
+ ID_NSUM(id,1) = max (1, ID_NSUM(id,1))
+ ID_NSUM(id,2) = max (1, ID_NSUM(id,2))
+ }
+ ID_MAXFEATURES(id) = clgeti ("maxfeatures")
+ ID_MINSEP(id) = clgetr ("minsep")
+ ID_MATCH(id) = clgetr ("match")
+ ID_ZWIDTH(id) = clgetr ("zwidth")
+ ID_FTYPE(id) = clgwrd ("ftype", Memc[str], SZ_LINE, FTYPES)
+ ID_FWIDTH(id) = clgetr ("fwidth")
+ ID_CRADIUS(id) = clgetr ("cradius")
+ ID_THRESHOLD(id) = clgetr ("threshold")
+ call clgstr ("database", Memc[ID_DATABASE(id)], SZ_FNAME)
+ call clgstr ("coordlist", Memc[ID_COORDLIST(id)], SZ_FNAME)
+ ID_LABELS(id) = 1
+
+ # Initialize features data structure.
+ ID_GT(id) = gt_init()
+ call gt_sets (ID_GT(id), GTTYPE, "line")
+ ID_CV(id) = NULL
+ ID_CURRENT(id) = 0
+ ID_SHIFT(id) = 0.
+ ID_REDSHIFT(id) = 0.
+
+ # Initialize ICFIT
+ call ic_open (ID_IC(id))
+ if (ID_TASK(id) == IDENTIFY) {
+ call clgstr ("function", Memc[str], SZ_LINE)
+ call ic_pstr (ID_IC(id), "function", Memc[str])
+ call ic_puti (ID_IC(id), "order", clgeti ("order"))
+ call clgstr ("sample", Memc[str], SZ_LINE)
+ call ic_pstr (ID_IC(id), "sample", Memc[str])
+ call ic_puti (ID_IC(id), "naverage", 1)
+ call ic_puti (ID_IC(id), "niterate", clgeti ("niterate"))
+ call ic_putr (ID_IC(id), "low", clgetr ("low_reject"))
+ call ic_putr (ID_IC(id), "high", clgetr ("high_reject"))
+ call ic_putr (ID_IC(id), "grow", clgetr ("grow"))
+ call ic_pstr (ID_IC(id), "xlabel", "Feature positions")
+ call ic_pstr (ID_IC(id), "xunits", "pixels")
+ call ic_pstr (ID_IC(id), "ylabel", "")
+ call ic_pkey (ID_IC(id), 1, 'y', 'x')
+ call ic_pkey (ID_IC(id), 2, 'y', 'v')
+ call ic_pkey (ID_IC(id), 3, 'y', 'r')
+ call ic_pkey (ID_IC(id), 4, 'y', 'd')
+ call ic_pkey (ID_IC(id), 5, 'y', 'n')
+ call ic_puti (ID_IC(id), "key", 3)
+ }
+
+ # Get the line list.
+ call id_mapll (id)
+
+ # Expand the image template and identify features in each image.
+ while (imtgetim (list, Memc[ID_IMAGE(id)], SZ_FNAME) != EOF)
+ call id_identify (id)
+
+ # Finish up.
+ call smw_daxis (NULL, NULL, 0, 0, 0)
+ call id_free (id)
+ call imtclose (list)
+ call sfree (sp)
+end
diff --git a/noao/rv/rvidlines/t_reidentify.x b/noao/rv/rvidlines/t_reidentify.x
new file mode 100644
index 00000000..5a724325
--- /dev/null
+++ b/noao/rv/rvidlines/t_reidentify.x
@@ -0,0 +1,1092 @@
+include <error.h>
+include <fset.h>
+include <gset.h>
+include <pkg/gtools.h>
+include <smw.h>
+include "identify.h"
+
+define ICFITHELP "noao$lib/scr/idicgfit.key"
+define VLIGHT 2.997925e5 # Speed of light, Km/sec
+
+
+# T_REIDENTIFY -- Reidentify features starting from reference features.
+
+procedure t_reidentify ()
+
+begin
+ call reiden (IDENTIFY)
+end
+
+
+# T_RVREIDLINES -- Reidentify lines for radial velocities
+
+procedure t_rvreidlines ()
+
+begin
+ call reiden (RVIDLINES)
+end
+
+
+# REIDEN -- Reidentify features starting from reference features.
+# A reference spectrum is specified and the same features are identified
+# in other images. Some lines may be lost due to bad centering. Additional
+# lines may be excluded from a new fit to the dispersion function. Instead
+# of refitting the dispersion function the user may elect to determine only
+# a shift in the reference dispersion function. Additional features may
+# be added given a coordinate list.
+#
+# In 2D images a starting line or column is selected. A number of lines
+# or columns may be averaged before identifying features. If a positive step
+# size is given then additional lines or columns may be reidentified in
+# the reference image. This may be done either by tracing or by reidentifying
+# starting from the same reference features. Reidentification between images
+# is done by taking the same line or column from the reference image.
+#
+# Multispec format images are matched by aperture number and the spectra
+# need not be in the same order in each image.
+
+procedure reiden (taskid)
+
+int taskid #I Task ID
+
+pointer reference # Reference image
+int list # List of images
+char ans[3] # Interactive?
+
+int i, fd, nlogfd
+pointer sp, logfile, str, id, logfd, pd
+
+int clscan(), clgeti(), clpopnu(), clgfil(), clgwrd(), btoi()
+int nscan(), open(), nowhite(), imtopenp(), imtgetim()
+bool clgetb()
+real clgetr()
+pointer gopen(), gt_init()
+
+begin
+ call smark (sp)
+ call salloc (reference, SZ_FNAME, TY_CHAR)
+ call salloc (logfile, SZ_FNAME, TY_CHAR)
+ call salloc (str, SZ_LINE, TY_CHAR)
+
+ # Allocate the basic data structures.
+ call id_init (id, taskid)
+ call ic_open (ID_IC(id))
+
+ # Get task parameters.
+ call clgstr ("reference", Memc[reference], SZ_FNAME)
+ list = imtopenp ("images")
+ i = nowhite (Memc[reference], Memc[reference], SZ_FNAME)
+
+ if (ID_TASK(id) == IDENTIFY)
+ ID_REFIT(id) = btoi (clgetb ("refit"))
+ else
+ ID_REFIT(id) = 3
+
+ if (clscan ("nsum") != EOF) {
+ call gargi (ID_NSUM(id,1))
+ call gargi (ID_NSUM(id,2))
+ if (nscan() == 0)
+ call error (1, "Error in 'nsum' parameter")
+ if (nscan() == 1)
+ ID_NSUM(id,2) = ID_NSUM(id,1)
+ ID_NSUM(id,1) = max (1, ID_NSUM(id,1))
+ ID_NSUM(id,2) = max (1, ID_NSUM(id,2))
+ }
+ ID_MAXFEATURES(id) = clgeti ("maxfeatures")
+ ID_MINSEP(id) = clgetr ("minsep")
+ ID_MATCH(id) = clgetr ("match")
+ if (ID_TASK(id) == IDENTIFY) {
+ ID_ZWIDTH(id) = clgetr ("identify.zwidth")
+ ID_FTYPE(id) = clgwrd ("identify.ftype", Memc[str], SZ_LINE, FTYPES)
+ ID_FWIDTH(id) = clgetr ("identify.fwidth")
+ } else {
+ ID_ZWIDTH(id) = clgetr ("rvidlines.zwidth")
+ ID_FTYPE(id) = clgwrd ("rvidlines.ftype", Memc[str], SZ_LINE,FTYPES)
+ ID_FWIDTH(id) = clgetr ("rvidlines.fwidth")
+ }
+ ID_CRADIUS(id) = clgetr ("cradius")
+ ID_THRESHOLD(id) = clgetr ("threshold")
+ call clgstr ("database", Memc[ID_DATABASE(id)], SZ_FNAME)
+ call clgstr ("coordlist", Memc[ID_COORDLIST(id)], SZ_FNAME)
+ ID_LABELS(id) = 1
+
+ call id_mapll (id)
+ ID_LOGFILES(id) = clpopnu ("logfiles")
+
+ switch (clgwrd ("interactive", ans, SZ_FNAME, "|no|yes|NO|YES|")) {
+ case 1, 3:
+ call strcpy ("NO", ans, 3)
+ ID_GP(id) = NULL
+ case 2, 4:
+ # Open graphics
+ call clgstr ("graphics", Memc[logfile], SZ_FNAME)
+ ID_GP(id) = gopen (Memc[logfile], NEW_FILE+AW_DEFER, STDGRAPH)
+ if (ID_TASK(id) == IDENTIFY) {
+ call ic_pstr (ID_IC(id), "help", ICFITHELP)
+ call ic_pstr (ID_IC(id), "xlabel", "Feature positions")
+ call ic_pstr (ID_IC(id), "xunits", "pixels")
+ call ic_pstr (ID_IC(id), "ylabel", "")
+ call ic_pkey (ID_IC(id), 1, 'y', 'x')
+ call ic_pkey (ID_IC(id), 2, 'y', 'v')
+ call ic_pkey (ID_IC(id), 3, 'y', 'r')
+ call ic_pkey (ID_IC(id), 4, 'y', 'd')
+ call ic_pkey (ID_IC(id), 5, 'y', 'n')
+ call ic_puti (ID_IC(id), "key", 3)
+ }
+ }
+
+ # Open log and plot files.
+ nlogfd = 0
+ if (clgetb ("verbose")) {
+ nlogfd = 1
+ call malloc (logfd, nlogfd, TY_INT)
+ Memi[logfd] = STDOUT
+ }
+ while (clgfil (ID_LOGFILES(id), Memc[logfile], SZ_FNAME) != EOF) {
+ fd = open (Memc[logfile], APPEND, TEXT_FILE)
+ call fseti (fd, F_FLUSHNL, YES)
+ nlogfd = nlogfd + 1
+ if (nlogfd == 1)
+ call malloc (logfd, nlogfd, TY_INT)
+ else
+ call realloc (logfd, nlogfd, TY_INT)
+ Memi[logfd+nlogfd-1] = fd
+ }
+ call ri_loghdr (id, Memc[reference], Memi[logfd], nlogfd, 1)
+
+ pd = NULL
+ if (ID_TASK(id) == IDENTIFY) {
+ call clgstr ("plotfile", Memc[logfile], SZ_FNAME)
+ if (nowhite (Memc[logfile], Memc[logfile], SZ_FNAME) > 0) {
+ fd = open (Memc[logfile], APPEND, BINARY_FILE)
+ pd = gopen ("stdvdm", NEW_FILE, fd)
+ }
+ }
+
+ ID_GT(id) = gt_init()
+ call gt_sets (ID_GT(id), GTTYPE, "line")
+
+ # Get and trace the reference solutions.
+ call ri_reference (id, Memc[reference], ans, Memi[logfd], nlogfd, pd)
+
+ # Expand the image template and reidentify features.
+ while (imtgetim (list, Memc[ID_IMAGE(id)], SZ_FNAME) != EOF)
+ call ri_image (id, Memc[reference], Memc[ID_IMAGE(id)], ans,
+ Memi[logfd], nlogfd, pd)
+
+ # Finish up.
+ if (nlogfd > 0) {
+ do i = 1, nlogfd
+ call close (Memi[logfd+i-1])
+ call mfree (logfd, TY_INT)
+ }
+
+ if (ID_GP(id) != NULL)
+ call gclose (ID_GP(id))
+ if (pd != NULL) {
+ call gclose (pd)
+ call close (fd)
+ }
+ call clpcls (ID_LOGFILES(id))
+ call imtclose (list)
+ call id_free (id)
+ call smw_daxis (NULL, NULL, 0, 0, 0)
+ call sfree (sp)
+end
+
+
+# RI_REFERENCE -- Set reference features. Trace if needed.
+
+procedure ri_reference (id, reference, ans, logfd, nlogfd, pd)
+
+pointer id # ID pointer
+char reference[ARB] # Reference image
+char ans[3] # Interactive?
+int logfd[ARB] # Logfiles
+int nlogfd # Number of logfiles
+pointer pd # Plot file pointer
+
+int step[2]
+double shift[2]
+int nreid
+bool override
+bool trace
+
+int i, start[2], line[2], loghdr
+double fit_shift[2]
+pointer ic, ic1
+bool clgetb()
+int clscan(), clgeti(), nscan(), id_gid(), id_dbcheck()
+errchk id_dbread
+
+begin
+ # Open the image and return if there is an error.
+ call strcpy (reference, Memc[ID_IMAGE(id)], SZ_FNAME)
+ iferr (call id_map (id)) {
+ call erract (EA_WARN)
+ return
+ }
+
+ # Set parameters
+ start[1] = ID_LINE(id,1)
+ start[2] = ID_LINE(id,2)
+ if (clscan ("step") != EOF) {
+ call gargi (step[1])
+ call gargi (step[2])
+ if (nscan() == 0)
+ call error (1, "Error in 'step' parameter")
+ if (nscan() == 1)
+ step[2] = step[1]
+ }
+ if (clscan ("shift") != EOF) {
+ call gargd (shift[1])
+ call gargd (shift[2])
+ if (nscan() == 0)
+ call error (1, "Error in 'shift' parameter")
+ if (nscan() == 1)
+ shift[2] = shift[1]
+ }
+ nreid = max (1, ID_NFEATURES(id) - clgeti ("nlost"))
+ override = clgetb ("override")
+ trace = clgetb ("trace")
+
+ # Get and save the reference database entry.
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO, NO)
+ call id_saveid (id, ID_LINE(id,1))
+
+ # Get and save other entries.
+ if (!override) {
+ for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) {
+ ID_LINE(id,2) = line[2]
+ ID_AP(id,2) = line[2]
+ for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+ if (line[1]==start[1] && line[2]==start[2])
+ next
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ ifnoerr (
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ NO, NO)) {
+ call id_saveid (id, ID_LINE(id,1))
+ }
+ }
+ for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+ line[1]=line[1]+step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)],
+ ID_AP(id,1), NO, NO)) {
+ call id_saveid (id, ID_LINE(id,1))
+ }
+ }
+ }
+ for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2);
+ line[2]=line[2]+step[2]) {
+ ID_LINE(id,2) = line[2]
+ ID_AP(id,2) = line[2]
+ for (line[1]=start[1]-step[1]; line[1]>0;
+ line[1]=line[1]-step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ ifnoerr (
+ call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1),
+ NO, NO)) {
+ call id_saveid (id, ID_LINE(id,1))
+ }
+ }
+ for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+ line[1]=line[1]+step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)],
+ ID_AP(id,1), NO, NO)) {
+ call id_saveid (id, ID_LINE(id,1))
+ }
+ }
+ }
+ }
+
+ # Reidentify.
+ loghdr = 2
+ ic = ID_IC(id)
+ if (ans[1] == 'N')
+ ic1 = ic
+ else {
+ call ic_open (ic1)
+ call ic_copy (ic, ic1)
+ }
+
+ fit_shift[2] = shift[2]
+ for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) {
+ ID_LINE(id,2) = line[2]
+ ID_AP(id,2) = line[2]
+ ID_IC(id) = ic
+
+ if (IS_INDEFD(shift[2]))
+ fit_shift[2] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[2] = fit_shift[2] - shift[2]
+ else
+ fit_shift[2] = -shift[2]
+ }
+
+ fit_shift[1] = fit_shift[2]
+ for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+ if (line[1]==start[1] && line[2]==start[2])
+ next
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ ID_IC(id) = ic
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ if (!trace) {
+ i = id_gid (id, start)
+ ID_LINE(id,1) = line[1]
+ ID_LINE(id,2) = line[2]
+ }
+
+ if (IS_INDEFD(shift[1]))
+ fit_shift[1] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[1] = fit_shift[1] - shift[1]
+ else
+ fit_shift[1] = -shift[1]
+ }
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) < nreid && trace) {
+ call ri_loghdr (id, reference, logfd, nlogfd, 3)
+ break
+ }
+
+ if (ID_NFEATURES(id) > 0) {
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ call id_saveid (id, line)
+ }
+ }
+
+ ID_IC(id) = ic
+ i = id_gid (id, start)
+ fit_shift[1] = fit_shift[2]
+ for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+ line[1]=line[1]+step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ ID_IC(id) = ic
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ if (!trace) {
+ i = id_gid (id, start)
+ ID_LINE(id,1) = line[1]
+ ID_LINE(id,2) = line[2]
+ }
+
+ if (IS_INDEFD(shift[1]))
+ fit_shift[1] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[1] = fit_shift[1] + shift[1]
+ else
+ fit_shift[1] = shift[1]
+ }
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) < nreid && trace) {
+ call ri_loghdr (id, reference, logfd, nlogfd, 3)
+ break
+ }
+
+ if (ID_NFEATURES(id) > 0) {
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ call id_saveid (id, line)
+ }
+ }
+ }
+
+
+ fit_shift[2] = 0.
+ for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2);
+ line[2]=line[2]+step[2]) {
+ ID_LINE(id,2) = line[2]
+ ID_AP(id,2) = line[2]
+ ID_IC(id) = ic
+
+ if (IS_INDEFD(shift[2]))
+ fit_shift[2] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[2] = fit_shift[2] + shift[2]
+ else
+ fit_shift[2] = shift[2]
+ }
+
+ fit_shift[1] = fit_shift[2]
+ for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ ID_IC(id) = ic
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ if (!trace) {
+ i = id_gid (id, start)
+ ID_LINE(id,1) = line[1]
+ ID_LINE(id,2) = line[2]
+ }
+
+ if (IS_INDEFD(shift[1]))
+ fit_shift[1] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[1] = fit_shift[1] - shift[1]
+ else
+ fit_shift[1] = -shift[1]
+ }
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) < nreid && trace) {
+ call ri_loghdr (id, reference, logfd, nlogfd, 3)
+ break
+ }
+
+ if (ID_NFEATURES(id) > 0) {
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ call id_saveid (id, line)
+ }
+ }
+
+ ID_IC(id) = ic
+ i = id_gid (id, start)
+ fit_shift[1] = fit_shift[2]
+ for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1);
+ line[1]=line[1]+step[1]) {
+ ID_LINE(id,1) = line[1]
+ ID_AP(id,1) = line[1]
+ ID_IC(id) = ic
+ if (ID_APS(id) != NULL)
+ ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1]
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ if (!trace) {
+ i = id_gid (id, start)
+ ID_LINE(id,1) = line[1]
+ ID_LINE(id,2) = line[2]
+ }
+
+ if (IS_INDEFD(shift[1]))
+ fit_shift[1] = INDEFD
+ else {
+ if (!trace)
+ fit_shift[1] = fit_shift[1] + shift[1]
+ else
+ fit_shift[1] = shift[1]
+ }
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) < nreid && trace) {
+ call ri_loghdr (id, reference, logfd, nlogfd, 3)
+ break
+ }
+
+ if (ID_NFEATURES(id) > 0) {
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ call id_saveid (id, line)
+ }
+ }
+ }
+
+ ID_IC(id) = ic
+ if (ic != ic1)
+ call ic_closed (ic1)
+
+ call smw_close (MW(ID_SH(id)))
+ call imunmap (IM(ID_SH(id)))
+ call shdr_close (ID_SH(id))
+end
+
+
+# RI_IMAGE -- Reidentify an image.
+
+procedure ri_image (id, reference, image, ans, logfd, nlogfd, pd)
+
+pointer id # ID pointer
+char reference[ARB] # Reference image
+char image[ARB] # Image to be reidentified
+char ans[3] # Interactive?
+int logfd[ARB] # Logfiles
+int nlogfd # Number of logfiles
+pointer pd # Plot file pointer
+
+bool newaps # Add new apertures not in reference?
+bool override # Override previous identifications?
+bool verbose # Verbose output?
+
+int i, j, ap, loghdr, id_getid(), id_dbcheck()
+double shift, fit_shift, clgetd()
+pointer ic, ic1
+bool clgetb()
+
+begin
+ # Open the image and return if there is an error.
+ call strcpy (image, Memc[ID_IMAGE(id)], SZ_FNAME)
+ iferr (call id_map (id)) {
+ call erract (EA_WARN)
+ return
+ }
+ call dtunmap (ID_DT(id))
+
+ newaps = clgetb ("newaps")
+ override = clgetb ("override")
+ verbose = clgetb ("verbose")
+
+ ic = ID_IC(id)
+ if (ans[1] == 'N')
+ ic1 = ic
+ else
+ call ic_open (ic1)
+
+ loghdr = 2
+ shift = clgetd ("shift")
+
+ # For MULTISPEC search the reference list of each aperture. If
+ # a reference of the same aperture is not found and the newaps
+ # flag is set use the initial reference and then add the
+ # reidentification to the reference list.
+ # For NDSPEC apply each reference to the image.
+
+ if (SMW_FORMAT(MW(ID_SH(id))) == SMW_ES ||
+ SMW_FORMAT(MW(ID_SH(id))) == SMW_MS) {
+ for (i=1; i<=ID_MAXLINE(id,1); i=i+1) {
+ ap = Memi[ID_APS(id)+i-1]
+ for (j=ID_NID(id); id_getid (id,j)!=EOF; j=j-1)
+ if (ap == ID_AP(id,1))
+ break
+
+ if (j == 0 && !newaps) {
+ if (verbose) {
+ call printf (
+ "%s: Reference for aperture %d not found\n")
+ call pargstr (image)
+ call pargi (ap)
+ }
+ next
+ }
+
+ ID_AP(id,1) = ap
+ ID_LINE(id,1) = i
+
+ if (i == 1 && ic != ic1)
+ call ic_copy (ic, ic1)
+
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+
+ fit_shift = shift
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) > 0) {
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ if (j == 0 && newaps) {
+ call id_sid (id, ID_NID(id)+1)
+ if (verbose) {
+ call printf (
+ "%s: New reference for aperture %d\n")
+ call pargstr (image)
+ call pargi (ap)
+ }
+ }
+ }
+ ID_IC(id) = ic
+ }
+
+ } else {
+ for (i=1; id_getid (id,i)!=EOF; i=i+1) {
+ if (i == 1 && ic != ic1)
+ call ic_copy (ic, ic1)
+
+ if (!override)
+ if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES)
+ next
+
+ ID_IC(id) = ic1
+ call id_gdata (id)
+ iferr (call id_fitdata (id))
+ ;
+
+ call ri_loghdr (id, reference, logfd, nlogfd, loghdr)
+ loghdr = 0
+
+ fit_shift = shift
+ call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+ if (ID_NFEATURES(id) > 0)
+ call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO)
+ ID_IC(id) = ic
+ }
+ }
+
+ ID_IC(id) = ic
+ if (ic != ic1)
+ call ic_closed (ic1)
+ call smw_close (MW(ID_SH(id)))
+ call imunmap (IM(ID_SH(id)))
+ call shdr_close (ID_SH(id))
+end
+
+
+# RI_REIDENTIFY -- Reidentify features using a reference image database entry.
+
+procedure ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd)
+
+pointer id # ID pointer
+double fit_shift # Shift in fit coords (input and output)
+char ans[3] # Interactive?
+int logfd[ARB] # Logfiles
+int nlogfd # Number of logfiles
+pointer pd # Plot file pointer
+
+int i, j, nfeatures1, nfeatures2, nfit, iden, mono, clgwrd()
+double shift, pix_shift, z_shift, v, vrms
+double id_fitpt(), fit_to_pix(), id_shift(), id_center(), id_rms(), id_zval()
+pointer sp, pix, fit
+bool clgetb()
+
+begin
+ nfeatures1 = ID_NFEATURES(id)
+ call smark (sp)
+ call salloc (pix, nfeatures1, TY_DOUBLE)
+ call salloc (fit, nfeatures1, TY_DOUBLE)
+ call amovd (PIX(id,1), Memd[pix], nfeatures1)
+ call amovd (FIT(id,1), Memd[fit], nfeatures1)
+
+ # If no initial shift is given then the procedure id_shift
+ # computes a shift between the reference features and the
+ # features in the image. The purpose of the shift is to get the
+ # reference feature positions close enough to those of the image
+ # being identified that the centering algorithm will determine
+ # the exact positions of the features. An initial shift of zero
+ # is used if the two images are very nearly aligned as in the
+ # case of tracing features in a two dimensional image or for a
+ # set of images taken with the same observing setup.
+
+ if (IS_INDEFD(fit_shift)) {
+ ID_FWIDTH(id) = FWIDTH(id,1)
+ ID_FTYPE(id) = FTYPE(id,1)
+ shift = id_shift (id)
+ } else
+ shift = fit_shift
+
+ # For each reference feature a shift is added to bring the pixel
+ # position near that for the image being identified and then the
+ # centering algorithm is used. If the centering algorithm fails
+ # the feature is discarded. A mean shift is computed for the
+ # features which have been reidentified.
+
+ do i = 1, ID_NFEATURES(id) {
+ PIX(id,i) = fit_to_pix (id, FIT(id,i) + shift)
+ PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i), FTYPE(id,i),
+ NO)
+ if (!IS_INDEFD(PIX(id,i)))
+ FIT(id,i) = id_fitpt (id, PIX(id,i))
+ }
+ for (i=1; i<ID_NFEATURES(id); i=i+1) {
+ if (IS_INDEFD(PIX(id,i)))
+ next
+ for (j=i+1; j<=ID_NFEATURES(id); j=j+1) {
+ if (IS_INDEFD(PIX(id,j)))
+ next
+ if (abs (PIX(id,i)-PIX(id,j)) < ID_MINSEP(id)) {
+ if (abs (FIT(id,i)-USER(id,i)) < abs (FIT(id,j)-USER(id,j)))
+ PIX(id,j) = INDEFD
+ else {
+ PIX(id,i) = INDEFD
+ break
+ }
+ }
+ }
+ }
+
+ pix_shift = 0.
+ fit_shift = 0.
+ z_shift = 0.
+ j = 0
+ do i = 1, ID_NFEATURES(id) {
+ if (IS_INDEFD(PIX(id,i)))
+ next
+
+ pix_shift = pix_shift + PIX(id,i) - Memd[pix+i-1]
+ fit_shift = fit_shift + FIT(id,i) - Memd[fit+i-1]
+ if (Memd[fit+i-1] != 0.)
+ z_shift = z_shift + id_zval (id, FIT(id,i), Memd[fit+i-1])
+
+ j = j + 1
+ PIX(id,j) = PIX(id,i)
+ FIT(id,j) = FIT(id,i)
+ USER(id,j) = USER(id,i)
+ WTS(id,j) = WTS(id,i)
+ FWIDTH(id,j) = FWIDTH(id,i)
+ FTYPE(id,j) = FTYPE(id,i)
+ }
+ ID_NFEATURES(id) = j
+
+ nfeatures2 = j
+ pix_shift = pix_shift / max (1, ID_NFEATURES(id))
+ fit_shift = fit_shift / max (1, ID_NFEATURES(id))
+ z_shift = z_shift / max (1, ID_NFEATURES(id))
+
+ # If refitting the coordinate function is requested and there is
+ # more than one feature and there is a previously defined
+ # coordinate function then refit. Otherwise compute a coordinate
+ # shift.
+
+ mono = YES
+ switch (ID_REFIT(id)) {
+ case 1:
+ if (ID_CV(id) != NULL) {
+ if (ID_NFEATURES(id)>1)
+ call id_dofit (id, NO)
+ else
+ call id_doshift (id, NO)
+ }
+ case 2:
+ if (ID_CV(id) != NULL)
+ call id_doshift (id, NO)
+ case 3:
+ call id_velocity (id, NO)
+ v = ID_REDSHIFT(id) * VLIGHT
+ vrms = ID_RMSRED(id) * VLIGHT
+ }
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata (id))
+ mono = NO
+ call id_fitfeatures (id)
+ }
+
+ if (clgetb ("addfeatures")) {
+ ID_FWIDTH(id) = FWIDTH(id,1)
+ ID_FTYPE(id) = FTYPE(id,1)
+ call id_linelist (id)
+ if (ID_NEWFEATURES(id) == YES) {
+ switch (ID_REFIT(id)) {
+ case 1:
+ if (ID_NFEATURES(id)>1)
+ call id_dofit (id, NO)
+ else
+ call id_doshift (id, NO)
+ case 2:
+ call id_doshift (id, NO)
+ case 3:
+ call id_velocity (id, NO)
+ v = ID_REDSHIFT(id) * VLIGHT
+ vrms = ID_RMSRED(id) * VLIGHT
+ }
+ if (ID_NEWCV(id) == YES) {
+ iferr (call id_fitdata (id))
+ mono = NO
+ call id_fitfeatures (id)
+ }
+ }
+ }
+
+ # Enter fitting interactively.
+ iden = NO
+ if ((ID_NFEATURES(id)>1) && (ID_CV(id)!=NULL || ID_REFIT(id) == 3)) {
+ if (ans[1] != 'N') {
+ if (ans[1] != 'Y') {
+ nfit = 0
+ for (j=1; j<=ID_NFEATURES(id); j=j+1)
+ if (WTS(id,j) > 0.)
+ nfit = nfit + 1
+ call printf (
+ "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargi (nfeatures2)
+ call pargi (nfeatures1)
+ call pargi (nfit)
+ call pargi (ID_NFEATURES(id))
+ call pargd (pix_shift)
+ call pargd (fit_shift)
+ if (ID_REFIT(id) == 3) {
+ call pargd (v)
+ call pargd (vrms)
+ } else {
+ call pargd (z_shift)
+ call pargd (id_rms(id))
+ }
+ call flush (STDOUT)
+ repeat {
+ ifnoerr (i = clgwrd ("answer", ans, SZ_FNAME,
+ "|no|yes|NO|YES|"))
+ break
+ }
+ call clpstr ("answer", ans)
+ }
+ switch (ans[1]) {
+ case 'y', 'Y':
+ mono = YES
+ i = ID_REFIT(id)
+ call reidentify (id)
+ ID_REFIT(id) = i
+ if (ID_REFIT(id) == 3) {
+ call id_velocity (id, NO)
+ v = ID_REDSHIFT(id) * VLIGHT
+ vrms = ID_RMSRED(id) * VLIGHT
+ }
+ iden = YES
+ }
+ if (ans[1] != 'Y')
+ call gdeactivate (ID_GP(id), 0)
+ }
+ }
+
+ # Record log information if a log file descriptor is given.
+ for (i = 1; i <= nlogfd; i = i + 1) {
+ if (ans[1] == 'n' && logfd[i] == STDOUT)
+ next
+ nfit = 0
+ for (j=1; j<=ID_NFEATURES(id); j=j+1)
+ if (WTS(id,j) > 0.)
+ nfit = nfit + 1
+ call fprintf (logfd[i],
+ "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargi (nfeatures2)
+ call pargi (nfeatures1)
+ call pargi (nfit)
+ call pargi (ID_NFEATURES(id))
+ call pargd (pix_shift)
+ call pargd (fit_shift)
+ if (ID_REFIT(id) == 3) {
+ call pargd (v)
+ call pargd (vrms)
+ } else {
+ call pargd (z_shift)
+ call pargd (id_rms(id))
+ }
+ if (ID_REFIT(id) == 3 && logfd[i] != STDOUT)
+ call id_log (id, "", logfd[i])
+ if (mono == NO)
+ call fprintf (logfd[i], "Non-monotonic dispersion function")
+ call flush (logfd[i])
+ if (logfd[i] == STDOUT)
+ iden = NO
+ }
+ # Print log if STDOUT is not used but if the IDENTIFY is done.
+ if (iden == YES) {
+ call printf (
+ "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n")
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargstr (Memc[ID_SECTION(id)])
+ call pargi (nfeatures2)
+ call pargi (nfeatures1)
+ call pargi (nfit)
+ call pargi (ID_NFEATURES(id))
+ call pargd (pix_shift)
+ call pargd (fit_shift)
+ if (ID_REFIT(id) == 3) {
+ call pargd (v)
+ call pargd (vrms)
+ } else {
+ call pargd (z_shift)
+ call pargd (id_rms(id))
+ }
+ if (mono == NO)
+ call printf ("Non-monotonic dispersion function")
+ call flush (STDOUT)
+ }
+
+ # Make log plot.
+ call ri_plot (id, pd)
+
+ call sfree (sp)
+end
+
+
+# RI_LOGHDR -- Print a log header in the log files.
+
+procedure ri_loghdr (id, reference, logfd, nlogfd, flag)
+
+pointer id # Identify structure
+char reference[ARB] # Reference image
+int logfd[ARB] # Log file descriptors
+int nlogfd # Number of log files
+int flag # Header type flag (1=banner, 2=Column labels, 3=Error)
+
+int i
+pointer str
+
+begin
+ for (i = 1; i <= nlogfd; i = i + 1) {
+ switch (flag) {
+ case 1: # Print ID
+ call malloc (str, SZ_LINE, TY_CHAR)
+ call sysid (Memc[str], SZ_LINE)
+ switch (ID_REFIT(id)) {
+ case 1, 2:
+ call fprintf (logfd[i], "\nREIDENTIFY: %s\n")
+ call pargstr (Memc[str])
+ case 3:
+ call fprintf (logfd[i], "\nRVREIDLINES: %s\n")
+ call pargstr (Memc[str])
+ }
+ call mfree (str, TY_CHAR)
+ case 2: # Print labels
+ switch (ID_REFIT(id)) {
+ case 1, 2:
+ call fprintf (logfd[i],
+ " Reference image = %s, New image = %s, Refit = %b\n")
+ call pargstr (reference)
+ call pargstr (Memc[ID_IMAGE(id)])
+ call pargi (ID_REFIT(id))
+ call fprintf (logfd[i],
+ "%20s %7s %7s %9s %10s %7s %7s\n")
+ call pargstr ("Image Data")
+ call pargstr ("Found")
+ call pargstr ("Fit")
+ call pargstr ("Pix Shift")
+ call pargstr ("User Shift")
+ call pargstr ("Z Shift")
+ call pargstr ("RMS")
+ case 3:
+ call fprintf (logfd[i],
+ " Reference image = %s, New image = %s\n")
+ call pargstr (reference)
+ call pargstr (Memc[ID_IMAGE(id)])
+ call fprintf (logfd[i],
+ "%20s %7s %7s %9s %10s %8s %7s\n")
+ call pargstr ("Image Data")
+ call pargstr ("Found")
+ call pargstr ("Fit")
+ call pargstr ("Pix Shift")
+ call pargstr ("User Shift")
+ call pargstr ("Velocity")
+ call pargstr ("RMS")
+ }
+ case 3: # Error
+ call fprintf (logfd[i], " ** Too many features lost **\n")
+ }
+ }
+end
+
+
+# RI_PLOT -- Plot residual graph of reidentified lines.
+
+procedure ri_plot (id, pd)
+
+pointer id # ID pointer
+pointer pd # GIO pointer
+
+int i, j
+pointer sp, str, x, y, gt, gt_init()
+double id_zshiftd()
+
+begin
+ # Check if there is anything to plot.
+ if (pd == NULL || ID_NFEATURES(id) == 0)
+ return
+
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (x, ID_NFEATURES(id), TY_REAL)
+ call salloc (y, ID_NFEATURES(id), TY_REAL)
+
+ # Set plot points.
+ j = 0
+ do i = 1, ID_NFEATURES(id) {
+ if (IS_INDEFD(USER(id,i)))
+ break
+
+ Memr[x+j] = USER(id,i)
+ Memr[y+j] = id_zshiftd (id, FIT(id,i), 0)
+ j = j + 1
+ }
+
+ if (j == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Make the plot.
+ call sprintf (Memc[str], SZ_LINE, "Reidentify: %s")
+ call pargstr (Memc[ID_IMAGE(id)])
+ gt = gt_init ()
+ call gt_sets (gt, GTTYPE, "mark")
+ call gt_sets (gt, GTXLABEL, "user coordinates")
+ call gt_sets (gt, GTYLABEL, "residuals (fit - user)")
+ call gt_sets (gt, GTTITLE, Memc[str])
+ call gclear (pd)
+ call gascale (pd, Memr[x], j, 1)
+ call gascale (pd, Memr[y], j, 2)
+ call gt_swind (pd, gt)
+ call gt_labax (pd, gt)
+ call gt_plot (pd, gt, Memr[x], Memr[y], j)
+ call gt_free (gt)
+
+ call sfree (sp)
+end