diff options
Diffstat (limited to 'noao/onedspec/ecidentify')
47 files changed, 4241 insertions, 0 deletions
diff --git a/noao/onedspec/ecidentify/eccenter.x b/noao/onedspec/ecidentify/eccenter.x new file mode 100644 index 00000000..730ad2a8 --- /dev/null +++ b/noao/onedspec/ecidentify/eccenter.x @@ -0,0 +1,34 @@ +include "ecidentify.h" + +# EC_CENTER -- Locate the center of a feature. + +double procedure ec_center (ec, x, width, type) + +pointer ec # EC pointer +double x # Initial guess +real width # Feature width +int type # Feature type + +double dvalue +real value + +real center1d() +double smw_c1trand() + +begin + if (IS_INDEFD(x)) + return (x) + + dvalue = smw_c1trand (EC_PL(ec), x) + if (IS_INDEFD(dvalue)) + return (dvalue) + + value = dvalue + value = center1d (value, IMDATA(ec,1), EC_NPTS(ec), width, + abs (type), EC_CRADIUS(ec), EC_THRESHOLD(ec)) + + if (IS_INDEF(value)) + return (INDEFD) + else + return (smw_c1trand (EC_LP(ec), double(value))) +end diff --git a/noao/onedspec/ecidentify/eccolon.x b/noao/onedspec/ecidentify/eccolon.x new file mode 100644 index 00000000..0fe22af5 --- /dev/null +++ b/noao/onedspec/ecidentify/eccolon.x @@ -0,0 +1,243 @@ +include <gset.h> +include <error.h> +include <pkg/center1d.h> +include "ecidentify.h" + +# List of colon commands. +define CMDS "|show|features|image|database|read|write|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 DATABASE 4 # Set new database +define READ 5 # Read database entry +define WRITE 6 # Write database entry +define COORDLIST 7 # Set new coordinate list +define MATCH 8 # Set coordinate list matching distance +# newline 9 +define MAXFEATURES 10 # Set maximum number of features for auto find +define MINSEP 11 # Set minimum separation distance +define ZWIDTH 12 # Set zoom window width +define LABEL 13 # Set label type +define WIDTH 14 # Set centering width +define TYPE 15 # Set centering type +define RADIUS 16 # Set centering radius +define THRESHOLD 17 # Set the centering threshold + +# EC_COLON -- Respond to colon command. + +procedure ec_colon (ec, cmdstr, newimage, prfeature) + +pointer ec # 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 +real rval[2] +pointer im + +int nscan(), strdic(), ec_next() +pointer immap() +errchk immap, ec_dbread, ec_dbwrite, ec_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 (EC_GP(ec), AW_CLEAR) + call ec_show (ec, "STDOUT") + call greactivate (EC_GP(ec), AW_PAUSE) + } else { + iferr (call ec_show (ec, cmd)) { + call erract (EA_WARN) + prfeature = NO + } + } + case FEATURES: # :features - list features + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call gdeactivate (EC_GP(ec), AW_CLEAR) + call ec_log (ec, "STDOUT") + call greactivate (EC_GP(ec), AW_PAUSE) + } else { + iferr (call ec_log (ec, cmd)) { + call erract (EA_WARN) + prfeature = NO + } + } + case IMAGE: # :image - set image + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("image %s\n") + call pargstr (Memc[EC_IMAGE(ec)]) + 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 DATABASE: # :database - set database + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("database %s\n") + call pargstr (Memc[EC_DATABASE(ec)]) + prfeature = NO + } else { + call strcpy (cmd, Memc[EC_DATABASE(ec)], SZ_FNAME) + EC_NEWDBENTRY(ec) = YES + } + case READ: # :read - read database entry + prfeature = NO + iferr { + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) + call ec_dbread (ec, Memc[EC_IMAGE(ec)], YES) + else { + call xt_stripwhite (cmd) + if (cmd[1] == EOS) + call ec_dbread (ec, Memc[EC_IMAGE(ec)], YES) + else + call ec_dbread (ec, cmd, YES) + } + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + } then + call erract (EA_WARN) + case WRITE: # :write - write database entry + prfeature = NO + iferr { + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) + call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], YES) + else { + call xt_stripwhite (cmd) + if (cmd[1] == EOS) + call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], YES) + else + call ec_dbwrite (ec, cmd, 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[EC_COORDLIST(ec)]) + prfeature = NO + } else { + call strcpy (cmd, Memc[EC_COORDLIST(ec)], SZ_FNAME) + call ec_unmapll (ec) + call ec_mapll (ec) + } + case MATCH: # :match - set matching distance for coordinate list + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("match %g\n") + call pargr (EC_MATCH(ec)) + prfeature = NO + } else + EC_MATCH(ec) = rval[1] + case MAXFEATURES: # :maxfeatures - set max num features for auto find + call gargi (ival) + if (nscan() == 1) { + call printf ("maxfeatures %d\n") + call pargi (EC_MAXFEATURES(ec)) + prfeature = NO + } else + EC_MAXFEATURES(ec) = ival + case MINSEP: # :minsep - set minimum feature separation allowed + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("minsep %g\n") + call pargr (EC_MINSEP(ec)) + prfeature = NO + } else + EC_MINSEP(ec) = rval[1] + case ZWIDTH: # :zwidth - set zoom window width + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("zwidth %g\n") + call pargr (EC_ZWIDTH(ec)) + prfeature = NO + } else { + EC_ZWIDTH(ec) = rval[1] + if (EC_GTYPE(ec) == 2) + EC_NEWGRAPH(ec) = YES + } + case LABEL: # :labels - set label type + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + switch (EC_LABELS(ec)) { + case 2: + call printf ("labels index\n") + case 3: + call printf ("labels pixel\n") + case 4: + call printf ("labels user\n") + default: + call printf ("labels none\n") + } + prfeature = NO + } else { + EC_LABELS(ec) = strdic (cmd, cmd, SZ_LINE, LABELS) + do i = 1, EC_NFEATURES(ec) { + if (APN(ec,i) == EC_AP(ec)) + call ec_mark (ec, i) + } + } + case WIDTH: # :fwidth - set centering width + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("fwidth %g\n") + call pargr (EC_FWIDTH(ec)) + prfeature = NO + } else + EC_FWIDTH(ec) = rval[1] + case TYPE: # :ftype - set centering type + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + switch (EC_FTYPE(ec)) { + case EMISSION: + call printf ("ftype emission\n") + case ABSORPTION: + call printf ("ftype absorption\n") + } + prfeature = NO + } else + EC_FTYPE(ec) = 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 (EC_CRADIUS(ec)) + prfeature = NO + } else + EC_CRADIUS(ec) = rval[1] + case THRESHOLD: # :threshold - set centering threshold + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("threshold %g\n") + call pargr (EC_THRESHOLD(ec)) + prfeature = NO + } else + EC_THRESHOLD(ec) = rval[1] + default: + call printf ("Unrecognized or ambiguous command\007") + prfeature = NO + } +end diff --git a/noao/onedspec/ecidentify/ecdb.x b/noao/onedspec/ecidentify/ecdb.x new file mode 100644 index 00000000..f6e02526 --- /dev/null +++ b/noao/onedspec/ecidentify/ecdb.x @@ -0,0 +1,268 @@ +include <math/gsurfit.h> +include <smw.h> +include <units.h> +include "ecidentify.h" + +# EC_DBREAD -- Read features data from the database. + +procedure ec_dbread (ec, name, verbose) + +pointer ec # ID pointer +char name[SZ_LINE] +int verbose + +pointer dt +int i, j, k, ncoeffs, rec, slope, offset, niterate +double shift, low, high +pointer sp, coeffs, line, cluster, un + +int ec_line() +int dtgeti(), dgsgeti(), dtlocate(), dtscan(), nscan() +real dtgetr() +bool un_compare() +double dgsgetd(), smw_c1trand() +pointer dtmap1(), un_open() + +errchk dtmap1, dtlocate, dtgeti, dtgad, un_open + +begin + call smark (sp) + call salloc (cluster, SZ_LINE, TY_CHAR) + call salloc (line, SZ_LINE, TY_CHAR) + + call imgcluster (name, Memc[cluster], SZ_LINE) + call sprintf (Memc[line], SZ_LINE, "ec%s") + call pargstr (Memc[cluster]) + dt = dtmap1 (Memc[EC_DATABASE(ec)], Memc[line], READ_ONLY) + + call sprintf (Memc[line], SZ_LINE, "ecidentify %s") + call pargstr (Memc[cluster]) + + rec = dtlocate (dt, Memc[line]) + if (rec == EOF) + call error (0, "Entry not found") + + i = dtgeti (dt, rec, "features") + + EC_NALLOC(ec) = i + call realloc (EC_APNUM(ec), i, TY_INT) + call realloc (EC_LINENUM(ec), i, TY_INT) + call realloc (EC_ORD(ec), i, TY_INT) + call realloc (EC_PIX(ec), i, TY_DOUBLE) + call realloc (EC_FIT(ec), i, TY_DOUBLE) + call realloc (EC_USER(ec), i, TY_DOUBLE) + call realloc (EC_FWIDTHS(ec), i, TY_REAL) + call realloc (EC_FTYPES(ec), i, TY_INT) + + j = 1 + do i = 1, EC_NALLOC(ec) { + k = dtscan (dt) + call gargi (APN(ec,j)) + call gargi (ORDER(ec,j)) + call gargd (PIX(ec,j)) + call gargd (FIT(ec,j)) + call gargd (USER(ec,j)) + call gargr (FWIDTH(ec,j)) + call gargi (FTYPE(ec,j)) + call gargi (k) + if (nscan() == 8 && k == 0) + FTYPE(ec,j) = -FTYPE(ec,j) + iferr (LINE(ec,j) = ec_line (ec, APN(ec,j))) + next + shift = smw_c1trand (EC_PL(ec), PIX(ec,j)) + low = 0.5 + high = SN(SH(ec,LINE(ec,j))) + 0.5 + if (shift < low || shift > high) + next + j = j + 1 + } + EC_NFEATURES(ec) = j - 1 + + iferr (shift = dtgetr (dt, rec, "shift")) + shift = 0. + iferr (offset = dtgeti (dt, rec, "offset")) + offset = 0 + iferr (slope = dtgeti (dt, rec, "slope")) + slope = 1 + call ecf_setd ("shift", shift) + call ecf_seti ("offset", offset) + call ecf_seti ("slope", slope) + + iferr { + ncoeffs = dtgeti (dt, rec, "coefficients") + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call dtgad (dt, rec, "coefficients", Memd[coeffs], ncoeffs, ncoeffs) + + if (EC_ECF(ec) != NULL) + call dgsfree (EC_ECF(ec)) + call dgsrestore (EC_ECF(ec), Memd[coeffs]) + + call ecf_setd ("xmin", dgsgetd (EC_ECF(ec), GSXMIN)) + call ecf_setd ("xmax", dgsgetd (EC_ECF(ec), GSXMAX)) + call ecf_setd ("ymin", dgsgetd (EC_ECF(ec), GSYMIN)) + call ecf_setd ("ymax", dgsgetd (EC_ECF(ec), GSYMAX)) + call ecf_seti ("xorder", dgsgeti (EC_ECF(ec), GSXORDER)) + call ecf_seti ("yorder", dgsgeti (EC_ECF(ec), GSYORDER)) + + switch (dgsgeti (EC_ECF(ec), GSTYPE)) { + case GS_LEGENDRE: + call ecf_sets ("function", "legendre") + case GS_CHEBYSHEV: + call ecf_sets ("function", "chebyshev") + } + + ifnoerr (niterate = dtgeti (dt, rec, "niterate")) + call ecf_seti ("niterate", niterate) + ifnoerr (low = dtgetr (dt, rec, "lowreject")) + call ecf_setd ("low", low) + ifnoerr (high = dtgeti (dt, rec, "highreject")) + call ecf_setd ("high", high) + + EC_NEWECF(ec) = YES + EC_CURRENT(ec) = min (1, EC_NFEATURES(ec)) + } then + ; + + ifnoerr (call dtgstr (dt, rec, "units", Memc[line], SZ_LINE)) { + if (EC_UN(ec) == NULL) + EC_UN(ec) = un_open (Memc[line]) + else { + un = un_open (Memc[line]) + if (!un_compare (un, EC_UN(ec))) { + call ec_unitsll (ec, Memc[line]) + call un_close (EC_UN(ec)) + EC_UN(ec) = un + } else + call un_close (un) + } + } + + call dtunmap (dt) + call sfree (sp) + + if (EC_NFEATURES(ec) > 0) { + EC_NEWGRAPH(ec) = YES + EC_NEWFEATURES(ec) = YES + EC_CURRENT(ec) = 1 + } else + EC_CURRENT(ec) = 0 + + if (verbose == YES) { + call printf ("ecidentify %s\n") + call pargstr (Memc[cluster]) + } +end + + +# EC_DBWRITE -- Write features data to the database. + +procedure ec_dbwrite (ec, name, verbose) + +pointer ec # ID pointer +char name[ARB] +int verbose + +int i, ncoeffs +pointer dt, sp, coeffs, root, cluster + +int dgsgeti(), ecf_geti() +double ecf_getd() +pointer dtmap1(), immap() + +errchk dtmap1, immap + +begin + call smark (sp) + call salloc (cluster, SZ_FNAME, TY_CHAR) + call salloc (root, SZ_FNAME, TY_CHAR) + + call imgcluster (name, Memc[cluster], SZ_FNAME) + call sprintf (Memc[root], SZ_FNAME, "ec%s") + call pargstr (Memc[cluster]) + dt = dtmap1 (Memc[EC_DATABASE(ec)], Memc[root], APPEND) + + call dtptime (dt) + call dtput (dt, "begin\tecidentify %s\n") + call pargstr (Memc[cluster]) + call dtput (dt, "\tid\t%s\n") + call pargstr (Memc[cluster]) + call dtput (dt, "\ttask\tecidentify\n") + call dtput (dt, "\timage\t%s\n") + call pargstr (Memc[EC_IMAGE(ec)]) + + if (EC_UN(ec) != NULL) { + call dtput (dt, "\tunits\t%s\n") + call pargstr (UN_UNITS(EC_UN(ec))) + } + call dtput (dt, "\tfeatures\t%d\n") + call pargi (EC_NFEATURES(ec)) + do i = 1, EC_NFEATURES(ec) { + call dtput (dt, + "\t\t%3d %3d %7.2f %10.9g %10.9g %4.1f %d %d\n") + call pargi (APN(ec,i)) + call pargi (ORDER(ec,i)) + call pargd (PIX(ec,i)) + call pargd (FIT(ec,i)) + call pargd (USER(ec,i)) + call pargr (FWIDTH(ec,i)) + call pargi (abs (FTYPE(ec,i))) + if (FTYPE(ec,i) > 0) + call pargi (1) + else + call pargi (0) + } + + if (ecf_getd ("shift") != 0.) { + call dtput (dt, "\tshift\t%g\n") + call pargd (ecf_getd ("shift")) + } + if (ecf_geti ("offset") != 0) { + call dtput (dt, "\toffset\t%d\n") + call pargi (ecf_geti ("offset")) + } + if (ecf_geti ("slope") != 1) { + call dtput (dt, "\tslope\t%d\n") + call pargi (ecf_geti ("slope")) + } + + if (EC_ECF(ec) != NULL) { + call dtput (dt, "\tniterate %d\n") + call pargi (ecf_geti ("niterate")) + call dtput (dt, "\tlowreject %g\n") + call pargd (ecf_getd ("low")) + call dtput (dt, "\thighreject %g\n") + call pargd (ecf_getd ("high")) + + ncoeffs = dgsgeti (EC_ECF(ec), GSNSAVE) + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call dgssave (EC_ECF(ec), 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") + call dtunmap (dt) + + EC_NEWFEATURES(ec) = NO + EC_NEWECF(ec) = NO + EC_NEWDBENTRY(ec) = NO + + if (verbose == YES) { + call printf ("ecidentify %s\n") + call pargstr (Memc[cluster]) + } + + # Enter reference spectrum name in image header. + call imgcluster (Memc[EC_IMAGE(ec)], Memc[root], SZ_FNAME) + dt = immap (Memc[root], READ_WRITE, 0) + call imastr (dt, "REFSPEC1", Memc[cluster]) + iferr (call imdelf (dt, "REFSPEC2")) + ; + call imunmap (dt) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecdelete.x b/noao/onedspec/ecidentify/ecdelete.x new file mode 100644 index 00000000..b729d326 --- /dev/null +++ b/noao/onedspec/ecidentify/ecdelete.x @@ -0,0 +1,28 @@ +include "ecidentify.h" + +# EC_DELETE -- Delete a feature. + +procedure ec_delete (ec, feature) + +pointer ec # ID pointer +int feature # Feature to be deleted + +int i + +begin + if (feature == 0) + return + + do i = feature + 1, EC_NFEATURES(ec) { + APN(ec,i-1) = APN(ec,i) + LINE(ec,i-1) = LINE(ec,i) + ORDER(ec,i-1) = ORDER(ec,i) + PIX(ec,i-1) = PIX(ec,i) + FIT(ec,i-1) = FIT(ec,i) + USER(ec,i-1) = USER(ec,i) + FWIDTH(ec,i-1) = FWIDTH(ec,i) + FTYPE(ec,i-1) = FTYPE(ec,i) + } + EC_NFEATURES(ec) = EC_NFEATURES(ec) - 1 + EC_NEWFEATURES(ec) = YES +end diff --git a/noao/onedspec/ecidentify/ecdofit.x b/noao/onedspec/ecidentify/ecdofit.x new file mode 100644 index 00000000..14dcea54 --- /dev/null +++ b/noao/onedspec/ecidentify/ecdofit.x @@ -0,0 +1,128 @@ +include <smw.h> +include "ecidentify.h" + +# EC_DOFIT -- Fit an echelle function to the features. Eliminate INDEF points. + +procedure ec_dofit (ec, interactive, fixedorder) + +pointer ec # EC pointer +int interactive # Interactive fit? +int fixedorder # Fixed order? + +int i, j, k, nfit +double xmin, xmax, ymin, ymax +pointer gt1, ecf +pointer sp, x, y, z, w, gt_init() +errchk ecf_fit + +begin + # Count number of points and determine the order range. + j = ORDER(ec,1) + k = ORDER(ec,1) + nfit = 0 + for (i=1; i<=EC_NFEATURES(ec); i=i+1) { + if (IS_INDEFD (PIX(ec,i)) || IS_INDEFD (USER(ec,i))) + next + j = min (j, ORDER(ec,i)) + k = max (k, ORDER(ec,i)) + nfit = nfit + 1 + } + + # Require at least 4 points and more than one order. + if (nfit < 4 || j == k) { + if (EC_ECF(ec) != NULL) { + call dgs_free (EC_ECF(ec)) + call ecf_setd ("shift", 0.D0) + EC_NEWGRAPH(ec) = YES + EC_NEWECF(ec) = YES + } + return + } + + # Allocate arrays for points to be fit and fill them in. + call smark (sp) + call salloc (x, nfit, TY_DOUBLE) + call salloc (y, nfit, TY_DOUBLE) + call salloc (z, nfit, TY_DOUBLE) + call salloc (w, nfit, TY_DOUBLE) + + nfit = 0 + do i = 1, EC_NFEATURES(ec) { + if (IS_INDEFD (PIX(ec,i)) || IS_INDEFD (USER(ec,i))) + next + Memd[x+nfit] = PIX(ec,i) + Memd[y+nfit] = APN(ec,i) + Memd[z+nfit] = USER(ec,i) + Memd[w+nfit] = 1. + nfit = nfit + 1 + } + + # Initialize fit limits. + ymin = APS(ec,1) + ymax = ymin + do i = 2, EC_NLINES(ec) { + xmin = APS(ec,i) + if (xmin < ymin) + ymin = xmin + if (xmin > ymax) + ymax = xmin + } + xmin = 1 + xmax = EC_NPTS(ec) + + call ecf_setd ("xmin", xmin) + call ecf_setd ("xmax", xmax) + call ecf_setd ("ymin", ymin) + call ecf_setd ("ymax", ymax) + + # Fit the echelle dispersion function. + ecf = EC_ECF(ec) + if (interactive == YES) { + gt1 = gt_init() + call ecf_fit (ecf, EC_GP(ec), gt1, Memd[x], Memd[y], + Memd[z], Memd[w], nfit, fixedorder) + call gt_free (gt1) + } else + call ecf_fit (ecf, NULL, NULL, Memd[x], Memd[y], Memd[z], + Memd[w], nfit, fixedorder) + EC_ECF(ec) = ecf + + # Remove any deleted points. + j = 0 + k = 0 + do i = 1, EC_NFEATURES(ec) { + if (IS_INDEFD (PIX(ec,i)) || IS_INDEFD (USER(ec,i))) { + j = j + 1 + APN(ec,j) = APN(ec,i) + LINE(ec,j) = LINE(ec,i) + ORDER(ec,j) = ORDER(ec,i) + PIX(ec,j) = PIX(ec,i) + FIT(ec,j) = FIT(ec,i) + USER(ec,j) = USER(ec,i) + FWIDTH(ec,j) = FWIDTH(ec,i) + FTYPE(ec,j) = abs (FTYPE(ec,i)) + } else { + if (Memd[w+k] != 0.) { + j = j + 1 + APN(ec,j) = APN(ec,i) + LINE(ec,j) = LINE(ec,i) + ORDER(ec,j) = ORDER(ec,i) + PIX(ec,j) = PIX(ec,i) + FIT(ec,j) = FIT(ec,i) + USER(ec,j) = USER(ec,i) + FWIDTH(ec,j) = FWIDTH(ec,i) + FTYPE(ec,j) = abs (FTYPE(ec,i)) + if (Memd[w+k] < 0.) + FTYPE(ec,j) = -FTYPE(ec,j) + } + k = k + 1 + } + } + EC_NFEATURES(ec) = j + + # Set flags. + EC_NEWECF(ec) = YES + EC_NEWGRAPH(ec) = YES + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecdoshift.x b/noao/onedspec/ecidentify/ecdoshift.x new file mode 100644 index 00000000..1689bc92 --- /dev/null +++ b/noao/onedspec/ecidentify/ecdoshift.x @@ -0,0 +1,44 @@ +include "ecidentify.h" + +# EC_DOSHIFT -- Minimize residuals by constant shift. + +procedure ec_doshift (ec, interactive) + +pointer ec # ID pointer +int interactive # Called interactively? + +int i, j +double shft, delta, rms, ec_fitpt(), ecf_getd() + +begin + shft = 0. + rms = 0. + j = 0 + for (i=1; i <= EC_NFEATURES(ec); i = i + 1) { + if (IS_INDEFD (USER(ec,i))) + next + delta = USER(ec,i) - ec_fitpt (ec, APN(ec,i), PIX(ec,i)) + delta = delta * ORDER(ec,i) + shft = shft + delta + rms = rms + delta * delta + j = j + 1 + } + + if (j > 0) { + shft = shft / j + rms = rms / j + if (interactive == YES) { + i = EC_ORDER(ec) + call printf ("Coordinate shift=%5f, rms=%5f") + call pargd (shft / i) + if (j == 1) + call pargd (INDEFD) + else + call pargd (sqrt (rms - shft ** 2) / i) + } + shft = shft + ecf_getd ("shift") + call ecf_setd ("shift", shft) + EC_NEWECF(ec) = YES + EC_NEWGRAPH(ec) = YES + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfcolon.x b/noao/onedspec/ecidentify/ecffit/ecfcolon.x new file mode 100644 index 00000000..4307335b --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfcolon.x @@ -0,0 +1,102 @@ +include <error.h> +include <gset.h> + +# List of colon commands +define CMDS "|show|function|xorder|yorder|niterate|lowreject|highreject|" + +define SHOW 1 # Show parameters +define FUNCTION 2 # Set or show function type +define XORDER 3 # Set or show x order of function +define YORDER 4 # Set or show y order of function +define NITERATE 5 # Set or show rejection iterations +define LOW 6 # Set or show low rejection threshold +define HIGH 7 # Set or show high rejection threshold + +# ECF_COLON -- Processes colon commands. + +procedure ecf_colon (cmdstr, gp) + +char cmdstr[ARB] # Command string +pointer gp # GIO pointer + +double dval +int ncmd, ival +int nscan(), strdic() +include "ecffit.com" + +begin + # Use formated scan to parse the command string. + # The first word is the command and it may be minimum match + # abbreviated with the list of commands. + + call sscan (cmdstr) + call gargwrd (ecfstr, SZ_LINE) + ncmd = strdic (ecfstr, ecfstr, SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: # :show - Show the values of the fitting parameters. + call gdeactivate (gp, AW_CLEAR) + call printf ("function %s\nxorder %d\nyorder %d\n") + call pargstr (function) + call pargi (xorder) + call pargi (yorder) + call printf ("niterate %d\nlowreject %g\nhighreject\nnreject %d\n") + call pargi (niterate) + call pargd (low) + call pargd (high) + call pargi (nreject) + call printf ("slope %d\noffset %d\nshift %g\n") + call pargi (slope) + call pargi (offset) + call pargd (shift) + call printf ("rms %g\n") + call pargd (rms) + call greactivate (gp, AW_PAUSE) + case FUNCTION: # :function - List or set the fitting function. + call gargwrd (ecfstr, SZ_LINE) + if (nscan() == 1) { + call printf ("function = %s\n") + call pargstr (function) + } else { + iferr (call ecf_sets ("function", ecfstr)) + call erract (EA_WARN) + } + case XORDER: # xorder: List or set the function order. + call gargi (ival) + if (nscan() == 1) { + call printf ("xorder %d\n") + call pargi (xorder) + } else + xorder = ival + case YORDER: # yorder: List or set the function order. + call gargi (ival) + if (nscan() == 1) { + call printf ("yorder %d\n") + call pargi (yorder) + } else + yorder = ival + case NITERATE: # niterate: List or set rejection iterations. + call gargi (ival) + if (nscan() == 1) { + call printf ("niterate %d\n") + call pargi (niterate) + } else + niterate = ival + case LOW: # low: List or set low rejection threshold. + call gargd (dval) + if (nscan() == 1) { + call printf ("lowreject %g\n") + call pargd (low) + } else + low = dval + case HIGH: # highreject: List or set high rejection threshold. + call gargd (dval) + if (nscan() == 1) { + call printf ("highreject %g\n") + call pargd (high) + } else + high = dval + default: + call printf ("Unrecognized or ambiguous command\007") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfeval.x b/noao/onedspec/ecidentify/ecffit/ecfeval.x new file mode 100644 index 00000000..1901522f --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfeval.x @@ -0,0 +1,68 @@ +# ECF_EVAL -- Evaluate wavelength at a given order and pixel position. + +double procedure ecf_eval (ecf, order, x) + +pointer ecf # GSURFIT pointer +int order # Order +double x # X point + +int ecf_oeval() +double y, dgseval() +include "ecffit.com" + +begin + y = ecf_oeval (ecf, order) + if (ecf == NULL) + return (x + shift / y) + else + return ((dgseval (ecf, x, y) + shift) / y) +end + + +# ECF_VECTOR -- Evaluate echelle dispersion function for a vector of points of +# the same order. + +procedure ecf_vector (ecf, order, x, fit, npts) + +pointer ecf # GSURFIT pointer +int order # Order +double x[npts] # X points +double fit[npts] # Fitted points +int npts # Number of points + +double yval +pointer sp, y +int ecf_oeval() +include "ecffit.com" + +begin + call smark (sp) + call salloc (y, npts, TY_DOUBLE) + + yval = ecf_oeval (ecf, order) + if (ecf == NULL) + call amovd (x, fit, npts) + else { + call amovkd (yval, Memd[y], npts) + call dgsvector (ecf, x, Memd[y], fit, npts) + call adivkd (fit, yval, fit, npts) + } + if (shift != 0.) + call aaddkd (fit, shift / yval, fit, npts) + + call sfree (sp) +end + + +# ECF_OEVAL -- Evaluate the fit order. + +int procedure ecf_oeval (ecf, order) + +pointer ecf # GSURFIT pointer +int order # User order + +include "ecffit.com" + +begin + return (slope * order + offset) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.com b/noao/onedspec/ecidentify/ecffit/ecffit.com new file mode 100644 index 00000000..61f3104a --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.com @@ -0,0 +1,23 @@ +# Common parameters. + +char function[SZ_FNAME] # Fitting function +char ecfstr[SZ_LINE] # Working char string +int gstype # Surface function type +int xorder # X order of surface function +int yorder # Y order of surface function +int niterate # Number of rejection iterations +int nreject # Number of rejected points +int xtype # X axis type +int ytype # Y axis type +int slope # Slope of order +int offset # Order offset of fit +double low, high # Low and high rejection thresholds +double xmin, xmax # X range +double ymin, ymax # Y range +double shift # First order shift +double rms # RMS of fit + + +common /ecfcom/ low, high, xmin, xmax, ymin, ymax, shift, rms, gstype, + xorder, yorder, niterate, nreject, xtype, ytype, slope, offset, + function, ecfstr diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.h b/noao/onedspec/ecidentify/ecffit/ecffit.h new file mode 100644 index 00000000..20825c71 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.h @@ -0,0 +1,20 @@ +define IGSPARAMS 7 + +define FEATURE 1 +define X 2 +define Y 3 +define Z 4 +define W 5 +define S 6 +define R 7 + +define IGS_FUNCTION 1 +define IGS_XORDER 2 +define IGS_YORDER 3 +define IGS_XMIN 4 +define IGS_XMAX 5 +define IGS_YMIN 6 +define IGS_YMAX 7 +define IGS_OFFSET 8 + +define SFTYPES "|chebyshev|legendre|" # Surface types diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.key b/noao/onedspec/ecidentify/ecffit/ecffit.key new file mode 100644 index 00000000..f24407b9 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.key @@ -0,0 +1,53 @@ + ECHELLE DISPERSION FUNCTION FITTING KEYS + + +CURSOR KEY SUMMARY + +? Help c Print coordinates d Delete point +f Fit dispersion o Fit with fixed order offset q Quit +r Redraw graph u Undelete point w Window graph +x Set ordinate y Set abscissa I Interrupt + + +COLON COMMAND SUMMARY + +:show :function [value] :highreject [value] :lowreject [value] +:niterate [value] :xorder [value] :yorder [value] + + +CURSOR KEYS + +? Print this list of cursor keys +c Print cursor coordinates +d Delete the nearest undeleted point to the cursor +f Fit dispersion function including determining the order offset +o Fit dispersion function with the order offset fixed +q Quit and return to the spectrum display +r Redraw the graph +u Undelete the nearest deleted point to the cursor (may be outside the window) +w Window the graph (type ? to the window prompt for more help) +x Set the quantity plotted along the ordinate (x axis) +y Set the quantity plotted along the abscissa (y axis) +I Interrupt the task immediately + + +COLON COMMANDS + +:show Print current function and orders +:function [value] Print or set the function type (chebyshev|legendre) +:highreject [value] Print or set high rejection limit +:lowreject [value] Print or set high rejection limit +:niterate [value] Print or set number of rejection iterations +:xorder [value] Print or set the order for the dispersion dependence +:yorder [value] Print or set the order for the echelle order dependence + + +The dispersion function fitted is given by a two dimensional function +(either chebyshev or legendre) of the pixel position along the +dispersion of an order (called x) and the order number (called y). The +order number is determined from the aperture number by an offset and +direction of increasing order number. The basic order dependence is +separated from the surface function as given below. + + y = offset +/- aperture + wavelength = f (x, y) / y diff --git a/noao/onedspec/ecidentify/ecffit/ecffit.x b/noao/onedspec/ecidentify/ecffit/ecffit.x new file mode 100644 index 00000000..408a1b77 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecffit.x @@ -0,0 +1,193 @@ +include <error.h> +include <pkg/gtools.h> + +define HELP "noao$onedspec/ecidentify/ecffit/ecffit.key" +define PROMPT "fitcoords surface fitting options" + +# EC_FIT -- Echelle dispersion fitting. +# +# X - Pixel coordinates along dispersion +# Y - Relative order number +# Z - Wavelength + +procedure ecf_fit (ecf, gp, gt, xd, yd, zd, wd, npts, fixedorder) + +pointer ecf # GSURFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +double xd[npts] # Pixel coordinates along dispersion +double yd[npts] # Order number +double zd[npts] # Wavelength +double wd[npts] # Weights +int npts # Number of points +int fixedorder # Fixed order? + +real wx, wy +int wcs, key +int i, newgraph +pointer sp, wd1, rd, xr, yr +char cmd[SZ_LINE] + +int ecf_nearest() +int clgcur(), scan(), nscan() +errchk ecf_solve() +include "ecffit.com" + +begin + # Allocate residuals and weights with rejected points arrays + call smark (sp) + call salloc (wd1, npts, TY_DOUBLE) + call salloc (rd, npts, TY_DOUBLE) + call amovd (wd, Memd[wd1], npts) + + # Compute a solution and return if not interactive. + if (gp == NULL) { + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + do i = 1, npts + if (Memd[wd1+i-1] != wd[i]) + wd[i] = -1. + call sfree (sp) + return + } + + # Allocate real graph vectors. + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + + # Read cursor commands. + key = 'f' + repeat { + switch (key) { + case 'o': + call printf ("Order offset (%d): ") + call pargi (offset) + call flush (STDOUT) + if (scan() != EOF) { + call gargi (i) + if (nscan() == 1) + offset = i + call amovd (wd, Memd[wd1], npts) + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + YES) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + YES) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + } + + case '?': # Print help text. + call gpagefile (gp, HELP, PROMPT) + + case ':': # List or set parameters + if (cmd[1] == '/') + call gt_colon (cmd, gp, gt, newgraph) + else + call ecf_colon (cmd, gp) + + case 'x': # Set ordinate + call printf ("Ordinate - ") + call printf ( + "(p)ixel, (o)rder, (w)avelength, (r)esidual, (v)elocity: ") + if (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + break + + if (key != xtype) { + if (key=='p'||key=='o'||key=='w'||key=='r'||key=='v') { + xtype = key + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_title (gt) + newgraph = YES + } else + call printf ("\007") + } + + case 'y': # Set abscissa + call printf ("Abscissa - ") + call printf ( + "(p)ixel, (o)rder, (w)avelength, (r)esidual, (v)elocity: ") + if(clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + break + + if (key != ytype) { + if (key=='p'||key=='o'||key=='w'||key=='r'||key=='v') { + ytype = key + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + } else + call printf ("\007") + } + + case 'r': # Redraw + newgraph = YES + + case 'c': # Cursor coordinates + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + call printf ("%10.2g %d %10.8g\n") + call pargd (xd[i]) + call pargd (yd[i]) + call pargd (zd[i]) + + case 'd': # Delete + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + if (i > 0) + Memd[wd1+i-1] = wd[i] + + case 'u': # Undelete + i = ecf_nearest (gp, gt, wx, wy, wcs, key, Memr[xr], Memr[yr], + wd, npts) + if (i > 0) + Memd[wd1+i-1] = wd[i] + + case 'f': # Fit + call amovd (wd, Memd[wd1], npts) + call ecf_solve (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_reject (ecf, xd, yd, zd, Memd[wd1], Memd[rd], npts, + fixedorder) + call ecf_gdata (ecf, xtype, xd, yd, zd, Memd[rd], + Memr[xr], npts) + call ecf_gdata (ecf, ytype, xd, yd, zd, Memd[rd], + Memr[yr], npts) + call ecf_title (gt) + newgraph = YES + + case 'w': # Window graph + call gt_window (gt, gp, "cursor", newgraph) + + case 'q': # Quit + break + + case 'I': # Interrupt + call fatal (0, "Interrupt") + + default: # Ring the bell. + call printf ("\07\n") + } + + if (newgraph == YES) { + call ecf_graph (gp, gt, Memr[xr], Memr[yr], wd, Memd[wd1], npts) + newgraph = NO + } + } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + + do i = 1, npts + if (Memd[wd1+i-1] != wd[i]) + wd[i] = -1. + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfgdata.x b/noao/onedspec/ecidentify/ecffit/ecfgdata.x new file mode 100644 index 00000000..eebb34d6 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfgdata.x @@ -0,0 +1,37 @@ +include <pkg/gtools.h> + +# ECF_GDATA -- Get graph data for the specified axis type from the fitting data. + +procedure ecf_gdata (ecf, type, x, y, z, r, data, npts) + +pointer ecf # GSURFIT pointer +int type # Axis type +double x[npts] # X fit data +double y[npts] # Y fit data +double z[npts] # Z fit data +double r[npts] # Residuals +real data[npts] # Graph data +int npts # Number of points + +pointer sp, v +include "ecffit.com" + +begin + switch (type) { + case 'p': + call achtdr (x, data, npts) + case 'o': + call achtdr (y, data, npts) + case 'w': + call achtdr (z, data, npts) + case 'r': + call achtdr (r, data, npts) + case 'v': + call smark (sp) + call salloc (v, npts, TY_DOUBLE) + call adivd (r, z, Memd[v], npts) + call amulkd (Memd[v], 300000.D0, Memd[v], npts) + call achtdr (Memd[v], data, npts) + call sfree (sp) + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfget.x b/noao/onedspec/ecidentify/ecffit/ecfget.x new file mode 100644 index 00000000..025059df --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfget.x @@ -0,0 +1,84 @@ +# ECF_GETI -- Get the value of an integer parameter. + +int procedure ecf_geti (param) + +char param[ARB] # ECF parameter + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|slope|offset|xorder|yorder|niterate|") + switch (i) { + case 1: + return (slope) + case 2: + return (offset) + case 3: + return (xorder) + case 4: + return (yorder) + case 5: + return (niterate) + default: + call error (0, "ecf_geti: Unknown parameter") + } +end + + +# ECF_GETS -- Get the value of a string parameter. + +procedure ecf_gets (param, str, maxchar) + +char param[ARB] # ECF parameter +char str[maxchar] # String +int maxchar # Maximum number of characters + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, "|function|") + switch (i) { + case 1: + call strcpy (function, str, maxchar) + default: + call error (0, "ecf_gets: Unknown parameter") + } +end + + +# ECF_GETD -- Get the values of double valued fitting parameters. + +double procedure ecf_getd (param) + +char param[ARB] # ECF parameter + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|xmin|xmax|ymin|ymax|shift|rms|low|high|") + switch (i) { + case 1: + return (xmin) + case 2: + return (xmax) + case 3: + return (ymin) + case 4: + return (ymax) + case 5: + return (shift) + case 6: + return (rms) + case 7: + return (low) + case 8: + return (high) + default: + call error (0, "ecf_gets: Unknown parameter") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfgraph.x b/noao/onedspec/ecidentify/ecffit/ecfgraph.x new file mode 100644 index 00000000..22749527 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfgraph.x @@ -0,0 +1,50 @@ +include <gset.h> +include <mach.h> +include <pkg/gtools.h> + +# ECF_GRAPH -- Graph the fitted data. + +procedure ecf_graph (gp, gt, x, y, w, rej, npts) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +real x[npts] # X data +real y[npts] # Y data +double w[npts] # Weights +double rej[npts] # Rejected points +int npts # Number of pts points + +int i +real xsize, ysize, ymin, ymax, gt_getr() + +begin + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gclear (gp) + + ymin = MAX_REAL + ymax = -MAX_REAL + do i = 1, npts + if (w[i] > 0.) { + ymin = min (ymin, y[i]) + ymax = max (ymax, y[i]) + } + + call gascale (gp, x, npts, 1) + call gswind (gp, INDEF, INDEF, ymin, ymax) + call gt_swind (gp, gt) + call gt_labax (gp, gt) + + do i = 1, npts { + if (rej[i] == 0.) { + if (y[i] >= ymin && y[i] <= ymax) { + if (w[i] == 0.) + call gmark (gp, x[i], y[i], GM_CROSS, xsize, ysize) + else + call gmark (gp, x[i], y[i], GM_DIAMOND, xsize, ysize) + } + } else + call gmark (gp, x[i], y[i], GM_PLUS, xsize, ysize) + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfnearest.x b/noao/onedspec/ecidentify/ecffit/ecfnearest.x new file mode 100644 index 00000000..af1b1f78 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfnearest.x @@ -0,0 +1,85 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> + +# ECF_NEAREST -- Find nearest point to the cursor. + +int procedure ecf_nearest (gp, gt, wx, wy, wcs, key, x, y, w, npts) + +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +real wx, wy # Cursor coordinates +int wcs # WCS +int key # Nearest key +real x[npts] # Data points +real y[npts] # Data points +double w[npts] # Weight +int npts # Number of data points + +int i, j +real r2, r2min, x0, y0, xsize, ysize, gt_getr() + +begin + call gctran (gp, wx, wy, wx, wy, wcs, 0) + r2min = MAX_REAL + j = 0 + + switch (key) { + case 'c': + do i = 1, npts { + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + call gscur (gp, x[j], y[j]) + case 'd': + do i = 1, npts { + if (w[i] == 0.) + next + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + if (j > 0) { + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, x[j], y[j], GM_PLUS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, x[j], y[j], GM_CROSS, xsize, ysize) + w[j] = 0. + call gscur (gp, x[j], y[j]) + } + case 'u': + do i = 1, npts { + if (w[i] != 0.) + next + call gctran (gp, x[i], y[i], x0, y0, wcs, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + if (j > 0) { + xsize = gt_getr (gt, GTXSIZE) + ysize = gt_getr (gt, GTYSIZE) + + call gseti (gp, G_PMLTYPE, 0) + call gmark (gp, x[j], y[j], GM_CROSS, xsize, ysize) + call gseti (gp, G_PMLTYPE, 1) + call gmark (gp, x[j], y[j], GM_PLUS, xsize, ysize) + w[j] = 1. + call gscur (gp, x[j], y[j]) + } + } + + return (j) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfreject.x b/noao/onedspec/ecidentify/ecffit/ecfreject.x new file mode 100644 index 00000000..a772069e --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfreject.x @@ -0,0 +1,53 @@ +include <mach.h> + +# ECF_REJECT -- Reject points with large residuals from the fit. + +procedure ecf_reject (ecf, x, y, z, w, r, npts, fixedorder) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points +int fixedorder # Fixed order? + +int i, j, newreject +double low_cut, high_cut +include "ecffit.com" + +begin + # Return if rejection is not desired. + nreject = 0 + if (niterate == 0 || (low == 0. && high == 0.)) + return + + # Reject points. + do i = 1, niterate { + if (low > 0.) + low_cut = -low * rms + else + low_cut = -MAX_REAL + if (high > 0.) + high_cut = high * rms + else + high_cut = MAX_REAL + + newreject = 0 + do j = 1, npts { + if (w[j] == 0.) + next + if ((r[j] > high_cut) || (r[j] < low_cut)) { + w[j] = 0. + newreject = newreject + 1 + } + } + + if (newreject == 0) + break + + call ecf_solve (ecf, x, y, z, w, r, npts, fixedorder) + nreject = nreject + newreject + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfrms.x b/noao/onedspec/ecidentify/ecffit/ecfrms.x new file mode 100644 index 00000000..1140dc29 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfrms.x @@ -0,0 +1,26 @@ +# ECF_RMS -- Compute the rms with deleted points ignored. + +double procedure ecf_rms (r, w, npts) + +double r[npts] # Residuals +double w[npts] # Weights +int npts # Number of points + +int i, n +double rms + +begin + n = 0 + rms = 0. + do i = 1, npts { + if (w[i] == 0.) + next + n = n + 1 + rms = rms + r[i] * r[i] + } + if (n > 0) + rms = sqrt (rms / n) + else + rms = INDEFD + return (rms) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfset.x b/noao/onedspec/ecidentify/ecffit/ecfset.x new file mode 100644 index 00000000..4b6402b1 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfset.x @@ -0,0 +1,92 @@ +# ECF_SETS -- Set the values of string valued fitting parameters. + +procedure ecf_sets (param, str) + +char param[ARB] # Parameter to be set +char str[ARB] # String value + +char temp[10] +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, temp, 10, "|function|") + switch (i) { + case 1: + i = strdic (str, str, SZ_FNAME, "|chebyshev|legendre|") + if (i == 0) + call error (0, "Unknown function type") + call strcpy (str, function, SZ_LINE) + gstype = i + default: + call error (0, "ecf_sets: Unknown parameter") + } +end + + +# ECF_SETI -- Set the values of integer valued fitting parameters. + +procedure ecf_seti (param, ival) + +char param[ARB] # Parameter to be set +int ival # Integer value + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|slope|offset|xorder|yorder|xtype|ytype|niterate|") + switch (i) { + case 1: + slope = ival + case 2: + offset = ival + case 3: + xorder = ival + case 4: + yorder = ival + case 5: + xtype = ival + case 6: + ytype = ival + case 7: + niterate = max (0, ival) + default: + call error (0, "ecf_seti: Unknown parameter") + } +end + + +# ECF_SETD -- Set the values of double valued fitting parameters. + +procedure ecf_setd (param, dval) + +char param[ARB] # Parameter to be set +double dval # Double value + +int i, strdic() +include "ecffit.com" + +begin + i = strdic (param, ecfstr, SZ_LINE, + "|xmin|xmax|ymin|ymax|shift|low|high|") + switch (i) { + case 1: + xmin = dval + case 2: + xmax = dval + case 3: + ymin = dval + case 4: + ymax = dval + case 5: + shift = dval + case 6: + low = max (0.D0, dval) + case 7: + high = max (0.D0, dval) + default: + call error (0, "ecf_setd: Unknown parameter") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfshift.x b/noao/onedspec/ecidentify/ecffit/ecfshift.x new file mode 100644 index 00000000..75655703 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfshift.x @@ -0,0 +1,55 @@ +# ECF_GSHIFT -- Return the shift for the given order. + +double procedure ecf_gshift (ecf, order) + +pointer ecf # GSURFIT pointer +int order # User order + +include "ecffit.com" + +begin + return (shift / (slope * order + offset)) +end + + +# ECF_PSHIFT -- Put the shift for the given order. + +procedure ecf_pshift (ecf, order, shft) + +pointer ecf # GSURFIT pointer +int order # User order +double shft # Shift at given order + +include "ecffit.com" + +begin + shift = shft * (slope * order + offset) +end + + +procedure ecf_vector (ecf, order, x, fit, npts) + +pointer ecf # GSURFIT pointer +int order # Order +double x[npts] # X points +double fit[npts] # Fitted points +int npts # Number of points + +double yval +pointer sp, y + +include "ecffit.com" + +begin + call smark (sp) + call salloc (y, npts, TY_DOUBLE) + + yval = slope * order + offset + call amovkd (yval, Memd[y], npts) + call dgsvector (ecf, x, Memd[y], fit, npts) + call adivkd (fit, yval, fit, npts) + if (shift != 0.) + call aaddkd (fit, shift / yval, fit, npts) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecfsolve.x b/noao/onedspec/ecidentify/ecffit/ecfsolve.x new file mode 100644 index 00000000..1c844e76 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecfsolve.x @@ -0,0 +1,196 @@ +include <mach.h> +include <math/gsurfit.h> + +define ECFTYPES "|chebyshev|legendre|" # Fit types + + +# ECF_SOLVE -- Fit +# +# f(x, slope*y+offset) = (y+slope*offset)*z +# +# with offset minimizing the RMS. + +procedure ecf_solve (ecf, x, y, z, w, r, npts, fixedorder) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points +int fixedorder # Fixed order? + +int i, j, k, err +double ya, yb, newrms, ecf_rms() +pointer sp, y1, ecf1 +errchk ecf_solve1 +include "ecffit.com" +define fit_ 99 + +begin + if (fixedorder == YES) { + call ecf_solve1 (ecf, x, y, z, w, r, npts) + return + } + + call smark (sp) + call salloc (y1, npts, TY_DOUBLE) + + # Determine if the orders are reversed. + j = 1 + k = 1 + do i = 1, npts { + if (z[i] < z[j]) + j = i + if (z[i] > z[k]) + k = i + } + if (y[j] >= y[k]) { + slope = 1 + offset = max (offset, int(1. - ymin)) + } else { + slope = -1 + offset = max (offset, int(1. + ymax)) + } + + call dgsfree (ecf) + shift = 0. + + rms = MAX_DOUBLE + j = 1 + k = 0 + + for (i=offset;;i=i+j) { + if (slope == 1) { + ya = i + ymin + yb = i + ymax + } else { + ya = i - ymax + yb = i - ymin + } + if (ya < 1.) + break + + call altmd (y, Memd[y1], npts, double(slope), double(i)) + call amuld (Memd[y1], z, r, npts) + +fit_ call dgsinit (ecf1, gstype, xorder, yorder, YES, xmin, xmax, ya, yb) + call dgsfit (ecf1, x, Memd[y1], r, w, npts, WTS_USER, err) + + if (err != OK) { + if (xorder > 2 || yorder > 2) { + call dgsfree (ecf) + xorder = max (2, xorder - 1) + yorder = max (2, yorder - 1) + goto fit_ + } + + switch (err) { + case SINGULAR: + call dgsfree (ecf) + ecf = ecf1 + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call sfree (sp) + call error (0, "No degrees of freedom") + } + } + + call dgsvector (ecf1, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + + newrms = ecf_rms (r, w, npts) + k = k + 1 + + if (newrms / rms < 0.999) { + call dgsfree (ecf) + ecf = ecf1 + offset = i + rms = newrms + } else { + call dgsfree (ecf1) + if (k > 2) + break + i = offset + j = -j + } + } + + call altmd (y, Memd[y1], npts, double(slope), double(offset)) + call dgsvector (ecf, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + + call sfree (sp) + +end + + +# ECF_SOLVE1 -- Fit f(x, y+offset) = (y+offset)*z with offset fixed. + +procedure ecf_solve1 (ecf, x, y, z, w, r, npts) + +pointer ecf # GSURFIT pointer +double x[npts] # X points +double y[npts] # Y points +double z[npts] # Z points +double w[npts] # Weights +double r[npts] # Residuals +int npts # Number of points + +int err +pointer sp, y1 +double ya, yb, ecf_rms() +include "ecffit.com" +define fit_ 99 + +begin + call smark (sp) + call salloc (y1, npts, TY_DOUBLE) + + call dgsfree (ecf) + shift = 0. + + if (slope == 1) { + offset = max (offset, int(1. - ymin)) + ya = offset + ymin + yb = offset + ymax + } else { + offset = max (offset, int(1. + ymax)) + ya = offset - ymax + yb = offset - ymin + } + + call altmd (y, Memd[y1], npts, double (slope), double (offset)) + call amuld (Memd[y1], z, r, npts) + +fit_ call dgsinit (ecf, gstype, xorder, yorder, YES, xmin, xmax, + min (ya, yb), max (ya, yb)) + call dgsfit (ecf, x, Memd[y1], r, w, npts, WTS_USER, err) + + if (err != OK) { + if (xorder > 2 || yorder > 2) { + call dgsfree (ecf) + xorder = max (2, xorder - 1) + yorder = max (2, yorder - 1) + goto fit_ + } + + switch (err) { + case SINGULAR: + call eprintf ("Singular solution\n") + case NO_DEG_FREEDOM: + call sfree (sp) + call error (0, "No degrees of freedom") + } + } + + call dgsvector (ecf, x, Memd[y1], r, npts) + call adivd (r, Memd[y1], r, npts) + call asubd (z, r, r, npts) + rms = ecf_rms (r, w, npts) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecffit/ecftitle.x b/noao/onedspec/ecidentify/ecffit/ecftitle.x new file mode 100644 index 00000000..3b754f31 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/ecftitle.x @@ -0,0 +1,48 @@ +include <pkg/gtools.h> + +# ECF_TITLE -- Set the GTOOLS parameter string. + +procedure ecf_title (gt) + +pointer gt # GTOOLS pointer + +include "ecffit.com" + +begin + call sprintf (ecfstr, SZ_LINE, + "Function=%s, xorder=%d, yorder=%d, slope=%d, offset=%d, rms=%6g") + call pargstr (function) + call pargi (xorder) + call pargi (yorder) + call pargi (slope) + call pargi (offset) + call pargd (rms) + call gt_sets (gt, GTPARAMS, ecfstr) + call gt_sets (gt, GTTITLE, "Echelle Dispersion Function Fitting") + + switch (xtype) { + case 'p': + call gt_sets (gt, GTXLABEL, "Pixel") + case 'o': + call gt_sets (gt, GTXLABEL, "Order") + case 'w': + call gt_sets (gt, GTXLABEL, "Wavelength") + case 'r': + call gt_sets (gt, GTXLABEL, "Residual") + case 'v': + call gt_sets (gt, GTXLABEL, "Velocity") + } + + switch (ytype) { + case 'p': + call gt_sets (gt, GTYLABEL, "Pixel") + case 'o': + call gt_sets (gt, GTYLABEL, "Order") + case 'w': + call gt_sets (gt, GTYLABEL, "Wavelength") + case 'r': + call gt_sets (gt, GTYLABEL, "Residual") + case 'v': + call gt_sets (gt, GTYLABEL, "Velocity") + } +end diff --git a/noao/onedspec/ecidentify/ecffit/mkpkg b/noao/onedspec/ecidentify/ecffit/mkpkg new file mode 100644 index 00000000..40324cb8 --- /dev/null +++ b/noao/onedspec/ecidentify/ecffit/mkpkg @@ -0,0 +1,21 @@ +# Echelle Dispersion Fitting Package + +$checkout libpkg.a ../../ +$update libpkg.a +$checkin libpkg.a ../../ +$exit + +libpkg.a: + ecfcolon.x ecffit.com <error.h> <gset.h> + ecfeval.x ecffit.com + ecffit.x ecffit.com <error.h> <pkg/gtools.h> + ecfgdata.x ecffit.com <pkg/gtools.h> + ecfget.x ecffit.com + ecfgraph.x <gset.h> <mach.h> <pkg/gtools.h> + ecfnearest.x <gset.h> <mach.h> <pkg/gtools.h> + ecfreject.x ecffit.com <mach.h> + ecfrms.x + ecfset.x ecffit.com + ecfsolve.x ecffit.com <mach.h> <math/gsurfit.h> + ecftitle.x ecffit.com <pkg/gtools.h> + ; diff --git a/noao/onedspec/ecidentify/ecfitdata.x b/noao/onedspec/ecidentify/ecfitdata.x new file mode 100644 index 00000000..998f5057 --- /dev/null +++ b/noao/onedspec/ecidentify/ecfitdata.x @@ -0,0 +1,146 @@ +include <pkg/gtools.h> +include <smw.h> +include <units.h> +include "ecidentify.h" + +# EC_FITDATA -- Compute fit coordinates from pixel coordinates. + +procedure ec_fitdata (ec) + +pointer ec # ID pointer + +int i, ecf_oeval() + +begin + call mfree (EC_FITDATA(ec), TY_DOUBLE) + call malloc (EC_FITDATA(ec), EC_NCOLS(ec)*EC_NLINES(ec), TY_DOUBLE) + + do i = 1, EC_NLINES(ec) { + call ec_gline (ec, i) + if (EC_ECF(ec) == NULL) { + if (DC(EC_SH(ec)) != DCNO && EC_UN(ec) != NULL) + iferr (call shdr_units (EC_SH(ec), UN_UNITS(EC_UN(ec)))) + ; + call achtrd (Memr[SX(EC_SH(ec))], FITDATA(ec,1), EC_NPTS(ec)) + call gt_sets (EC_GT(ec), GTXLABEL, LABEL(EC_SH(ec))) + call gt_sets (EC_GT(ec), GTXUNITS, UNITS(EC_SH(ec))) + } else { + ORDERS(ec,i) = ecf_oeval (EC_ECF(ec), APS(ec,i)) + call ecf_vector (EC_ECF(ec), APS(ec,i), PIXDATA(ec,1), + FITDATA(ec,1), EC_NPTS(ec)) + if (EC_UN(ec) == NULL) { + call gt_sets (EC_GT(ec), GTXLABEL, LABEL(EC_SH(ec))) + call gt_sets (EC_GT(ec), GTXUNITS, UNITS(EC_SH(ec))) + } else { + call gt_sets (EC_GT(ec), GTXLABEL, UN_LABEL(EC_UN(ec))) + call gt_sets (EC_GT(ec), GTXUNITS, UN_UNITS(EC_UN(ec))) + } + } + } + + call ec_gline (ec, EC_LINE(ec)) + EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec)) +end + + +# EC_FITFEATURES -- Compute fit coordinates for features. + +procedure ec_fitfeatures (ec) + +pointer ec # ID pointer + +int i, ec_line() +double ec_fitpt() + +begin + if (EC_NFEATURES(ec) < 1) + return + + do i = 1, EC_NFEATURES(ec) { + LINE(ec,i) = ec_line (ec, APN(ec,i)) + ORDER(ec,i) = ORDERS(ec,LINE(ec,i)) + FIT(ec,i) = ec_fitpt (ec, APN(ec,i), PIX(ec,i)) + } +end + + +# EC_FITPT -- Compute fit coordinates from pixel coordinates. + +double procedure ec_fitpt (ec, order, pix) + +pointer ec # ID pointer +int order # Order +double pix # Pixel coordinate + +double fit, ecf_eval(), smw_c1trand(), shdr_lw() + +begin + if (EC_ECF(ec) == NULL) { + fit = smw_c1trand (EC_PL(ec), pix) + fit = shdr_lw (EC_SH(ec), fit) + } else + fit = ecf_eval (EC_ECF(ec), order, pix) + + return (fit) +end + + +# EC_FITTOPIX -- Transform fit coordinate to pixel coordinate. + +define DXMIN .01 + +double procedure ec_fittopix (ec, fitcoord) + +pointer ec # ID pointer +double fitcoord # Fit coordinate to be transformed +double pixcoord # Pixel coordinate returned + +int i, n +double dx + +double ec_fitpt(), smw_c1trand() + +begin + n = EC_NPTS(ec) + if (FITDATA(ec,1) < FITDATA(ec,n)) { + if ((fitcoord<FITDATA(ec,1)) || (fitcoord>FITDATA(ec,n))) + return (INDEFD) + + for (i = 1; fitcoord > FITDATA(ec,i); i = i + 1) + ; + + if (FITDATA(ec,i) == fitcoord) + return (double (i)) + + pixcoord = smw_c1trand (EC_LP(ec), double(i-.5)) + dx = smw_c1trand (EC_LP(ec), double(i+.5)) - pixcoord + while (dx > DXMIN) { + dx = dx / 2 + if (ec_fitpt (ec, EC_AP(ec), pixcoord) < fitcoord) + pixcoord = pixcoord + dx + else + pixcoord = pixcoord - dx + } + } else { + if ((fitcoord<FITDATA(ec,n)) || (fitcoord>FITDATA(ec,1))) + return (INDEFD) + + for (i = 1; fitcoord < FITDATA(ec,i); i = i + 1) + ; + + if (FITDATA(ec,i) == fitcoord) + return (double (i)) + + pixcoord = smw_c1trand (EC_LP(ec), double(i-.5)) + dx = smw_c1trand (EC_LP(ec), double(i+.5)) - pixcoord + while (dx > DXMIN) { + dx = dx / 2 + if (ec_fitpt (ec, EC_AP(ec), pixcoord) < fitcoord) + pixcoord = pixcoord - dx + else + pixcoord = pixcoord + dx + } + } + + return (pixcoord) +end diff --git a/noao/onedspec/ecidentify/ecgdata.x b/noao/onedspec/ecidentify/ecgdata.x new file mode 100644 index 00000000..1087d38c --- /dev/null +++ b/noao/onedspec/ecidentify/ecgdata.x @@ -0,0 +1,74 @@ +include <imhdr.h> +include <imio.h> +include <pkg/gtools.h> +include <smw.h> +include <units.h> +include "ecidentify.h" + +# EC_GDATA -- Get image data. + +procedure ec_gdata (ec) + +pointer ec # ID pointer + +int i, j +pointer im, mw, sh, sp, str1, str2 + +double smw_c1trand() +pointer immap(), smw_openim(), smw_sctran() +errchk immap, smw_openim, shdr_open + +begin + # Map the image. + im = immap (Memc[EC_IMAGE(ec)], READ_ONLY, 0) + + # Free previous data + do i = 1, EC_NLINES(ec) + call shdr_close (SH(ec,i)) + call mfree (EC_SHS(ec), TY_POINTER) + call mfree (EC_PIXDATA(ec), TY_DOUBLE) + + # Set MWCS + mw = smw_openim (im) + EC_LP(ec) = smw_sctran (mw, "logical", "physical", 1) + EC_PL(ec) = smw_sctran (mw, "physical", "logical", 1) + + # Allocate new vectors. + EC_NCOLS(ec) = IM_LEN(im, 1) + EC_NLINES(ec) = IM_LEN(im, 2) + call calloc (EC_SHS(ec), EC_NLINES(ec), TY_POINTER) + call malloc (EC_PIXDATA(ec), EC_NCOLS(ec)*EC_NLINES(ec), TY_DOUBLE) + + # Set the coordinates. + sh = NULL + do j = 1, EC_NLINES(ec) { + call shdr_open (im, mw, j, 1, INDEFI, SHDATA, sh) + if (EC_UN(ec) != NULL) + iferr (call shdr_units (sh, UN_UNITS(EC_UN(ec)))) + ; + if (j != EC_NLINES(ec)) + call shdr_copy (sh, SH(ec,j), NO) + else + SH(ec,j) = sh + call ec_gline (ec, j) + do i = 1, EC_NPTS(ec) + PIXDATA(ec,i) = smw_c1trand (EC_LP(ec), double(i)) + } + EC_LINE(ec) = 1 + call ec_gline (ec, EC_LINE(ec)) + EC_AP(ec) = APS(ec,EC_LINE(ec)) + EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec)) + + # Set graph title. + call smark (sp) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + + call sprintf (Memc[str1], SZ_LINE, "ecidentify %s: %s") + call pargstr (Memc[EC_IMAGE(ec)]) + call pargstr (IM_TITLE(im)) + call gt_sets (EC_GT(ec), GTTITLE, Memc[str1]) + + call imunmap (im) + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecgetim.x b/noao/onedspec/ecidentify/ecgetim.x new file mode 100644 index 00000000..cbcb244e --- /dev/null +++ b/noao/onedspec/ecidentify/ecgetim.x @@ -0,0 +1,17 @@ +# EC_GETIM -- Get next image name with standard image extensions removed. + +int procedure ec_getim (list, image, maxchar) + +int list # Image list +char image[maxchar] # Image name +int maxchar # Maximum number of chars in image name + +int stat, imtgetim() + +begin + stat = imtgetim (list, image, maxchar) + if (stat != EOF) + call xt_imroot (image, image, maxchar) + + return (stat) +end diff --git a/noao/onedspec/ecidentify/ecgline.x b/noao/onedspec/ecidentify/ecgline.x new file mode 100644 index 00000000..7d3f9e16 --- /dev/null +++ b/noao/onedspec/ecidentify/ecgline.x @@ -0,0 +1,20 @@ +include <smw.h> +include "ecidentify.h" + +# EC_GLINE -- Get line of data. + +procedure ec_gline (ec, line) + +pointer ec # EC pointer +int line # Image line + +begin + if (IS_INDEFI(line)) + return + + EC_SH(ec) = SH(ec,line) + EC_NPTS(ec) = SN(EC_SH(ec)) + EC_IMLINE(ec) = SY(EC_SH(ec)) + EC_PIXLINE(ec) = EC_PIXDATA(ec) + (line - 1) * EC_NCOLS(ec) + EC_FITLINE(ec) = EC_FITDATA(ec) + (line - 1) * EC_NCOLS(ec) +end diff --git a/noao/onedspec/ecidentify/ecgraph.x b/noao/onedspec/ecidentify/ecgraph.x new file mode 100644 index 00000000..9eaeaa5f --- /dev/null +++ b/noao/onedspec/ecidentify/ecgraph.x @@ -0,0 +1,155 @@ +include <gset.h> +include <pkg/gtools.h> +include "ecidentify.h" + +# EC_GRAPH -- Graph image vector in which features are to be ecentified. + +procedure ec_graph (ec, gtype) + +pointer ec # ID pointer +int gtype # Graph type + +begin + switch (gtype) { + case 1: + if (IS_INDEFI (EC_AP(ec))) + call ec_graph3(ec) + else + call ec_graph1 (ec) + case 2: + call ec_graph2 (ec) + default: + call ec_graph1 (ec) + } +end + + +procedure ec_graph1 (ec) + +pointer ec # ID pointer + +int i +real xmin, xmax, ymin, ymax, dy +pointer sp, str, x, y + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (x, EC_NPTS(ec), TY_REAL) + y = EC_IMLINE(ec) + + call sprintf (Memc[str], SZ_LINE, + "Aperture %d, Image line %d, Order %d") + call pargi (EC_AP(ec)) + call pargi (EC_LINE(ec)) + call pargi (EC_ORDER(ec)) + call gt_sets (EC_GT(ec), GTPARAMS, Memc[str]) + call achtdr (FITDATA(ec,1), Memr[x], EC_NPTS(ec)) + + call gclear (EC_GP(ec)) + xmin = min (Memr[x], Memr[x+EC_NPTS(ec)-1]) + xmax = max (Memr[x], Memr[x+EC_NPTS(ec)-1]) + call alimr (Memr[y], EC_NPTS(ec), ymin, ymax) + dy = ymax - ymin + call gswind (EC_GP(ec), xmin, xmax, ymin - .2 * dy, ymax + .2 * dy) + call gt_swind (EC_GP(ec), EC_GT(ec)) + call gt_labax (EC_GP(ec), EC_GT(ec)) + call gt_plot (EC_GP(ec), EC_GT(ec), Memr[x], Memr[y], EC_NPTS(ec)) + + do i = 1, EC_NFEATURES(ec) + if (APN(ec,i) == EC_AP(ec)) + call ec_mark (ec, i) + + call sfree (sp) +end + + +# EC_GRAPH2 -- Make review graph for current feature. + +procedure ec_graph2 (ec) + +pointer ec # ID pointer + +int i, j, k +real xmin, xmax, ymin, ymax, dy +pointer sp, str, x, y + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (x, EC_NPTS(ec), TY_REAL) + y = EC_IMLINE(ec) + + call sprintf (Memc[str], SZ_LINE, "Order %d") + call pargi (EC_AP(ec)) + call gt_sets (EC_GT(ec), GTPARAMS, Memc[str]) + call achtdr (FITDATA(ec,1), Memr[x], EC_NPTS(ec)) + + xmin = real (FIT(ec,EC_CURRENT(ec))) - EC_ZWIDTH(ec) / 2. + xmax = real (FIT(ec,EC_CURRENT(ec))) + EC_ZWIDTH(ec) / 2. + + i = 0 + do k = 1, EC_NPTS(ec) { + 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 (EC_GP(ec)) + call gswind (EC_GP(ec), xmin, xmax, ymin - .2 * dy, ymax + .2 * dy) + call gt_labax (EC_GP(ec), EC_GT(ec)) + call gt_plot (EC_GP(ec), EC_GT(ec), Memr[x], Memr[y], EC_NPTS(ec)) + + do i = 1, EC_NFEATURES(ec) + if (APN(ec,i) == EC_AP(ec)) + call ec_mark (ec, i) + + call sfree (sp) +end + + +procedure ec_graph3 (ec) + +pointer ec # ID pointer + +int i, npts +real xmin, xmax, ymin, ymax, dy +pointer sp, str, x, y + +begin + npts = EC_NPTS(ec) * EC_NLINES(ec) + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (x, npts, TY_REAL) + y = EC_IMLINE(ec) + + call sprintf (Memc[str], SZ_LINE, "All orders") + call gt_sets (EC_GT(ec), GTPARAMS, Memc[str]) + call achtdr (Memd[EC_FITDATA(ec)], Memr[x], npts) + + call gclear (EC_GP(ec)) + xmin = min (Memr[x], Memr[x+npts-1]) + xmax = max (Memr[x], Memr[x+npts-1]) + call alimr (Memr[y], npts, ymin, ymax) + dy = ymax - ymin + call gswind (EC_GP(ec), xmin, xmax, ymin - .2 * dy, ymax + .2 * dy) + call gt_swind (EC_GP(ec), EC_GT(ec)) + call gt_labax (EC_GP(ec), EC_GT(ec)) + do i = 1, EC_NLINES(ec) { + call gt_plot (EC_GP(ec), EC_GT(ec), Memr[x], Memr[y], EC_NPTS(ec)) + x = x + EC_NPTS(ec) + y = y + EC_NPTS(ec) + } + + do i = 1, EC_NFEATURES(ec) + call ec_mark (ec, i) + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/ecidentify.h b/noao/onedspec/ecidentify/ecidentify.h new file mode 100644 index 00000000..63e4c6bd --- /dev/null +++ b/noao/onedspec/ecidentify/ecidentify.h @@ -0,0 +1,94 @@ +# Task parameters + +define LEN_EC 52 # Length ID structure + +define EC_IMAGE Memi[$1] # Image name (pointer) +define EC_MAXFEATURES Memi[$1+1] # Maximum number of features +define EC_FTYPE Memi[$1+2] # Feature type +define EC_MINSEP Memr[P2R($1+3)] # Minimum pixel separation +define EC_MATCH Memr[P2R($1+4)] # Maximum matching separation +define EC_FWIDTH Memr[P2R($1+5)] # Feature width in pixels +define EC_CRADIUS Memr[P2R($1+6)] # Centering radius in pixels +define EC_THRESHOLD Memr[P2R($1+7)] # Centering threshold +define EC_ZWIDTH Memr[P2R($1+8)] # Zoom window width in fit units +define EC_DATABASE Memi[$1+9] # Name of database (pointer) +define EC_COORDLIST Memi[$1+10] # Name of coordinate list (pointer) +define EC_LABELS Memi[$1+11] # Type of feature labels +define EC_LOGFILES Memi[$1+12] # List of logfiles + +# Common image data + +define EC_NCOLS Memi[$1+13] # Number of columns +define EC_NLINES Memi[$1+14] # Number of lines/apertures/orders +define EC_SHS Memi[$1+15] # Pointer to SHDR pointers +define EC_PIXDATA Memi[$1+16] # Pixel coordinates (pointer) +define EC_FITDATA Memi[$1+17] # Fit coordinates (pointer) + +define EC_IMLINE Memi[$1+18] # Image data (pointer) +define EC_PIXLINE Memi[$1+19] # Pixel coordinates (pointer) +define EC_FITLINE Memi[$1+20] # Fit coordinates (pointer) +define EC_NPTS Memi[$1+21] # Number of points + +define EC_SHIFT Memd[P2D($1+22)]# Wavelength shift + +# Features + +define EC_NFEATURES Memi[$1+24] # Number of features +define EC_NALLOC Memi[$1+25] # Length of allocated feature arrays +define EC_APNUM Memi[$1+26] # Aperture number (pointer) +define EC_LINENUM Memi[$1+27] # Image line number (pointer) +define EC_ORD Memi[$1+28] # Feature order number (pointer) +define EC_PIX Memi[$1+29] # Feature pixel coordinates (pointer) +define EC_FIT Memi[$1+30] # Feature fit coordinates (pointer) +define EC_USER Memi[$1+31] # Feature user coordinates (pointer) +define EC_FWIDTHS Memi[$1+32] # Feature width (pointer) +define EC_FTYPES Memi[$1+33] # Feature type (pointer) + +# Current status + +define EC_CURRENT Memi[$1+34] # Current feature +define EC_SH Memi[$1+35] # Current SHDR pointer +define EC_AP Memi[$1+36] # Current aperture +define EC_LINE Memi[$1+37] # Current line +define EC_ORDER Memi[$1+38] # Current order + +# Pointers for other packages + +define EC_LP Memi[$1+39] # Logical to physical transformation +define EC_PL Memi[$1+40] # Physical to logical transformation +define EC_LL Memi[$1+41] # Linelist pointer +define EC_ECF Memi[$1+42] # Curfit pointer +define EC_GP Memi[$1+43] # GIO pointer +define EC_GT Memi[$1+44] # Gtools pointer +define EC_UN Memi[$1+45] # Units pointer + +# Flags + +define EC_NEWFEATURES Memi[$1+46] # Has feature list changed? +define EC_NEWECF Memi[$1+47] # Has fitting function changed? +define EC_NEWGRAPH Memi[$1+48] # Has graph changed? +define EC_NEWDBENTRY Memi[$1+49] # Has database entry changed? +define EC_REFIT Memi[$1+50] # Refit feature data? +define EC_GTYPE Memi[$1+51] # Graph type + +# End of structure ---------------------------------------------------------- + +define LABELS "|none|index|pixel|user|" +define FTYPES "|emission|absorption|" + +define IMDATA Memr[EC_IMLINE($1)+$2-1] +define PIXDATA Memd[EC_PIXLINE($1)+$2-1] +define FITDATA Memd[EC_FITLINE($1)+$2-1] + +define SH Memi[EC_SHS($1)+$2-1] +define APS AP(SH($1,$2)) +define ORDERS BEAM(SH($1,$2)) + +define APN Memi[EC_APNUM($1)+$2-1] +define LINE Memi[EC_LINENUM($1)+$2-1] +define ORDER Memi[EC_ORD($1)+$2-1] +define PIX Memd[EC_PIX($1)+$2-1] +define FIT Memd[EC_FIT($1)+$2-1] +define USER Memd[EC_USER($1)+$2-1] +define FWIDTH Memr[EC_FWIDTHS($1)+$2-1] +define FTYPE Memi[EC_FTYPES($1)+$2-1] diff --git a/noao/onedspec/ecidentify/ecidentify.key b/noao/onedspec/ecidentify/ecidentify.key new file mode 100644 index 00000000..c19698ef --- /dev/null +++ b/noao/onedspec/ecidentify/ecidentify.key @@ -0,0 +1,76 @@ +1. ECIDENTIFY CURSOR KEY SUMMARY + +? Help a Affect all features c Center feature(s) +d Delete feature(s) f Fit dispersion g Fit zero point shift +i Initialize j Go to previous order k Go to next order +l Match coordinate list m Mark feature n Next feature +o Go to specified order p Pan graph q Quit +r Redraw graph s Shift feature t Reset position +u Enter user coordinate w Window graph x Crosscorrelate peaks +y Find peaks z Zoom graph . Nearest feature ++ Next feature - Previous feature I Interrupt + + +2. ECIDENTIFY 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] :write [image] +:zwidth [value] + + +3. ECIDENTIFY CURSOR KEYS + +? Clear the screen and print menu of options +a Apply next (c)enter or (d)elete operation to (a)ll features +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 +i (I)nitialize (delete features and coordinate fit) +j Go to the previous order +k Go to the next order +l Match coordinates in the coordinate (l)ist to features in the data +m (M)ark a new feature near the cursor +n Move the cursor or zoom to the (n)ext feature (same as +) +o Go to the specified (o)rder +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 (move) the position of a feature without centering +u Enter a new (u)ser coordinate for the current feature +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 (also space bar) ++ Move the cursor or zoom to the next feature +- Move the cursor or zoom to the previous feature +I Interrupt task immediately. Database information is not saved. + + +4. ECIDENTIFY 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 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, or user) +:match value Coordinate list matching distance +:maxfeatures value Maximum number of features automatically found +:minsep value Minimum separation allowed between features +:read name Read a record from the database (name defaults to image) +:write name Write a record to the database (name defaults to image) +:zwidth value Zoom width in user units diff --git a/noao/onedspec/ecidentify/ecidentify.x b/noao/onedspec/ecidentify/ecidentify.x new file mode 100644 index 00000000..827568d1 --- /dev/null +++ b/noao/onedspec/ecidentify/ecidentify.x @@ -0,0 +1,535 @@ +include <error.h> +include <imhdr.h> +include <gset.h> +include <smw.h> +include "ecidentify.h" + +define HELP "noao$onedspec/ecidentify/ecidentify.key" +define PROMPT "ecidentify options" + +define PAN 1 # Pan graph +define ZOOM 2 # Zoom graph + +# EC_IDENTIFY -- Identify echelle features in an image. +# This is the basic interactive loop. + +procedure ec_identify (ec) + +pointer ec # EC pointer + +real wx, wy +int wcs, key +char cmd[SZ_LINE] + +char newimage[SZ_FNAME] +int i, j, last, all, prfeature, nfeatures1, npeaks +bool answer +double pix, fit, user, shift, pix_shift, z_shift +pointer peaks + +bool clgetb() +int clgcur(), scan(), nscan(), find_peaks(), ec_next(), ec_previous() +int ec_line() +double ec_center(), ec_fittopix(), ec_fitpt(), ec_shift(), ec_rms() +double ecf_getd() +errchk ec_gdata(), ec_graph(), ec_dbread(), xt_mk1d(), ec_line() + +define newim_ 10 +define newkey_ 20 +define beep_ 99 + +begin +newim_ # Start here for each new image. + + # Get the image data. Return if there is an error. + iferr (call ec_gdata (ec)) { + call erract (EA_WARN) + return + } + + # Look for a database entry for the image. + iferr { + call ec_dbread (ec, Memc[EC_IMAGE(ec)], NO) + EC_NEWDBENTRY(ec) = NO + } then + if ((EC_NFEATURES(ec) > 0) || (EC_ECF(ec) != NULL)) + EC_NEWDBENTRY(ec) = YES + + # Set the coordinate array and the feature data. + iferr (call ec_fitdata (ec)) + call erract (EA_WARN) + call ec_fitfeatures (ec) + + # Begin with the first image line. + EC_LINE(ec) = 1 + EC_AP(ec) = APS(ec,EC_LINE(ec)) + EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec)) + call ec_gline (ec, EC_LINE(ec)) + + # Initialize. + EC_GTYPE(ec) = PAN + EC_REFIT(ec) = NO + EC_NEWFEATURES(ec) = NO + EC_NEWECF(ec) = NO + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + last = EC_CURRENT(ec) + all = 0 + newimage[1] = EOS + key = 'r' + + repeat { + prfeature = YES + if (all != 0) + all = mod (all + 1, 3) + + switch (key) { + case '?': # Page help + call gpagefile (EC_GP(ec), HELP, PROMPT) + case ':': # Execute colon commands + if (cmd[1] == '/') + call gt_colon (cmd, EC_GP(ec), EC_GT(ec), EC_NEWGRAPH(ec)) + else + call ec_colon (ec, cmd, newimage, prfeature) + case ' ': # Go to the current feature + case '.': # Go to the nearest feature + if (EC_NFEATURES(ec) == 0) + goto beep_ + call ec_nearest (ec, double (wx)) + case '-': # Go to the previous feature + if (ec_previous (ec, EC_CURRENT(ec)) == EOF) + goto beep_ + case '+', 'n': # Go to the next feature + if (ec_next (ec, EC_CURRENT(ec)) == EOF) + goto beep_ + case 'a': # Set the all flag for the next key + all = 1 + case 'c': # Center features on data + if (all != 0) { + call eprintf ("Recentering features ...\n") + for (i = 1; i <= EC_NFEATURES(ec); i = i + 1) { + call ec_gline (ec, LINE(ec,i)) + call gseti (EC_GP(ec), G_PLTYPE, 0) + call ec_mark (ec, i) + call gseti (EC_GP(ec), G_PLTYPE, 1) + FWIDTH(ec,i) = EC_FWIDTH(ec) + PIX(ec,i) = ec_center (ec, PIX(ec,i), FWIDTH(ec,i), + FTYPE(ec,i)) + if (!IS_INDEFD (PIX(ec,i))) { + FIT(ec,i) = ec_fitpt (ec, APN(ec,i), PIX(ec,i)) + call ec_mark (ec, i) + } else { + call ec_delete (ec, i) + i = i - 1 + } + } + call ec_gline (ec, EC_LINE(ec)) + EC_NEWFEATURES(ec) = YES + } else { + if (EC_NFEATURES(ec) == 0) + goto beep_ + + call ec_nearest (ec, double (wx)) + pix = PIX(ec,EC_CURRENT(ec)) + pix = ec_center (ec, pix, EC_FWIDTH(ec), + FTYPE(ec,EC_CURRENT(ec))) + if (!IS_INDEFD (pix)) { + call gseti (EC_GP(ec), G_PLTYPE, 0) + call ec_mark (ec, EC_CURRENT(ec)) + PIX(ec,EC_CURRENT(ec)) = pix + FWIDTH(ec,EC_CURRENT(ec)) = EC_FWIDTH(ec) + FIT(ec,EC_CURRENT(ec)) = + ec_fitpt (ec, APN(ec,EC_CURRENT(ec)), pix) + call gseti (EC_GP(ec), G_PLTYPE, 1) + call ec_mark (ec, EC_CURRENT(ec)) + EC_NEWFEATURES(ec) = YES + } else { + call eprintf ("Centering failed\n") + prfeature = NO + } + } + case 'd': # Delete features + if (all != 0) { + EC_NFEATURES(ec) = 0 + EC_CURRENT(ec) = 0 + EC_NEWFEATURES(ec) = YES + EC_NEWGRAPH(ec) = YES + } else { + if (EC_NFEATURES(ec) == 0) + goto beep_ + + call ec_nearest (ec, double (wx)) + call gseti (EC_GP(ec), G_PLTYPE, 0) + call ec_mark (ec, EC_CURRENT(ec)) + call gseti (EC_GP(ec), G_PLTYPE, 1) + call ec_delete (ec, EC_CURRENT(ec)) + call ec_nearest (ec, double (wx)) + last = 0 + } + case 'f': # Fit dispersion function + iferr (call ec_dofit (ec, YES, NO)) { + call erract (EA_WARN) + prfeature = NO + goto beep_ + } + case 'g': # Fit shift + call ec_doshift (ec, YES) + prfeature = NO + case 'i': # Initialize + call dgsfree (EC_ECF(ec)) + call ecf_setd ("shift", 0.D0) + EC_NEWECF(ec) = YES + EC_NFEATURES(ec) = 0 + EC_CURRENT(ec) = 0 + EC_NEWFEATURES(ec) = YES + EC_NEWGRAPH(ec) = YES + case 'j': # Go to the previous order + EC_LINE(ec) = + mod (EC_LINE(ec)+EC_NLINES(ec)-2, EC_NLINES(ec)) + 1 + EC_AP(ec) = APS(ec,EC_LINE(ec)) + EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec)) + call ec_gline (ec, EC_LINE(ec)) + EC_NEWGRAPH(ec) = YES + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + case 'k': # Go to the next order + EC_LINE(ec) = mod (EC_LINE(ec), EC_NLINES(ec)) + 1 + EC_AP(ec) = APS(ec,EC_LINE(ec)) + EC_ORDER(ec) = ORDERS(ec,EC_LINE(ec)) + call ec_gline (ec, EC_LINE(ec)) + EC_NEWGRAPH(ec) = YES + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + case 'l': # Find features using a line list + if (EC_ECF(ec) == NULL) { + call eprintf ("Doing initial fit ...\n") + iferr (call ec_dofit (ec, NO, NO)) { + call erract (EA_WARN) + prfeature = NO + goto beep_ + } + if (EC_NEWECF(ec) == YES) { + iferr (call ec_fitdata (ec)) { + call erract (EA_WARN) + prfeature = NO + } + call ec_fitfeatures (ec) + EC_NEWECF(ec) = NO + } + } + + call eprintf ("Searching coordinate list ...\n") + call ec_linelist (ec) + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + if (EC_NEWFEATURES(ec) == YES) + EC_NEWGRAPH(ec) = YES + case 'm': # Mark a new feature + fit = wx + pix = ec_fittopix (ec, fit) + pix = ec_center (ec, pix, EC_FWIDTH(ec), EC_FTYPE(ec)) + if (IS_INDEFD (pix)) + goto beep_ + fit = ec_fitpt (ec, EC_AP(ec), pix) + user = fit + call ec_newfeature (ec, EC_AP(ec), pix, fit, user, + EC_FWIDTH(ec), EC_FTYPE(ec)) + USER(ec,EC_CURRENT(ec)) = INDEFD + call ec_match (ec, FIT(ec,EC_CURRENT(ec)), + USER(ec,EC_CURRENT(ec))) + call ec_mark (ec, EC_CURRENT(ec)) + call printf ("%3d %10.2f %10.8g (%10.8g): ") + call pargi (APN(ec,EC_CURRENT(ec))) + call pargd (PIX(ec,EC_CURRENT(ec))) + call pargd (FIT(ec,EC_CURRENT(ec))) + call pargd (USER(ec,EC_CURRENT(ec))) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() == 1) { + USER(ec,EC_CURRENT(ec)) = user + call ec_match (ec, user, USER(ec,EC_CURRENT(ec))) + } + } + case 'o': # Go to a specified order + call printf ("Aperture (%d): ") + call pargi (EC_AP(ec)) + call flush (STDOUT) + if (scan() != EOF) { + call gargi (j) + if (nscan() == 1) { + if (j != EC_AP(ec)) { + iferr { + i = ec_line (ec, j) + call ec_gline (ec, i) + EC_LINE(ec) = i + EC_AP(ec) = j + EC_ORDER(ec) = ORDERS(ec,i) + EC_NEWGRAPH(ec) = YES + EC_CURRENT(ec) = 0 + i = ec_next (ec, EC_CURRENT(ec)) + } then + goto beep_ + } + } + } + case 'p': # Go to pan graph mode + if (EC_GTYPE(ec) == PAN) + goto beep_ + + EC_GTYPE(ec) = PAN + EC_NEWGRAPH(ec) = YES + case 'q': # Quit + break + case 'r': # Redraw the current graph + EC_NEWGRAPH(ec) = YES + case 's', 'x': # Shift or cross correlate features + # 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) * EC_ORDER(ec) + } else + shift = 0. + case 'x': + if (EC_NFEATURES(ec) > 5) { + call eprintf ( + "Cross correlating features with peaks ...\n") + shift = ec_shift (ec) + } else + goto beep_ + } + + EC_NEWFEATURES(ec) = YES + EC_NEWECF(ec) = YES + EC_NEWGRAPH(ec) = YES + prfeature = NO + + if (EC_NFEATURES(ec) < 1) { + call printf ("User coordinate shift=%5f") + call pargd (shift / EC_ORDER(ec)) + call ecf_setd ("shift", ecf_getd ("shift") - shift) + goto newkey_ + } + + # Recenter features. + call eprintf ("Recentering features ...\n") + pix_shift = 0. + z_shift = 0. + nfeatures1 = EC_NFEATURES(ec) + + j = 0. + do i = 1, EC_NFEATURES(ec) { + call ec_gline (ec, LINE(ec,i)) + pix = ec_fittopix (ec, FIT(ec,i) + shift/ORDER(ec,i)) + pix = ec_center (ec, pix, FWIDTH(ec,i), FTYPE(ec,i)) + if (IS_INDEFD (pix)) { + if (EC_CURRENT(ec) == i) + EC_CURRENT(ec) = 0 + next + } + fit = ec_fitpt (ec, APN(ec,i), pix) + + pix_shift = pix_shift + pix - PIX(ec,i) + if (FIT(ec,i) != 0.) + z_shift = z_shift + (fit - FIT(ec,i)) / FIT(ec,i) + + j = j + 1 + APN(ec,j) = APN(ec,i) + LINE(ec,j) = LINE(ec,i) + ORDER(ec,j) = ORDER(ec,i) + PIX(ec,j) = pix + FIT(ec,j) = FIT(ec,i) + USER(ec,j) = USER(ec,i) + FWIDTH(ec,j) = FWIDTH(ec,i) + FTYPE(ec,j) = FTYPE(ec,i) + if (EC_CURRENT(ec) == i) + EC_CURRENT(ec) = j + } + call ec_gline (ec, EC_LINE(ec)) + EC_NFEATURES(ec) = j + if (EC_CURRENT(ec) == 0) + i = ec_next (ec, EC_CURRENT(ec)) + + if (EC_NFEATURES(ec) < 1) { + call printf ("User coordinate shift=%5f") + call pargd (shift / EC_ORDER(ec)) + call printf (", No features found during recentering") + call ecf_setd ("shift", ecf_getd ("shift") - shift) + goto newkey_ + } + + # Adjust shift. + pix = ecf_getd ("shift") + call ec_doshift (ec, NO) + call ec_fitfeatures (ec) + + # Print results. + call printf ("Recentered=%d/%d") + call pargi (EC_NFEATURES(ec)) + call pargi (nfeatures1) + call printf ( + ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g") + call pargd (pix_shift / EC_NFEATURES(ec)) + call pargd ((pix - ecf_getd ("shift")) / EC_ORDER(ec)) + call pargd (z_shift / EC_NFEATURES(ec)) + call pargd (ec_rms(ec)) + case 't': # Move current feature + if (EC_CURRENT(ec) == 0) + goto beep_ + + call gseti (EC_GP(ec), G_PLTYPE, 0) + call ec_mark (ec, EC_CURRENT(ec)) + pix = ec_fittopix (ec, double (wx)) + PIX(ec,EC_CURRENT(ec)) = pix + FIT(ec,EC_CURRENT(ec)) = + ec_fitpt (ec, APN(ec,EC_CURRENT(ec)), pix) + call gseti (EC_GP(ec), G_PLTYPE, 1) + call ec_mark (ec, EC_CURRENT(ec)) + EC_NEWFEATURES(ec) = YES + case 'u': # Set uesr coordinate value + if (EC_NFEATURES(ec) == 0) + goto beep_ + + call printf ("%3d %10.2f %10.8g (%10.8g): ") + call pargi (APN(ec,EC_CURRENT(ec))) + call pargd (PIX(ec,EC_CURRENT(ec))) + call pargd (FIT(ec,EC_CURRENT(ec))) + call pargd (USER(ec,EC_CURRENT(ec))) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() == 1) { + USER(ec,EC_CURRENT(ec)) = user + EC_NEWFEATURES(ec) = YES + } + } + case 'w': # Window graph + call gt_window (EC_GT(ec), EC_GP(ec), "cursor", EC_NEWGRAPH(ec)) + case 'y': # Find peaks in order + call malloc (peaks, EC_NPTS(ec), TY_REAL) + npeaks = find_peaks (IMDATA(ec,1), Memr[peaks], + EC_NPTS(ec), 0., int (EC_MINSEP(ec)), 0, EC_MAXFEATURES(ec), + 0., false) + for (j = 1; j <= EC_NFEATURES(ec); j = j + 1) { + for (i = 1; i <= npeaks; i = i + 1) { + if (!IS_INDEF(pix)) { + pix = Memr[peaks + i - 1] + if (abs (pix - PIX(ec,j)) < EC_MINSEP(ec)) + Memr[peaks + i - 1] = INDEF + } + } + } + for (i = 1; i <= npeaks; i = i + 1) { + pix = Memr[peaks+i-1] + pix = ec_center (ec, pix, EC_FWIDTH(ec), EC_FTYPE(ec)) + if (IS_INDEFD (pix)) + next + fit = ec_fitpt (ec, EC_AP(ec), pix) + user = INDEFD + call ec_match (ec, fit, user) + call ec_newfeature (ec, EC_AP(ec), pix, fit, user, + EC_FWIDTH(ec), EC_FTYPE(ec)) + call ec_mark (ec, EC_CURRENT(ec)) + } + call mfree (peaks, TY_REAL) + case 'z': # Go to zoom mode + if (EC_CURRENT(ec) == 0) + goto beep_ + + if (EC_GTYPE(ec) == PAN) + EC_NEWGRAPH(ec) = YES + EC_GTYPE(ec) = ZOOM + call ec_nearest (ec, double (wx)) + case 'I': # Interrupt + call fatal (0, "Interrupt") + default: # Beep +beep_ call printf ("\007\n") + } + +newkey_ + # Set database update flag if there has been a change. + if ((EC_NEWFEATURES(ec) == YES) || (EC_NEWECF(ec) == YES)) + EC_NEWDBENTRY(ec) = YES + + # Exit loop and then start new image. + if (newimage[1] != EOS) + break + + # Refit the dispersion function if needed. + if (EC_REFIT(ec) == YES) { + iferr (call ec_dofit (ec, NO, NO)) { + call erract (EA_WARN) + prfeature = NO + } + EC_REFIT(ec) = NO + } + + # Recompute the coordinate information. + if (EC_NEWECF(ec) == YES) { + iferr (call ec_fitdata (ec)) { + call erract (EA_WARN) + prfeature = NO + } + call ec_fitfeatures (ec) + EC_NEWECF(ec) = NO + } + + # Redraw new feature in zoom mode. + if ((EC_GTYPE(ec) == ZOOM) && (last != EC_CURRENT(ec))) + EC_NEWGRAPH(ec) = YES + + # Redraw graph. + if (EC_NEWGRAPH(ec) == YES) { + call ec_graph (ec, EC_GTYPE(ec)) + EC_NEWGRAPH(ec) = NO + } + + # Set cursor and print current feature on status (unless canceled). + if (EC_CURRENT(ec) > 0) { + call gscur (EC_GP(ec), real (FIT(ec,EC_CURRENT(ec))), wy) + if (prfeature == YES) { + call printf ("%d %10.2f %10.8g %10.8g\n") + call pargi (APN(ec,EC_CURRENT(ec))) + call pargd (PIX(ec,EC_CURRENT(ec))) + call pargd (FIT(ec,EC_CURRENT(ec))) + call pargd (USER(ec,EC_CURRENT(ec))) + } + } + + last = EC_CURRENT(ec) + } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + + # Warn user that feature data is newer than database entry. + if (EC_NEWDBENTRY(ec) == YES) { + answer = true + if (!clgetb ("autowrite")) { + call printf ("Write feature data to the database (yes)? ") + call flush (STDOUT) + if (scan() != EOF) + call gargb (answer) + } + if (answer) + call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], NO) + } + + call flush (STDOUT) + + # Free image data and MWCS + call mfree (EC_PIXDATA(ec), TY_DOUBLE) + call mfree (EC_FITDATA(ec), TY_DOUBLE) + call smw_close (MW(EC_SH(ec))) + do i = 1, EC_NLINES(ec) + MW(SH(ec,i)) = NULL + + # If a new image was specified by a colon command don't return. + if (newimage[1] != EOS) { + call strcpy (newimage, Memc[EC_IMAGE(ec)], SZ_FNAME) + goto newim_ + } +end diff --git a/noao/onedspec/ecidentify/ecinit.x b/noao/onedspec/ecidentify/ecinit.x new file mode 100644 index 00000000..8b3b7b62 --- /dev/null +++ b/noao/onedspec/ecidentify/ecinit.x @@ -0,0 +1,64 @@ +include <gset.h> +include "ecidentify.h" + +# EC_INIT -- Allocate and initialize the identify structure. + +procedure ec_init (ec) + +pointer ec # ID pointer + +begin + call calloc (ec, LEN_EC, TY_STRUCT) + + EC_NALLOC(ec) = 20 + EC_NFEATURES(ec) = 0 + EC_CURRENT(ec) = 0 + EC_NLINES(ec) = 0 + EC_LL(ec) = NULL + EC_ECF(ec) = NULL + EC_LABELS(ec) = 1 + + call malloc (EC_IMAGE(ec), SZ_FNAME, TY_CHAR) + call malloc (EC_DATABASE(ec), SZ_FNAME, TY_CHAR) + call malloc (EC_COORDLIST(ec), SZ_FNAME, TY_CHAR) + + call malloc (EC_APNUM(ec), EC_NALLOC(ec), TY_INT) + call malloc (EC_LINENUM(ec), EC_NALLOC(ec), TY_INT) + call malloc (EC_PIX(ec), EC_NALLOC(ec), TY_DOUBLE) + call malloc (EC_ORD(ec), EC_NALLOC(ec), TY_INT) + call malloc (EC_FIT(ec), EC_NALLOC(ec), TY_DOUBLE) + call malloc (EC_USER(ec), EC_NALLOC(ec), TY_DOUBLE) + call malloc (EC_FWIDTHS(ec), EC_NALLOC(ec), TY_REAL) + call malloc (EC_FTYPES(ec), EC_NALLOC(ec), TY_INT) +end + + +# EC_FREE -- Free identify structure. + +procedure ec_free (ec) + +pointer ec # ID pointer +int i + +begin + if (EC_UN(ec) != NULL) + call un_close (EC_UN(ec)) + do i = 1, EC_NLINES(ec) + call shdr_close (SH(ec,i)) + call mfree (EC_SHS(ec), TY_POINTER) + + call mfree (EC_IMAGE(ec), TY_CHAR) + call mfree (EC_DATABASE(ec), TY_CHAR) + call mfree (EC_COORDLIST(ec), TY_CHAR) + + call mfree (EC_APNUM(ec), TY_INT) + call mfree (EC_LINENUM(ec), TY_INT) + call mfree (EC_PIX(ec), TY_DOUBLE) + call mfree (EC_ORD(ec), TY_INT) + call mfree (EC_FIT(ec), TY_DOUBLE) + call mfree (EC_USER(ec), TY_DOUBLE) + call mfree (EC_FWIDTHS(ec), TY_REAL) + call mfree (EC_FTYPES(ec), TY_INT) + + call mfree (ec, TY_STRUCT) +end diff --git a/noao/onedspec/ecidentify/ecline.x b/noao/onedspec/ecidentify/ecline.x new file mode 100644 index 00000000..63e55072 --- /dev/null +++ b/noao/onedspec/ecidentify/ecline.x @@ -0,0 +1,22 @@ +include <smw.h> +include "ecidentify.h" + +# EC_LINE -- Get line corresponding to aperture. + +int procedure ec_line (ec, ap) + +pointer ec # EC pointer +int ap # Aperture + +int i + +begin + if (IS_INDEFI (ap)) + return (INDEFI) + + do i = 1, EC_NLINES(ec) + if (ap == APS(ec,i)) + return (i) + + call error (0, "Image line for aperture number not found") +end diff --git a/noao/onedspec/ecidentify/eclinelist.x b/noao/onedspec/ecidentify/eclinelist.x new file mode 100644 index 00000000..6653dd4b --- /dev/null +++ b/noao/onedspec/ecidentify/eclinelist.x @@ -0,0 +1,281 @@ +include <error.h> +include <mach.h> +include <smw.h> +include <units.h> +include "ecidentify.h" + +# EC_MAPLL -- Read the line list into memory. + +procedure ec_mapll (ec) + +pointer ec # Echelle pointer + +int fd, nalloc, nlines, open(), fscan(), nscan() +double value, lastval +pointer ec_ll +pointer sp, str, units, un_open() +bool streq() +errchk open, fscan, malloc, realloc, un_open + +begin + EC_LL(ec) = NULL + + call xt_stripwhite (Memc[EC_COORDLIST(ec)]) + if (Memc[EC_COORDLIST(ec)] == EOS) + return + iferr (fd = open (Memc[EC_COORDLIST(ec)], READ_ONLY, TEXT_FILE)) + return + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (units, SZ_LINE, TY_CHAR) + call strcpy ("Angstroms", Memc[units], SZ_LINE) + + lastval = -MAX_DOUBLE + nalloc = 0 + nlines = 0 + while (fscan (fd) != EOF) { + call gargwrd (Memc[str], SZ_LINE) + if (nscan() != 1) + next + if (Memc[str] == '#') { + call gargwrd (Memc[str], SZ_LINE) + call strlwr (Memc[str]) + if (streq (Memc[str], "units")) { + call gargstr (Memc[units], SZ_LINE) + call xt_stripwhite (Memc[units]) + } + next + } + call reset_scan () + + call gargd (value) + if (nscan() != 1) + next + + if (nalloc == 0) { + nalloc = 100 + call malloc (ec_ll, nalloc, TY_DOUBLE) + } else if (nlines == nalloc) { + nalloc = nalloc + 100 + call realloc (ec_ll, nalloc, TY_DOUBLE) + } + + if (value < lastval) { + call close (fd) + call mfree (ec_ll, TY_DOUBLE) + call error (0, "Line list not sorted in increasing order") + } + + Memd[ec_ll+nlines] = value + nlines = nlines + 1 + } + call close (fd) + + if (nlines > 0) { + call realloc (ec_ll, nlines + 1, TY_DOUBLE) + Memd[ec_ll+nlines] = INDEFD + EC_LL(ec) = ec_ll + + if (EC_UN(ec) == NULL && Memc[units] != EOS) + EC_UN(ec) = un_open (Memc[units]) + call ec_unitsll (ec, Memc[units]) + } + + call sfree (sp) +end + + +# EC_UNMAPLL -- Unmap the linelist. + +procedure ec_unmapll (ec) + +pointer ec # Line list pointer + +begin + call mfree (EC_LL(ec), TY_DOUBLE) +end + + +# EC_UNITSLL -- Change the line list units from the input units to the +# units given by EC_UN. This may involve reversing the order of the list. + +procedure ec_unitsll (ec, units) + +pointer ec # Identify structure +char units[ARB] # Input units + +int i, nll +double value +pointer un, ll, llend, un_open() +bool un_compare() +errchk un_open + +begin + if (EC_LL(ec) == NULL) + return + if (IS_INDEFD(Memd[EC_LL(ec)])) + return + if (units[1] == EOS || EC_UN(ec) == NULL) + return + if (UN_CLASS(EC_UN(ec)) == UN_UNKNOWN) + return + + un = un_open (units) + if (un_compare (un, EC_UN(ec))) { + call un_close (un) + return + } + + ll = EC_LL(ec) + do i = 0, ARB + if (IS_INDEFD(Memd[ll+i])) { + nll = i + break + } + call un_ctrand (un, EC_UN(ec), Memd[ll], Memd[ll], nll) + call un_close (un) + + if (Memd[ll] > Memd[ll+nll-1]) { + llend = ll + nll - 1 + do i = 0, nll / 2 - 1 { + value = Memd[ll+i] + Memd[ll+i] = Memd[llend-i] + Memd[llend-i] = value + } + } +end + + + +# EC_MATCH -- Match current feature against a line list. + +procedure ec_match (ec, in, out) + +pointer ec # Echelle pointer +double in # Coordinate to be matched +double out # Matched coordinate + +double match, alpha, delta, delta1, delta2, out1 +pointer ll + +begin + if (EC_LL(ec) == NULL) { + out = in + return + } + + match = EC_MATCH(ec) + alpha = 1.25 + delta1 = MAX_REAL + + # Find nearest match. + for (ll=EC_LL(ec); !IS_INDEFD(Memd[ll]); ll = ll + 1) { + delta = abs (in - Memd[ll]) + if (delta < delta1) { + delta2 = delta1 + delta1 = delta + if (delta1 <= match) + out1 = Memd[ll] + } + } + + # Only return match if no other candidate is also possible. + if (delta1 > match) + return + if (delta2 < alpha * delta1) + return + + out = out1 +end + +# EC_LINELIST -- Add features from a line list. + +procedure ec_linelist (ec) + +pointer ec # Echelle pointer + +int i, line, ap, nfound, nextpix +double pix, fit, user, peak, minval, match, fit1, fit2 +pointer sp, aps, pixes, fits, users, peaks, ll + +double ec_center(), ec_fittopix(), ec_fitpt(), ec_peak() + +begin + if (EC_LL(ec) == NULL) + return + + call smark (sp) + call salloc (aps, EC_MAXFEATURES(ec), TY_INT) + call salloc (pixes, EC_MAXFEATURES(ec), TY_DOUBLE) + call salloc (fits, EC_MAXFEATURES(ec), TY_DOUBLE) + call salloc (users, EC_MAXFEATURES(ec), TY_DOUBLE) + call salloc (peaks, EC_MAXFEATURES(ec), TY_DOUBLE) + + nfound = 0 + minval = MAX_REAL + + do line = 1, EC_NLINES(ec) { + call ec_gline (ec, line) + ap = APS(ec,line) + fit1 = min (FITDATA(ec,1), FITDATA(ec,EC_NPTS(ec))) + fit2 = max (FITDATA(ec,1), FITDATA(ec,EC_NPTS(ec))) + for (ll=EC_LL(ec); !IS_INDEFD(Memd[ll]); ll = ll + 1) { + user = Memd[ll] + if (user < fit1) + next + if (user > fit2) + break + + pix = ec_center (ec, ec_fittopix (ec, user), EC_FWIDTH(ec), + EC_FTYPE(ec)) + if (!IS_INDEFD(pix)) { + fit = ec_fitpt (ec, ap, pix) + match = abs (fit - user) + if (match > EC_MATCH(ec)) + next + + peak = abs (ec_peak (ec, pix)) + if (nfound < EC_MAXFEATURES(ec)) { + nfound = nfound + 1 + if (peak < minval) { + nextpix = nfound + minval = peak + } + Memi[aps+nfound-1] = ap + Memd[pixes+nfound-1] = pix + Memd[fits+nfound-1] = fit + Memd[users+nfound-1] = user + Memd[peaks+nfound-1] = peak + } else if (peak > minval) { + Memi[aps+nextpix-1] = ap + Memd[pixes+nextpix-1] = pix + Memd[fits+nextpix-1] = fit + Memd[users+nextpix-1] = user + Memd[peaks+nextpix-1] = peak + + minval = MAX_REAL + do i = 1, nfound { + peak = Memd[peaks+i-1] + if (peak < minval) { + nextpix = i + minval = peak + } + } + } + } + } + } + call ec_gline (ec, EC_LINE(ec)) + + do i = 1, nfound { + ap = Memi[aps+i-1] + pix = Memd[pixes+i-1] + fit = Memd[fits+i-1] + user = Memd[users+i-1] + call ec_newfeature (ec, ap, pix, fit, user, EC_FWIDTH(ec), + EC_FTYPE(ec)) + } + + call sfree (sp) +end diff --git a/noao/onedspec/ecidentify/eclog.x b/noao/onedspec/ecidentify/eclog.x new file mode 100644 index 00000000..e2730ca0 --- /dev/null +++ b/noao/onedspec/ecidentify/eclog.x @@ -0,0 +1,77 @@ +include <time.h> +include "ecidentify.h" + +# EC_LOG -- Write log + +procedure ec_log (ec, file) + +pointer ec # ID pointer +char file[ARB] # Log file + +char str[SZ_TIME] +int i, fd, nrms +double resid, rms + +int open() +long clktime() +errchk open() + +begin + if (EC_NFEATURES(ec) == 0) + return + + fd = open (file, APPEND, TEXT_FILE) + + call cnvtime (clktime (0), str, SZ_TIME) + call fprintf (fd, "\n%s\n") + call pargstr (str) + call fprintf (fd, "Features identified in image %s.\n") + call pargstr (Memc[EC_IMAGE(ec)]) + + call fprintf (fd, " %3s %4s %5s %8s %10s %10s %10s %6s %6d\n") + call pargstr ("Ap") + call pargstr ("Line") + call pargstr ("Order") + call pargstr ("Pixel") + call pargstr ("Fit") + call pargstr ("User") + call pargstr ("Residual") + call pargstr ("Fwidth") + call pargstr ("Reject") + + rms = 0. + nrms = 0 + do i = 1, EC_NFEATURES(ec) { + call fprintf (fd, + "%5d %3d %4d %5d %8.2f %10.8g %10.8g %10.8g %6.2f %6b\n") + call pargi (i) + call pargi (APN(ec,i)) + call pargi (LINE(ec,i)) + call pargi (ORDER(ec,i)) + call pargd (PIX(ec,i)) + call pargd (FIT(ec,i)) + call pargd (USER(ec,i)) + if (IS_INDEFD (USER(ec,i))) + call pargd (USER(ec,i)) + else { + resid = FIT(ec,i) - USER(ec,i) + call pargd (resid) + if (FTYPE(ec,i) > 0) { + rms = rms + resid ** 2 + nrms = nrms + 1 + } + } + call pargr (FWIDTH(ec,i)) + if (FTYPE(ec,i) > 0) + call pargb (false) + else + call pargb (true) + } + + if (nrms > 1) { + call fprintf (fd, "RMS = %0.8g\n") + call pargd (sqrt (rms / nrms)) + } + + call close (fd) +end diff --git a/noao/onedspec/ecidentify/ecmark.x b/noao/onedspec/ecidentify/ecmark.x new file mode 100644 index 00000000..58b02d0f --- /dev/null +++ b/noao/onedspec/ecidentify/ecmark.x @@ -0,0 +1,71 @@ +include <gset.h> +include <pkg/center1d.h> +include "ecidentify.h" + +procedure ec_mark (ec, feature) + +pointer ec # ID pointer +int feature + +int pix +real x, y +real mx, my, x1, x2, y1, y2, tick, gap +pointer sp, format, label +double smw_c1trand() + +define TICK .03 # Tick size in NDC +define GAP .02 # Gap size in NDC + +begin + call ggwind (EC_GP(ec), x1, x2, y1, y2) + + x = FIT(ec,feature) + + if ((x < min (x1, x2)) || (x > max (x1, x2))) + return + + pix = smw_c1trand (EC_PL(ec), PIX(ec,feature)) + pix = max (1, min (pix, EC_NPTS(ec) - 1)) + + call smark (sp) + call salloc (format, SZ_LINE, TY_CHAR) + call salloc (label, SZ_LINE, TY_CHAR) + switch (EC_FTYPE(ec)) { + case EMISSION: + y = max (IMDATA(ec,pix), IMDATA(ec,pix+1)) + tick = TICK + gap = GAP + call strcpy ("u=180;h=c;v=b;s=0.5", Memc[format], SZ_LINE) + case ABSORPTION: + y = min (IMDATA(ec,pix), IMDATA(ec,pix+1)) + tick = -TICK + gap = -GAP + call strcpy ("u=0;h=c;v=t;s=0.5", Memc[format], SZ_LINE) + } + + call gctran (EC_GP(ec), x, y, mx, my, 1, 0) + call gctran (EC_GP(ec), mx, my + gap, x1, y1, 0, 1) + call gctran (EC_GP(ec), mx, my + gap + tick, x1, y2, 0, 1) + call gline (EC_GP(ec), x1, y1, x1, y2) + + call gctran (EC_GP(ec), mx, my + tick + 2 * gap, x1, y2, 0, 1) + switch (EC_LABELS(ec)) { + case 2: + call sprintf (Memc[label], SZ_LINE, "%d") + call pargi (feature) + call gtext (EC_GP(ec), x1, y2, Memc[label], Memc[format]) + case 3: + call sprintf (Memc[label], SZ_LINE, "%0.2f") + call pargd (PIX(ec,feature)) + call gtext (EC_GP(ec), x1, y2, Memc[label], Memc[format]) + case 4: + if (!IS_INDEFD (USER(ec,feature))) { + call sprintf (Memc[label], SZ_LINE, "%0.4f") + call pargd (USER(ec,feature)) + call gtext (EC_GP(ec), x1, y2, Memc[label], Memc[format]) + } + } + + call sfree (sp) + call gflush (EC_GP(ec)) +end diff --git a/noao/onedspec/ecidentify/ecnearest.x b/noao/onedspec/ecidentify/ecnearest.x new file mode 100644 index 00000000..7b061472 --- /dev/null +++ b/noao/onedspec/ecidentify/ecnearest.x @@ -0,0 +1,26 @@ +include <mach.h> +include "ecidentify.h" + +# EC_NEAREST -- Find the nearest feature to a given coordinate. + +procedure ec_nearest (ec, fitnear) + +pointer ec # ID pointer +double fitnear # Coordinate to find nearest feature + +int i, ec_next() +double delta, delta1 + +begin + EC_CURRENT(ec) = 0 + + i = 0 + delta = MAX_REAL + while (ec_next (ec, i) != EOF) { + delta1 = abs (FIT(ec,i) - fitnear) + if (delta1 < delta) { + EC_CURRENT(ec) = i + delta = delta1 + } + } +end diff --git a/noao/onedspec/ecidentify/ecnewfeature.x b/noao/onedspec/ecidentify/ecnewfeature.x new file mode 100644 index 00000000..525c034a --- /dev/null +++ b/noao/onedspec/ecidentify/ecnewfeature.x @@ -0,0 +1,91 @@ +include <mach.h> +include <smw.h> +include "ecidentify.h" + +# EC_NEWFEATURE -- Allocate and initialize memory for a new feature. + +procedure ec_newfeature (ec, ap, pix, fit, user, width, type) + +pointer ec # ID pointer +int ap # Order +double pix # Pixel coordinate +double fit # Fit coordinate +double user # User coordinate +real width # Feature width +int type # Feature type + +int i, j, ec_line() +double delta + +define NALLOC 20 # Length of additional allocations + +begin + if (IS_INDEFD (pix)) + return + + delta = MAX_REAL + do i = 1, EC_NFEATURES(ec) { + if (APN(ec,i) != ap) + next + if (abs (pix - PIX(ec,i)) < delta) { + delta = abs (pix - PIX(ec,i)) + j = i + } + } + + if (delta >= EC_MINSEP(ec)) { + EC_NFEATURES(ec) = EC_NFEATURES(ec) + 1 + if (EC_NALLOC(ec) < EC_NFEATURES(ec)) { + EC_NALLOC(ec) = EC_NALLOC(ec) + NALLOC + call realloc (EC_APNUM(ec), EC_NALLOC(ec), TY_INT) + call realloc (EC_LINENUM(ec), EC_NALLOC(ec), TY_INT) + call realloc (EC_ORD(ec), EC_NALLOC(ec), TY_INT) + call realloc (EC_PIX(ec), EC_NALLOC(ec), TY_DOUBLE) + call realloc (EC_FIT(ec), EC_NALLOC(ec), TY_DOUBLE) + call realloc (EC_USER(ec), EC_NALLOC(ec), TY_DOUBLE) + call realloc (EC_FWIDTHS(ec), EC_NALLOC(ec), TY_REAL) + call realloc (EC_FTYPES(ec), EC_NALLOC(ec), TY_INT) + } + for (j=EC_NFEATURES(ec); (j>1)&&(ap<APN(ec,j-1)); j=j-1) { + APN(ec,j) = APN(ec,j-1) + LINE(ec,j) = LINE(ec,j-1) + ORDER(ec,j) = ORDER(ec,j-1) + PIX(ec,j) = PIX(ec,j-1) + FIT(ec,j) = FIT(ec,j-1) + USER(ec,j) = USER(ec,j-1) + FWIDTH(ec,j) = FWIDTH(ec,j-1) + FTYPE(ec,j) = FTYPE(ec,j-1) + } + for (; (j>1)&&(ap==APN(ec,j-1))&&(pix<PIX(ec,j-1)); j=j-1) { + APN(ec,j) = APN(ec,j-1) + LINE(ec,j) = LINE(ec,j-1) + ORDER(ec,j) = ORDER(ec,j-1) + PIX(ec,j) = PIX(ec,j-1) + FIT(ec,j) = FIT(ec,j-1) + USER(ec,j) = USER(ec,j-1) + FWIDTH(ec,j) = FWIDTH(ec,j-1) + FTYPE(ec,j) = FTYPE(ec,j-1) + } + APN(ec,j) = ap + LINE(ec,j) = ec_line (ec, ap) + ORDER(ec,j) = ORDERS(ec,LINE(ec,j)) + PIX(ec,j) = pix + FIT(ec,j) = fit + USER(ec,j) = user + FWIDTH(ec,j) = width + FTYPE(ec,j) = type + EC_NEWFEATURES(ec) = YES + } else if (abs (fit-user) < abs (FIT(ec,j)-USER(ec,j))) { + APN(ec,j) = ap + LINE(ec,j) = ec_line (ec, ap) + ORDER(ec,j) = ORDERS(ec,LINE(ec,j)) + PIX(ec,j) = pix + FIT(ec,j) = fit + USER(ec,j) = user + FWIDTH(ec,j) = width + FTYPE(ec,j) = type + EC_NEWFEATURES(ec) = YES + } + + EC_CURRENT(ec) = j +end diff --git a/noao/onedspec/ecidentify/ecnext.x b/noao/onedspec/ecidentify/ecnext.x new file mode 100644 index 00000000..1028371d --- /dev/null +++ b/noao/onedspec/ecidentify/ecnext.x @@ -0,0 +1,23 @@ +include "ecidentify.h" + +# EC_NEXT -- Return the next feature. + +int procedure ec_next (ec, feature) + +pointer ec # ID pointer +int feature # Starting feature (input), next feature (returned) + +int i + +begin + for (i=feature+1; i<=EC_NFEATURES(ec); i=i+1) + if (APN(ec,i) == EC_AP(ec)) + break + + if (i <= EC_NFEATURES(ec)) + feature = i + else + i = EOF + + return (i) +end diff --git a/noao/onedspec/ecidentify/ecpeak.x b/noao/onedspec/ecidentify/ecpeak.x new file mode 100644 index 00000000..f797fbac --- /dev/null +++ b/noao/onedspec/ecidentify/ecpeak.x @@ -0,0 +1,24 @@ +include "ecidentify.h" + +# EC_PEAK -- Find the peak value above continuum. + +double procedure ec_peak (ec, pix) + +pointer ec # 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 - EC_FWIDTH(ec))) + u = min (EC_NPTS(ec), nint (pix + EC_FWIDTH(ec))) + peak = IMDATA(ec,c) - (IMDATA(ec,l) + + IMDATA(ec,u)) / 2. + + return (peak) +end diff --git a/noao/onedspec/ecidentify/ecprevious.x b/noao/onedspec/ecidentify/ecprevious.x new file mode 100644 index 00000000..4301b722 --- /dev/null +++ b/noao/onedspec/ecidentify/ecprevious.x @@ -0,0 +1,23 @@ +include "ecidentify.h" + +# EC_PREVIOUS -- Return the previous feature. + +int procedure ec_previous (ec, feature) + +pointer ec # ID pointer +int feature # Starting feature (input), previous feature (returned) + +int i + +begin + for (i=feature-1; i>0; i=i-1) + if (APN(ec,i) == EC_AP(ec)) + break + + if (i > 0) + feature = i + else + i = EOF + + return (i) +end diff --git a/noao/onedspec/ecidentify/ecrms.x b/noao/onedspec/ecidentify/ecrms.x new file mode 100644 index 00000000..de84ae26 --- /dev/null +++ b/noao/onedspec/ecidentify/ecrms.x @@ -0,0 +1,28 @@ +include "ecidentify.h" + +# EC_RMS -- Compute RMS of fit about the user coordinates + +double procedure ec_rms (ec) + +pointer ec # ID pointer + +int i, nrms +double rms + +begin + rms = 0. + nrms = 0 + for (i=1; i<=EC_NFEATURES(ec); i=i+1) { + if (!IS_INDEFD (USER(ec,i)) && FTYPE(ec,i) > 0) { + rms = rms + (FIT(ec,i) - USER(ec,i)) ** 2 + nrms = nrms + 1 + } + } + + if (nrms > 0) + rms = sqrt (rms / nrms) + else + rms = INDEFD + + return (rms) +end diff --git a/noao/onedspec/ecidentify/ecshift.x b/noao/onedspec/ecidentify/ecshift.x new file mode 100644 index 00000000..22b050a7 --- /dev/null +++ b/noao/onedspec/ecidentify/ecshift.x @@ -0,0 +1,77 @@ +include <smw.h> +include "ecidentify.h" + +define NBIN 10 # Bin parameter for mode determination + +# EC_SHIFT -- Determine a shift by correlating feature user positions +# with peaks in the image data. + +double procedure ec_shift (ec) + +pointer ec # EC pointer + +int i, j, k, ap, order, nx, ndiff, find_peaks() +real d, dmin +double pix, ec_center(), ec_fitpt() +pointer x, y, diff +errchk malloc, realloc, find_peaks + +begin + ndiff = 0 + call malloc (x, EC_NCOLS(ec), TY_REAL) + call malloc (y, EC_NCOLS(ec), TY_DOUBLE) + do k = 1, EC_NLINES(ec) { + call ec_gline (ec, k) + ap = APS(ec,k) + order = ORDERS(ec,k) + + # Find the peaks in the image data. + i = max (5, EC_MAXFEATURES(ec) / EC_NLINES(ec)) + nx = find_peaks (IMDATA(ec,1), Memr[x], EC_NPTS(ec), 0., + int (EC_MINSEP(ec)), 0, i, 0., false) + + # Center the peaks and convert to user coordinates. + j = 0 + do i = 1, nx { + pix = Memr[x+i-1] + pix = ec_center (ec, pix, EC_FWIDTH(ec), EC_FTYPE(ec)) + if (!IS_INDEFD (pix)) { + Memd[y+j] = ec_fitpt (ec, ap, pix) + j = j + 1 + } + } + nx = j + + # Compute differences with feature list. + do i = 1, EC_NFEATURES(ec) { + if (APN(ec,i) != ap) + next + if (ndiff == 0) + call malloc (diff, nx, TY_REAL) + else + call realloc (diff, ndiff+nx, TY_REAL) + do j = 1, nx { + Memr[diff+ndiff] = (Memd[y+j-1] - FIT(ec,i)) * order + 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/onedspec/ecidentify/ecshow.x b/noao/onedspec/ecidentify/ecshow.x new file mode 100644 index 00000000..e8fb5acc --- /dev/null +++ b/noao/onedspec/ecidentify/ecshow.x @@ -0,0 +1,78 @@ +include <pkg/center1d.h> +include "ecidentify.h" + +# EC_SHOW -- Show parameter information. + +procedure ec_show (ec, file) + +pointer ec # ID pointer +char file[ARB] # File + +char line[SZ_LINE] +int fd + +int open(), ecf_geti() +double ecf_getd() +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[EC_IMAGE(ec)]) + switch (EC_FTYPE(ec)) { + case EMISSION: + call fprintf (fd, "ftype emission\n") + case ABSORPTION: + call fprintf (fd, "ftype absorption\n") + } + switch (EC_LABELS(ec)) { + case 2: + call fprintf (fd, "labels index\n") + case 3: + call fprintf (fd, "labels pixel\n") + case 4: + call fprintf (fd, "labels user\n") + default: + call fprintf (fd, "labels none\n") + } + call fprintf (fd, "maxfeatures %d\n") + call pargi (EC_MAXFEATURES(ec)) + call fprintf (fd, "match %g\n") + call pargr (EC_MATCH(ec)) + call fprintf (fd, "zwidth %g\n") + call pargr (EC_ZWIDTH(ec)) + call fprintf (fd, "fwidth %g\n") + call pargr (EC_FWIDTH(ec)) + call fprintf (fd, "database %s\n") + call pargstr (Memc[EC_DATABASE(ec)]) + call fprintf (fd, "coordlist %s\n") + call pargstr (Memc[EC_COORDLIST(ec)]) + call fprintf (fd, "cradius %g\n") + call pargr (EC_CRADIUS(ec)) + call fprintf (fd, "threshold %g\n") + call pargr (EC_THRESHOLD(ec)) + call fprintf (fd, "minsep %g\n") + call pargr (EC_MINSEP(ec)) + if (EC_ECF(ec) != NULL) { + call fprintf (fd, "function = %s\n") + call ecf_gets ("function", line, SZ_LINE) + call pargstr (line) + call fprintf (fd, "xorder = %d, yorder = %d\n") + call pargi (ecf_geti ("xorder")) + call pargi (ecf_geti ("yorder")) + call fprintf (fd, + "niterate = %d, lowreject = %g, highreject = %g\n") + call pargi (ecf_geti ("niterate")) + call pargd (ecf_getd ("low")) + call pargd (ecf_getd ("high")) + call fprintf (fd, "Fit at first pixel = %0.8g\n") + call pargd (Memd[EC_FITDATA(ec)]) + } + + call close (fd) +end diff --git a/noao/onedspec/ecidentify/mkpkg b/noao/onedspec/ecidentify/mkpkg new file mode 100644 index 00000000..1c8664a7 --- /dev/null +++ b/noao/onedspec/ecidentify/mkpkg @@ -0,0 +1,39 @@ +# ECIDENTIFY Task + +$checkout libpkg.a .. +$update libpkg.a +$checkin libpkg.a .. +$exit + +libpkg.a: + @ecffit + + eccenter.x ecidentify.h + eccolon.x ecidentify.h <error.h> <gset.h> <pkg/center1d.h> + ecdb.x ecidentify.h <math/gsurfit.h> <smw.h> <units.h> + ecdelete.x ecidentify.h + ecdofit.x ecidentify.h <smw.h> + ecdoshift.x ecidentify.h + ecfitdata.x ecidentify.h <pkg/gtools.h> <smw.h> <units.h> + ecgdata.x ecidentify.h <imhdr.h> <imio.h> <pkg/gtools.h> <smw.h>\ + <units.h> + ecgetim.x + ecgline.x ecidentify.h <smw.h> + ecgraph.x ecidentify.h <gset.h> <pkg/gtools.h> + ecidentify.x ecidentify.h <error.h> <gset.h> <imhdr.h> <smw.h> + ecinit.x ecidentify.h <gset.h> + ecline.x ecidentify.h <smw.h> + eclinelist.x ecidentify.h <error.h> <mach.h> <smw.h> <units.h> + eclog.x ecidentify.h <time.h> + ecmark.x ecidentify.h <gset.h> <pkg/center1d.h> + ecnearest.x ecidentify.h <mach.h> + ecnewfeature.x ecidentify.h <mach.h> <smw.h> + ecnext.x ecidentify.h + ecpeak.x ecidentify.h + ecprevious.x ecidentify.h + ecrms.x ecidentify.h + ecshift.x ecidentify.h <smw.h> + ecshow.x ecidentify.h <pkg/center1d.h> + t_eciden.x ecidentify.h <mach.h> <pkg/center1d.h> <pkg/gtools.h> + t_ecreid.x ecidentify.h <error.h> <smw.h> + ; diff --git a/noao/onedspec/ecidentify/t_eciden.x b/noao/onedspec/ecidentify/t_eciden.x new file mode 100644 index 00000000..7590dc17 --- /dev/null +++ b/noao/onedspec/ecidentify/t_eciden.x @@ -0,0 +1,68 @@ +include <mach.h> +include <pkg/gtools.h> +include <pkg/center1d.h> +include "ecidentify.h" + +# T_ECIDENTIFY -- Identify features in echelle format data. +# +# The input data must be in the echelle format produced by APEXTRACT. + +procedure t_ecidentify () + +int images +pointer ec, gopen(), gt_init(), un_open() +int clgeti(), clgwrd(), imtopenp(), ec_getim() +real clgetr() +double clgetd() + +begin + # Allocate the basic data structure. + call ec_init (ec) + + # Get task parameters. + images = imtopenp ("images") + EC_MAXFEATURES(ec) = clgeti ("maxfeatures") + EC_MINSEP(ec) = clgetr ("minsep") + EC_MATCH(ec) = clgetr ("match") + EC_ZWIDTH(ec) = clgetr ("zwidth") + EC_FTYPE(ec) = clgwrd ("ftype", Memc[EC_IMAGE(ec)], SZ_FNAME, FTYPES) + EC_FWIDTH(ec) = clgetr ("fwidth") + EC_CRADIUS(ec) = clgetr ("cradius") + EC_THRESHOLD(ec) = clgetr ("threshold") + call clgstr ("database", Memc[EC_DATABASE(ec)], SZ_FNAME) + call clgstr ("coordlist", Memc[EC_COORDLIST(ec)], SZ_FNAME) + + # Get the line list. + call clgstr ("units", Memc[EC_IMAGE(ec)], SZ_FNAME) + call xt_stripwhite (Memc[EC_IMAGE(ec)]) + if (Memc[EC_IMAGE(ec)] != EOS) + EC_UN(ec) = un_open (Memc[EC_IMAGE(ec)]) + call ec_mapll (ec) + + # Initialize graphics and fitting. + call clgstr ("function", Memc[EC_IMAGE(ec)], SZ_FNAME) + call ecf_sets ("function", Memc[EC_IMAGE(ec)]) + call ecf_seti ("xorder", clgeti ("xorder")) + call ecf_seti ("yorder", clgeti ("yorder")) + call ecf_seti ("niterate", clgeti ("niterate")) + call ecf_setd ("low", clgetd ("lowreject")) + call ecf_setd ("high", clgetd ("highreject")) + call ecf_seti ("xtype", 'p') + call ecf_seti ("ytype", 'r') + call clgstr ("graphics", Memc[EC_IMAGE(ec)], SZ_FNAME) + EC_GP(ec) = gopen (Memc[EC_IMAGE(ec)], NEW_FILE, STDGRAPH) + EC_GT(ec) = gt_init() + call gt_sets (EC_GT(ec), GTTYPE, "line") + + # Identify features in each image. + while (ec_getim (images, Memc[EC_IMAGE(ec)], SZ_FNAME) != EOF) + call ec_identify (ec) + + # Finish up. + call gclose (EC_GP(ec)) + call gt_free (EC_GT(ec)) + call dgsfree (EC_ECF(ec)) + call imtclose (images) + call ec_unmapll (ec) + call ec_free (ec) +end diff --git a/noao/onedspec/ecidentify/t_ecreid.x b/noao/onedspec/ecidentify/t_ecreid.x new file mode 100644 index 00000000..5e9769ff --- /dev/null +++ b/noao/onedspec/ecidentify/t_ecreid.x @@ -0,0 +1,181 @@ +include <error.h> +include <smw.h> +include "ecidentify.h" + +# T_ECREIDENTIFY -- Reidentify echelle features starting from reference. +# If no initial shift is given then the procedure ec_shift computes a +# shift between the reference features and the features in the image. +# The purpose of the shift is to get the feature positions from the +# reference image close enough to those of the image being identified +# that the centering algorithm will determine the exact positions of the +# features. The recentered features are then fit with either a shift +# of a full echelle function and written to database. + +procedure t_ecreidentify () + +int images # List of images +pointer ref # Reference image +double shift # Initial shift + +int i, j, fd, nfeatures1, nfeatures2 +double shift1, pix, fit, pix_shift, fit_shift, z_shift +pointer sp, log, ec + +int imtopenp(), ec_getim(), clpopnu(), clgfil(), open(), btoi() +double ec_fitpt(), ec_fittopix(), ec_shift(), ec_center(), ec_rms() +double clgetd() +bool clgetb() +real clgetr() +errchk ec_dbread(), ec_gdata(), ec_fitdata() + +begin + call smark (sp) + call salloc (ref, SZ_FNAME, TY_CHAR) + call salloc (log, SZ_FNAME, TY_CHAR) + + # Allocate the basic data structure. + call ec_init (ec) + + # Initialize fitting + call ecf_seti ("niterate", 0) + call ecf_setd ("low", 3.D0) + call ecf_setd ("high", 3.D0) + + # Get task parameters. + images = imtopenp ("images") + call clgstr ("reference", Memc[ref], SZ_FNAME) + shift = clgetd ("shift") + call clgstr ("database", Memc[EC_DATABASE(ec)], SZ_FNAME) + EC_CRADIUS(ec) = clgetr ("cradius") + EC_THRESHOLD(ec) = clgetr ("threshold") + EC_LOGFILES(ec) = clpopnu ("logfiles") + EC_REFIT(ec) = btoi (clgetb ("refit")) + + # Write logfile header. + while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) { + iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) { + call erract (EA_WARN) + next + } + call sysid (Memc[log], SZ_LINE) + call fprintf (fd, "\nECREIDENTIFY: %s\n") + call pargstr (Memc[log]) + call fprintf (fd, + " Reference image = %s, Refit = %b\n") + call pargstr (Memc[ref]) + call pargb (EC_REFIT(ec) == YES) + call fprintf (fd, "%20s %7s %7s %9s %10s %7s %7s\n") + call pargstr ("Image") + call pargstr ("Found") + call pargstr ("Fit") + call pargstr ("Pix Shift") + call pargstr ("User Shift") + call pargstr ("Z Shift") + call pargstr ("RMS") + call close (fd) + } + + # Reidentify features in each spectrum. + while (ec_getim (images, Memc[EC_IMAGE(ec)], SZ_FNAME) != EOF) { + call ec_gdata (ec) + call ec_dbread (ec, Memc[ref], NO) + call ec_fitdata (ec) + call ec_fitfeatures (ec) + + if (IS_INDEFD (shift)) { + EC_FWIDTH(ec) = FWIDTH(ec,1) + EC_FTYPE(ec) = abs (FTYPE(ec,1)) + EC_MINSEP(ec) = 1. + EC_MAXFEATURES(ec) = 20 + shift1 = ec_shift (ec) + } else + shift1 = shift + + # Recenter features. + pix_shift = 0. + fit_shift = 0. + z_shift = 0. + nfeatures1 = EC_NFEATURES(ec) + + j = 0. + do i = 1, EC_NFEATURES(ec) { + call ec_gline (ec, LINE(ec,i)) + pix = ec_fittopix (ec, FIT(ec,i) + shift1/ORDER(ec,i)) + pix = ec_center (ec, pix, FWIDTH(ec,i), FTYPE(ec,i)) + if (IS_INDEFD (pix)) + next + fit = ec_fitpt (ec, APN(ec,i), pix) + + pix_shift = pix_shift + pix - PIX(ec,i) + fit_shift = fit_shift + (fit - FIT(ec,i)) * ORDER(ec,i) + if (FIT(ec,i) != 0.) + z_shift = z_shift + (fit - FIT(ec,i)) / FIT(ec,i) + + j = j + 1 + APN(ec,j) = APN(ec,i) + LINE(ec,j) = LINE(ec,i) + ORDER(ec,j) = ORDER(ec,i) + PIX(ec,j) = pix + FIT(ec,j) = FIT(ec,i) + USER(ec,j) = USER(ec,i) + FWIDTH(ec,j) = FWIDTH(ec,i) + FTYPE(ec,j) = abs (FTYPE(ec,i)) + } + EC_NFEATURES(ec) = j + + # If refitting the coordinate function is requested and there + # is more than one feature and there is a previously defined + # coordinate function then refit. Otherwise compute a coordinate + # shift. + + if ((EC_REFIT(ec)==YES)&&(EC_NFEATURES(ec)>1)&&(EC_ECF(ec)!=NULL)) { + iferr (call ec_dofit (ec, NO, YES)) { + call erract (EA_WARN) + next + } + } else + call ec_doshift (ec, NO) + if (EC_NEWECF(ec) == YES) + call ec_fitfeatures (ec) + + nfeatures2 = 0 + do i = 1, EC_NFEATURES(ec) + if (FTYPE(ec,i) > 0) + nfeatures2 = nfeatures2 + 1 + + # Write a database entry for the reidentified image. + if (EC_NFEATURES(ec) > 0) + call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], NO) + + # Record log information if a log file descriptor is given. + call clprew (EC_LOGFILES(ec)) + while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) { + iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) { + call erract (EA_WARN) + next + } + call fprintf (fd, + "%20s %3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n") + call pargstr (Memc[EC_IMAGE(ec)]) + call pargi (EC_NFEATURES(ec)) + call pargi (nfeatures1) + call pargi (nfeatures2) + call pargi (EC_NFEATURES(ec)) + call pargd (pix_shift / max (1, EC_NFEATURES(ec))) + call pargd (fit_shift / max (1, EC_NFEATURES(ec))) + call pargd (z_shift / max (1, EC_NFEATURES(ec))) + call pargd (ec_rms(ec)) + call close (fd) + } + + call smw_close (MW(EC_SH(ec))) + do i = 1, EC_NLINES(ec) + MW(SH(ec,i)) = NULL + } + + call dgsfree (EC_ECF(ec)) + call clpcls (EC_LOGFILES(ec)) + call ec_free (ec) + call imtclose (images) + call sfree (sp) +end |