diff options
Diffstat (limited to 'noao/astcat/src/agetcat')
| -rw-r--r-- | noao/astcat/src/agetcat/atcatinit.x | 167 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atfcat.x | 1986 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/athedit.x | 614 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atincat.x | 70 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atoutcat.x | 72 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atrcquery.x | 522 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atrcrd.x | 314 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atrcsym.x | 29 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/attquery.x | 183 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atwcat.x | 197 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/atwedit.x | 83 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/mkpkg | 31 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_aclist.x | 112 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_afiltcat.x | 211 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_agetcat.x | 251 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_agetim.x | 247 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_ahedit.x | 175 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_aimfind.x | 318 | ||||
| -rw-r--r-- | noao/astcat/src/agetcat/t_aslist.x | 102 |
19 files changed, 5684 insertions, 0 deletions
diff --git a/noao/astcat/src/agetcat/atcatinit.x b/noao/astcat/src/agetcat/atcatinit.x new file mode 100644 index 00000000..04e9dff6 --- /dev/null +++ b/noao/astcat/src/agetcat/atcatinit.x @@ -0,0 +1,167 @@ + +# AT_AGINIT -- Inititialize the AGETCAT task structure. + +procedure at_aginit (at) + +pointer at #O the pointer to the astrometry descriptor + +begin + # Initialize the astrometry structure. + call at_ainit (at) + + # Initialize the i/o structure. + call at_ioinit (at) + + # Initialize the region definition structure. + call at_rcinit (at) + + # Initialize the filtering / selection structure. + call at_fsinit (at) +end + + +# AT_AGFREE -- Free the AGETCAT task structure. + +procedure at_agfree (at) + +pointer at #U the pointer to the astrometry descriptor + +begin + # Free the filtering / selection structure. + call at_fsfree (at) + + # Free the field center structure. + call at_rcfree (at) + + # Free the i/o structure + call at_iofree (at) + + # Free the astrometry structure. + call at_afree (at) +end + + +# AT_AFINIT -- Inititialize the AFILTCAT task structure. + +procedure at_afinit (at) + +pointer at #O the pointer to the astrometry descriptor + +begin + # Initialize the astrometry structure. + call at_ainit (at) + + # Initialize the i/o structure. + call at_ioinit (at) + + # Initialize the filtering / selection structure. + call at_fsinit (at) +end + + +# AT_AFFREE -- Free the AFILTCAT task structure. + +procedure at_affree (at) + +pointer at #U the pointer to the astrometry descriptor + +begin + # Free the filtering / selection structure. + call at_fsfree (at) + + # Free the i/o structure + call at_iofree (at) + + # Free the astrometry structure. + call at_afree (at) +end + + +# AT_AIGINIT -- Inititialize the agetim task structure. + +procedure at_aiginit (at) + +pointer at #O the pointer to the astrometry descriptor + +begin + # Initialize the astrometry structure. + call at_ainit (at) + + # Initialize the i/o structure. + call at_ioinit (at) + + # Initialize the region definition structure. + call at_rcinit (at) + + # Initialize the default wcs structure. + #call at_wcinit (at) + + # Initialize the default image data structure. + #call at_iminit (at) +end + + +# AT_AIGFREE -- Free the agetim task structure. + +procedure at_aigfree (at) + +pointer at #U the pointer to the astrometry descriptor + +begin + # Free the default image data structure. + #call at_imfree (at) + + # Free the default wcs structure. + #call at_wcfree (at) + + # Free the field center structure. + call at_rcfree (at) + + # Free the i/o structure + call at_iofree (at) + + # Free the astrometry structure. + call at_afree (at) +end + + +# AT_AHINIT -- Inititialize the AHEDIT task structure. + +procedure at_ahinit (at) + +pointer at #O the pointer to the astrometry descriptor + +begin + # Initialize the astrometry structure. + call at_ainit (at) + + # Initialize the i/o structure. + call at_ioinit (at) + + # Initialize the default wcs structure. + call at_wcinit (at) + + # Initialize the default image data structure. + call at_iminit (at) +end + + +# AT_AHFREE -- Free the AHEDIT task structure. + +procedure at_ahfree (at) + +pointer at #U the pointer to the astrometry descriptor + +begin + # Free the default image data structure. + call at_imfree (at) + + # Free the default wcs structure. + call at_wcfree (at) + + # Free the i/o structure + call at_iofree (at) + + # Free the astrometry structure. + call at_afree (at) +end diff --git a/noao/astcat/src/agetcat/atfcat.x b/noao/astcat/src/agetcat/atfcat.x new file mode 100644 index 00000000..b2265688 --- /dev/null +++ b/noao/astcat/src/agetcat/atfcat.x @@ -0,0 +1,1986 @@ +include <ctotok.h> +include <ctype.h> +include <evvexpr.h> +include <imhdr.h> +include <pkg/cq.h> +include <pkg/skywcs.h> +include "../../lib/astrom.h" +include "../../lib/acatalog.h" + +############################################################################## + +# Create a small data structure to describe the field list. Decide whether +# this should be part of the main astrometry package structure later ... + +define FL_FLENGTH 12 + +define FL_NEXPR Memi[$1] # The number of input expressions +define FL_NFIELDS Memi[$1+1] # The number of output fields + +# The field list decription +define FL_FLIST Memi[$1+2] # The list of field expressions +define FL_FRANGES Memi[$1+3] # The list of field ranges + +# New quantities to be written in the header (could be a symbol table ...) +define FL_FNAMES Memi[$1+4] # The list of field names +define FL_FOFFSETS Memi[$1+5] # The list of field offsets +define FL_FSIZES Memi[$1+6] # The list of field sizes +define FL_FTYPES Memi[$1+7] # The list of field types +define FL_FUNITS Memi[$1+8] # The list of field units +define FL_FFMTS Memi[$1+9] # The list of field formats + +# Useful constants +define FL_MAX_NEXPR 20 +define FL_MAX_NFIELDS 100 +define FL_SZ_EXPR SZ_LINE + +############################################################################## + +# AT_WIFILRECS -- Filter and write the output catalog. + +procedure at_wifilrecs (fd, im, at, res, standard) + +int fd #I the output file descriptor +pointer im #I the associated input image descriptor +pointer at #I the astrometry package descriptor +pointer res #I results descriptor +bool standard #I write a standard catalog header + +double raval, decval, oraval, odecval, iraval, idecval, xpval, ypval +pointer sp, sexpr, sfield, record, raname, decname, sindex +pointer flist, o, catcoo, outcoo, imcoo, mwim, ct +int i, nlines, nrecs, rafield, decfield, xpfield, ypfield, xcfield, ycfield +pointer at_flinit(), evvexpr(), locpr(), mw_sctran() +int at_wfcathdr(), cq_rstati(), at_srtcat(), at_flnexpr(), cq_grecord() +int cq_setrecord(), at_wcathdr(), at_mkrecord(), cq_gvald() +bool streq() +extern at_getop() + +int nchars + +begin + call smark (sp) + call salloc (sexpr, FL_SZ_EXPR, TY_CHAR) + call salloc (sfield, FL_SZ_EXPR, TY_CHAR) + call salloc (record, SZ_LINE, TY_CHAR) + call salloc (raname, FL_SZ_EXPR, TY_CHAR) + call salloc (decname, FL_SZ_EXPR, TY_CHAR) + + # Initialize the catalog, output, and image coordinate systems. + # and set up the image world to logical coordinate transformation. + call at_cowcs (at, res, im, catcoo, outcoo, imcoo, mwim) + if (imcoo != NULL) + ct = mw_sctran (mwim, "world", "logical", 03B) + + # Determine whether it is necessary to reformat. + call at_stats (at, FIELDS, Memc[record], SZ_LINE) + if (streq (Memc[record], "f[*]") && outcoo == NULL) { + + # The field list is NULL. + flist = NULL + + # Write the filtered catalog header. + if (standard) + nlines = at_wcathdr (fd, res) + + # Coordinate fields are not modified. + rafield = 0 + decfield = 0 + xpfield = 0 + ypfield = 0 + xcfield = 0 + ycfield = 0 + + } else { + + # Decode the output field list. + flist = at_flinit (at, res) + + # Write the filtered catalog header. + if (standard) + nlines = at_wfcathdr (fd, at, res, flist) + + # Get the offsets for the ra, dec, xp, yp, xc, and yc output fields. + call at_coofields (at, res, flist, Memc[raname], Memc[decname], + rafield, decfield, xpfield, ypfield, xcfield, ycfield) + if (outcoo == NULL) { + rafield = 0 + decfield = 0 + } + if (imcoo == NULL) { + xpfield = 0 + ypfield = 0 + } + xcfield = 0 + ycfield = 0 + } + + # Compute the sort index. + nrecs = cq_rstati (res, CQRNRECS) + call malloc (sindex, nrecs, TY_INT) + nrecs = at_srtcat (at, res, Memi[sindex], nrecs) + + # Get the selection expression and replace generic selection expression + # field names with their catalog equivalents. + call at_stats (at, FEXPR, Memc[sexpr], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[sexpr], i, Memc[sfield], FL_SZ_EXPR) == EOF) + Memc[sfield] = EOS + + # Loop over the sorted records. Note that any reference to + # coordinates in the selection expression refers to the original + # not the transformed coordinates. + + o = NULL + do i = 1, nrecs { + + # Reject every record. + if (streq (Memc[sfield], "no")) + next + + # Evaluate the selection expression. + if (! streq (Memc[sfield], "yes")) { + if (cq_setrecord (res, Memi[sindex+i-1]) != Memi[sindex+i-1]) + next + if (o != NULL) + call evvfree (o) + o = evvexpr (Memc[sfield], locpr (at_getop), res, 0, res, 0) + if (O_TYPE(o) != TY_BOOL) + next + if (O_VALI(o) == NO) + next + } + + # Write the record. + if (flist == NULL) { + + # Copy the record. + nchars = cq_grecord (res, Memc[record], SZ_LINE, + Memi[sindex+i-1]) + + } else { + + # Get the ra and dec fields. + raval = INDEFD + decval = INDEFD + if (outcoo != NULL || imcoo != NULL) { + if (cq_gvald (res, Memi[sindex+i-1], Memc[raname], + raval) <= 0) + raval = INDEFD + if (cq_gvald (res, Memi[sindex+i-1], Memc[decname], + decval) <= 0) + decval = INDEFD + } + + # Transform the catalog coordinates to the output coordinate + # system. + oraval = INDEFD + odecval = INDEFD + if (outcoo != NULL && (rafield > 0 || decfield > 0)) { + if (! IS_INDEFD(raval) && ! IS_INDEFD(decval)) + call sk_ultran (catcoo, outcoo, raval, decval, oraval, + odecval, 1) + } + + # Transform the catalog coordinates to the image coordinate + # system and then to the image pixel coordinate system. + xpval = INDEFD + ypval = INDEFD + if (imcoo != NULL && (xpfield > 0 || ypfield > 0)) { + if (! IS_INDEFD(raval) && ! IS_INDEFD(decval)) { + call sk_ultran (catcoo, imcoo, raval, decval, iraval, + idecval, 1) + call mw_c2trand (ct, iraval, idecval, xpval, ypval) + if (xpval < 0.5d0 || xpval > (IM_LEN(im,1)+0.5d0) || + ypval < 0.5d0 || ypval > (IM_LEN(im,2)+0.5d0)) + next + } else { + xpval = INDEFD + ypval = INDEFD + } + } + + # Reformat the record. + nchars = at_mkrecord (flist, res, Memc[record], SZ_LINE, + Memi[sindex+i-1], rafield, decfield, oraval, odecval, + xpfield, ypfield, xpval, ypval) + + } + + # Write the new record. + if (nchars > 0) { + call fprintf (fd, "%s") + call pargstr (Memc[record]) + } + } + + # Free the selection expression descriptor. + if (o != NULL) + call evvfree (o) + + # Free the catalog, output, and image coordinate system descriptors. + if (catcoo != NULL) + call sk_close (catcoo) + if (outcoo != NULL) + call sk_close (outcoo) + if (imcoo != NULL) + call sk_close (imcoo) + if (mwim != NULL) + call mw_close (mwim) + + # Free output field list. + if (flist != NULL) + call at_flfree (flist) + + # Free thesort index descriptor. + call mfree (sindex, TY_INT) + + call sfree (sp) +end + + +# AT_WFILRECS -- Filter and write the output catalog. + +procedure at_wfilrecs (fd, at, res, standard) + +int fd #I the output file descriptor +pointer at #I the astrometry package descriptor +pointer res #I results descriptor +bool standard #I write a standard catalog header + +double raval, decval, oraval, odecval, iraval, idecval, xpval, ypval +pointer sp, sexpr, sfield, record, raname, decname, sindex +pointer flist, o, catcoo, outcoo, imcoo, mwim, ct +int i, nlines, nrecs, rafield, decfield, xpfield, ypfield, xcfield, ycfield +pointer at_flinit(), evvexpr(), locpr(), mw_sctran() +int at_wfcathdr(), cq_rstati(), at_srtcat(), at_flnexpr(), cq_grecord() +int cq_setrecord(), at_wcathdr(), at_mkrecord(), cq_gvald() +bool streq() +extern at_getop() + +int nchars + +begin + call smark (sp) + call salloc (sexpr, FL_SZ_EXPR, TY_CHAR) + call salloc (sfield, FL_SZ_EXPR, TY_CHAR) + call salloc (record, SZ_LINE, TY_CHAR) + call salloc (raname, FL_SZ_EXPR, TY_CHAR) + call salloc (decname, FL_SZ_EXPR, TY_CHAR) + + # Initialize the catalog, output, and field coordinate systems + # and set up the image world to logical coordinate transformation. + call at_cowcs (at, res, NULL, catcoo, outcoo, imcoo, mwim) + if (imcoo != NULL) + ct = mw_sctran (mwim, "world", "logical", 03B) + + # Determine whether it is necessary to reformat. + call at_stats (at, FIELDS, Memc[record], SZ_LINE) + if (streq (Memc[record], "f[*]") && outcoo == NULL) { + + # The field list is NULL. + flist = NULL + + # Write the filtered catalog header. + if (standard) + nlines = at_wcathdr (fd, res) + + # Coordinate fields are not altered. + rafield = 0 + decfield = 0 + xpfield = 0 + ypfield = 0 + xcfield = 0 + ycfield = 0 + + } else { + + # Decode the output field list. + flist = at_flinit (at, res) + + # Write the filtered catalog header. + if (standard) + nlines = at_wfcathdr (fd, at, res, flist) + + # Get the offsets for the ra, dec, xp, yp, xc, and yc output fields. + call at_coofields (at, res, flist, Memc[raname], Memc[decname], + rafield, decfield, xpfield, ypfield, xcfield, ycfield) + if (outcoo == NULL) { + rafield = 0 + decfield = 0 + } + if (imcoo == NULL) { + xpfield = 0 + ypfield = 0 + } + xcfield = 0 + ycfield = 0 + } + + # Compute the sort index. + nrecs = cq_rstati (res, CQRNRECS) + call malloc (sindex, nrecs, TY_INT) + nrecs = at_srtcat (at, res, Memi[sindex], nrecs) + + # Get the selection expression and replace generic selection expression + # field names with their catalog equivalents. + call at_stats (at, FEXPR, Memc[sexpr], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[sexpr], i, Memc[sfield], FL_SZ_EXPR) == EOF) + Memc[sfield] = EOS + + # Loop over the sorted records. Note that any reference to + # coordinates in the selection expression refers to the original + # not the transformed coordinates. + + o = NULL + do i = 1, nrecs { + + # Reject every record. + if (streq (Memc[sfield], "no")) + next + + # Evaluate the selection expression. + if (! streq (Memc[sfield], "yes")) { + if (cq_setrecord (res, Memi[sindex+i-1]) != Memi[sindex+i-1]) + next + if (o != NULL) + call evvfree (o) + o = evvexpr (Memc[sfield], locpr (at_getop), res, 0, res, 0) + if (O_TYPE(o) != TY_BOOL) + next + if (O_VALI(o) == NO) + next + } + + # Write the record. + if (flist == NULL) { + + # Copy the record. + nchars = cq_grecord (res, Memc[record], SZ_LINE, + Memi[sindex+i-1]) + + } else { + + # Get the ra and dec fields. + raval = INDEFD + decval = INDEFD + if (outcoo != NULL || imcoo != NULL) { + if (cq_gvald (res, Memi[sindex+i-1], Memc[raname], + raval) <= 0) + raval = INDEFD + if (cq_gvald (res, Memi[sindex+i-1], Memc[decname], + decval) <= 0) + decval = INDEFD + } + + # Transform the catalog coordinates to the output coordinate + # system. + oraval = INDEFD + odecval = INDEFD + if (outcoo != NULL && (rafield > 0 || decfield > 0)) { + if (! IS_INDEFD(raval) && ! IS_INDEFD(decval)) + call sk_ultran (catcoo, outcoo, raval, decval, oraval, + odecval, 1) + } + + # Transform the catalog coordinates to the image coordinate + # system and then to the image pixel coordinate system. + xpval = INDEFD + ypval = INDEFD + if (imcoo != NULL && (xpfield > 0 || ypfield > 0)) { + if (! IS_INDEFD(raval) && ! IS_INDEFD(decval)) { + call sk_ultran (catcoo, imcoo, raval, decval, iraval, + idecval, 1) + call mw_c2trand (ct, iraval, idecval, xpval, ypval) + } else { + xpval = INDEFD + ypval = INDEFD + } + } + + # Reformat the record. + nchars = at_mkrecord (flist, res, Memc[record], SZ_LINE, + Memi[sindex+i-1], rafield, decfield, oraval, odecval, + xpfield, ypfield, xpval, ypval) + + } + + # Write the new record. + if (nchars > 0) { + call fprintf (fd, "%s") + call pargstr (Memc[record]) + } + } + + # Free the selection expression descriptor. + if (o != NULL) + call evvfree (o) + + # Free the catalog, output, and field coordinate system descriptors. + if (catcoo != NULL) + call sk_close (catcoo) + if (outcoo != NULL) + call sk_close (outcoo) + if (imcoo != NULL) + call sk_close (imcoo) + if (mwim != NULL) + call mw_close (mwim) + + # Free output field list. + if (flist != NULL) + call at_flfree (flist) + + # Free thesort index descriptor. + call mfree (sindex, TY_INT) + + call sfree (sp) +end + + +# AT_FCATHDR -- Write the filtered catalog header + +int procedure at_wfcathdr (fd, at, res, fl) + +int fd #I the output file descriptor +pointer at #I the astrometry pacakge descriptor +pointer res #I the results descriptor descriptor +pointer fl #I the output field list descriptor + +pointer sp, catname, qpnames, qpvalues, qpunits, fname, fvalue, funits +int i, nlines, nfields +int cq_rstati(), at_wrdstr(), cq_hinfon() +char cq_itype() +bool streq(), strne() + +begin + nlines = 0 + + # Allocate working space. + call smark (sp) + call salloc (catname, SZ_FNAME, TY_CHAR) + call salloc (qpnames, SZ_LINE, TY_CHAR) + call salloc (qpvalues, SZ_LINE, TY_CHAR) + call salloc (qpunits, SZ_LINE, TY_CHAR) + call salloc (fname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (fvalue, CQ_SZ_QPVALUE, TY_CHAR) + call salloc (funits, CQ_SZ_QPUNITS, TY_CHAR) + + # Write the header banner. + call fprintf (fd, "# BEGIN CATALOG HEADER\n") + nlines = nlines + 1 + + # Write the catalog database and id. + call cq_rstats (res, CQRCATDB, Memc[catname], SZ_FNAME) + call fprintf (fd, "# catdb %s\n") + call pargstr (Memc[catname]) + nlines = nlines + 1 + call cq_rstats (res, CQRCATNAME, Memc[catname], SZ_FNAME) + call fprintf (fd, "# catname %s\n") + call pargstr (Memc[catname]) + nlines = nlines + 1 + + # Write out the query parameter names, values, and units used + # to generate the catalog. + call cq_rstats (res, CQRQPNAMES, Memc[qpnames], SZ_LINE) + call cq_rstats (res, CQRQPVALUES, Memc[qpvalues], SZ_LINE) + call cq_rstats (res, CQRQPUNITS, Memc[qpunits], SZ_LINE) + nfields = cq_rstati (res, CQRNQPARS) + call fprintf (fd, "# nquery %d\n") + call pargi (nfields) + nlines = nlines + 1 + do i = 1, nfields { + if (at_wrdstr (i, Memc[fname], CQ_SZ_QPNAME, Memc[qpnames]) != i) + ; + if (at_wrdstr (i, Memc[fvalue], CQ_SZ_QPVALUE, Memc[qpvalues]) != i) + ; + if (at_wrdstr (i, Memc[funits], CQ_SZ_QPUNITS, Memc[qpunits]) != i) + ; + call fprintf (fd, "# %s %s %s\n") + call pargstr (Memc[fname]) + call pargstr (Memc[fvalue]) + call pargstr (Memc[funits]) + nlines = nlines + 1 + } + + # Write out the results format type. + if (at_wrdstr (cq_rstati(res, CQRTYPE), Memc[fvalue], CQ_SZ_QPVALUE, + CQ_RTYPESTR) <= 0) + call strcpy ("stext", Memc[fvalue], CQ_SZ_QPVALUE) + call fprintf (fd, "# type %s\n") + call pargstr (Memc[fvalue]) + nlines = nlines + 1 + + # Write out the header parameters, + nfields = cq_rstati (res, CQNHEADER) + call fprintf (fd, "# nheader %d\n") + call pargi (nfields) + nlines = nlines + 1 + do i = 1, nfields { + if (cq_hinfon (res, i, Memc[fname], CQ_SZ_QPNAME, Memc[fvalue], + CQ_SZ_QPVALUE) != i) + next + + # Check for a changed coordinate system here + if (streq ("csystem", Memc[fname])) { + call at_stats (at, FOSYSTEM, Memc[qpvalues], SZ_LINE) + if (Memc[qpvalues] != EOS && strne (Memc[qpvalues], + Memc[fvalue])) + call strcpy (Memc[qpvalues], Memc[fvalue], CQ_SZ_QPVALUE) + } + + call fprintf (fd, "# %s %s\n") + call pargstr (Memc[fname]) + call pargstr (Memc[fvalue]) + nlines = nlines + 1 + } + + # Write out the field desription. + nfields = FL_NFIELDS(fl) + call fprintf (fd, "# nfields %d\n") + call pargi (nfields) + do i = 0, nfields - 1 { + call fprintf (fd, "# %s %d %d %c %s %s\n") + call pargstr (Memc[FL_FNAMES(fl)+i*(CQ_SZ_QPNAME+1)]) + call pargi (Memi[FL_FOFFSETS(fl)+i]) + call pargi (Memi[FL_FSIZES(fl)+i]) + call pargc (cq_itype (Memi[FL_FTYPES(fl)+i])) + call pargstr (Memc[FL_FUNITS(fl)+i*(CQ_SZ_QPUNITS+1)]) + call pargstr (Memc[FL_FFMTS(fl)+i*(CQ_SZ_QPFMTS+1)]) + nlines = nlines + 1 + } + + # Write the header trailer. + call fprintf (fd, "# END CATALOG HEADER\n#\n") + nlines = nlines + 1 + + call sfree (sp) + + return (nlines) +end + + +# AT_SRTCAT -- Sort the catalog on the user specified field. + +int procedure at_srtcat (at, res, sindex, max_nrecs) + +pointer at #I the astrometry package descriptor +pointer res #I results descriptor +int sindex[ARB] #O the output sort index +int max_nrecs #I the maximum number of records + +double dval +pointer sp, sexpr, sfield, sname, sval, darray, carray, o +int i, ip, nrecs, stype, snum, nchars, sz_carray +pointer evvexpr(), locpr() +int ctotok(), cq_fnumber(), cq_ftype(), cq_fname(), ctoi(), cq_gvald() +int cq_gvalc(), gstrcpy(), cq_setrecord(), at_flnexpr(), at_stati() +bool streq() +extern at_getop() + +begin + call smark (sp) + call salloc (sexpr, FL_SZ_EXPR, TY_CHAR) + call salloc (sfield, FL_SZ_EXPR, TY_CHAR) + call salloc (sname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (sval, SZ_LINE, TY_CHAR) + + # Get the sort expression. + call at_stats (at, FSORT, Memc[sexpr], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[sexpr], i, Memc[sfield], FL_SZ_EXPR) == EOF) + Memc[sfield] = EOS + + # Return initialized index array if the sort expression is undefined. + if (Memc[sfield] == EOS) { + do i = 1, max_nrecs + sindex[i] = i + call sfree (sp) + return (max_nrecs) + } + + # Determine the type of sort. If sfield and sname are identical + # sort expression is a field, otherwise it is an expression which + # must be evaluated. + ip = 1 + if (ctotok (Memc[sfield], ip, Memc[sname], CQ_SZ_QPNAME) == + TOK_IDENTIFIER) + ; + + # Initialize the sort index array. + do i = 1, max_nrecs + sindex[i] = i + + # The sort expression is a simple field. + if (streq (Memc[sfield], Memc[sname])) { + + if (cq_fnumber (res, Memc[sfield]) > 0) { # Catalog field name + stype = cq_ftype (res, Memc[sfield]) + } else if (Memc[sfield] == 'f') { # Generic f# name + ip = 2 + if (ctoi (Memc[sfield], ip, snum) <= 0) + stype = INDEFI + else if (cq_fname (res, snum, Memc[sname], CQ_SZ_FNAME) <= 0) + stype = INDEFI + else + stype = cq_ftype (res, Memc[sname]) + } else { # Unknown name. + stype = INDEFI + } + + # Do the sort. + if (IS_INDEFI(stype)) { # Field is undecodable. + nrecs = max_nrecs + + } else if (stype == TY_CHAR) { # Character sort. + sz_carray = 10 * SZ_LINE + call malloc (carray, sz_carray, TY_CHAR) + ip = 1 + do i = 1, max_nrecs { + nchars = cq_gvalc (res, i, Memc[sname], Memc[sval], + SZ_LINE) + if (nchars > sz_carray - ip + 1) { + sz_carray = sz_carray + 10 * SZ_LINE + call realloc (carray, sz_carray, TY_CHAR) + } + sindex[i] = ip + ip = ip + gstrcpy (Memc[sval], Memc[carray+ip-1], nchars) + Memc[carray+ip-1] = EOS + ip = ip + 1 + } + call at_ssquick (Memc[carray], sindex, sindex, max_nrecs) + call mfree (carray, TY_CHAR) + nrecs = max_nrecs + + } else { # Numeric sort. + call malloc (darray, max_nrecs, TY_DOUBLE) + do i = 1, max_nrecs { + nchars = cq_gvald (res, i, Memc[sname], dval) + if (nchars <= 0) + Memd[darray+i-1] = INDEFD + else + Memd[darray+i-1] = dval + } + call at_qsortd (Memd[darray], sindex, sindex, max_nrecs) + call mfree (darray, TY_DOUBLE) + nrecs = max_nrecs + } + + # The sort field is an expression which must be evaluated. + } else { + + # Determine the data type of the output from the first record. + if (cq_setrecord (res, 1) != 1) + ; + o = evvexpr (Memc[sfield], locpr (at_getop), res, 0, res, 0) + stype = O_TYPE(o) + call evvfree (o) + + if (stype == 0) # Expression cannot be decoded. + nrecs = max_nrecs + else if (stype == TY_CHAR || stype == TY_BOOL) { + sz_carray = 10 * SZ_LINE + call malloc (carray, sz_carray, TY_CHAR) + ip = 1 + do i = 1, max_nrecs { + if (cq_setrecord (res, i) != i) + break + o = evvexpr (Memc[sfield], locpr (at_getop), res, 0, res, 0) + if (O_LEN(o) > sz_carray - ip + 1) { + sz_carray = sz_carray + 10 * SZ_LINE + call realloc (carray, sz_carray, TY_CHAR) + } + sindex[i] = ip + ip = ip + gstrcpy (O_VALC(o), Memc[carray+ip-1], O_LEN(o)) + Memc[carray+ip-1] = EOS + ip = ip + 1 + call evvfree (o) + } + + call at_ssquick (Memc[carray], sindex, sindex, max_nrecs) + call mfree (carray, TY_CHAR) + nrecs = max_nrecs + } else { + call malloc (darray, max_nrecs, TY_DOUBLE) + do i = 1, max_nrecs { + if (cq_setrecord (res, i) != i) + break + o = evvexpr (Memc[sfield], locpr (at_getop), res, 0, res, 0) + switch (O_TYPE(o)) { + case TY_SHORT: + dval = O_VALS(o) + case TY_INT: + dval = O_VALI(o) + case TY_LONG: + dval = O_VALL(o) + case TY_REAL: + dval = O_VALR(o) + case TY_DOUBLE: + dval = O_VALD(o) + default: + dval = INDEFD + } + Memd[darray+i-1] = dval + call evvfree (o) + } + call at_qsortd (Memd[darray], sindex, sindex, max_nrecs) + call mfree (darray, TY_DOUBLE) + nrecs = max_nrecs + } + + } + + # Flip the index array if the sense of the sort is reversed. + if (at_stati (at, FREVERSE) == YES) { + do i = 1, nrecs / 2 { + ip = sindex[i] + sindex[i] = sindex[nrecs-i+1] + sindex[nrecs-i+1] = ip + } + } + + call sfree (sp) + + return (nrecs) +end + + +# AT_COWCS -- Initialize the catalog and output coordinate system +# descriptors. + +procedure at_cowcs (at, res, im, catcoo, outcoo, imcoo, mwim) + +pointer at #I the astrometry package descriptor +pointer res #I the catalog results descriptor +pointer im #I the associated image descriptor +pointer catcoo #O the output catalog system descriptor +pointer outcoo #O the output output system descriptor +pointer imcoo #O the output image system descriptor +pointer mwim #O the output image mwcs descriptor + +pointer sp, csystem, cfield, fra, fdec, mw +int i, catstat, outstat, imstat +int cq_hinfo(), sk_decwcs(), strdic(), at_wrdstr(), sk_stati() +int at_stati(), sk_decim() + +begin + call smark (sp) + call salloc (csystem, SZ_LINE, TY_CHAR) + call salloc (cfield, SZ_LINE, TY_CHAR) + call salloc (fra, SZ_FNAME, TY_CHAR) + call salloc (fdec, SZ_FNAME, TY_CHAR) + + # Get the catalog system. + if (cq_hinfo (res, "csystem", Memc[csystem], SZ_LINE) <= 0) + call strcpy ("", Memc[csystem], SZ_LINE) + + # Open the catalog system. + catstat = sk_decwcs (Memc[csystem], mw, catcoo, NULL) + if (catstat == ERR || mw != NULL) { + #call eprintf ( + #"Error: Cannot decode the catalog coordinate system\n") + if (mw != NULL) + call mw_close (mw) + call sk_close (catcoo) + catcoo = NULL + imcoo = NULL + outcoo = NULL + call sfree (sp) + return + } + + + # Get and set the ra catalog coordinate units. + call at_stats (at, FIRA, Memc[fra], SZ_FNAME) + call cq_funits (res, Memc[fra], Memc[cfield], SZ_LINE) + i = strdic (Memc[cfield], Memc[cfield], SZ_LINE, SKY_LNG_UNITLIST) + if (i > 0) + call sk_seti (catcoo, S_NLNGUNITS, i) + + # Get and set the dec catalog coordinate units. + call at_stats (at, FIDEC, Memc[fdec], SZ_FNAME) + call cq_funits (res, Memc[fdec], Memc[cfield], SZ_LINE) + i = strdic (Memc[cfield], Memc[cfield], SZ_LINE, SKY_LAT_UNITLIST) + if (i > 0) + call sk_seti (catcoo, S_NLATUNITS, i) + + # Open the output coordinate system if the output coordinate system is + # different from the catalog coordinate system or the units are + # different. + call at_stats (at, FOSYSTEM, Memc[csystem], SZ_LINE) + if (Memc[csystem] != EOS || at_stati(at,FORAUNITS) > 0 || at_stati(at, + FODECUNITS) > 0) { + + if (Memc[csystem] == EOS) + outstat = sk_decwcs (Memc[csystem], mw, outcoo, catcoo) + else + outstat = sk_decwcs (Memc[csystem], mw, outcoo, NULL) + + if (outstat == ERR || mw != NULL) { + #call eprintf ( + #"Error: Cannot decode the output coordinate system\n") + if (mw != NULL) + call mw_close (mw) + call sk_close (outcoo) + outcoo = NULL + } else { + + # Set the output catalog ra units. + i = at_stati (at, FORAUNITS) + if (i <= 0) { + Memc[cfield] = EOS + } else if (at_wrdstr (i, Memc[cfield], SZ_LINE, + AT_RA_UNITS) <= 0) { + Memc[cfield] = EOS + } + if (Memc[cfield] == EOS) { + call sk_seti (outcoo, S_NLNGUNITS, sk_stati (catcoo, + S_NLNGUNITS)) + } else { + i = strdic (Memc[cfield], Memc[cfield], FL_SZ_EXPR, + SKY_LNG_UNITLIST) + if (i > 0) + call sk_seti (outcoo, S_NLNGUNITS, i) + else + call sk_seti (outcoo, S_NLNGUNITS, sk_stati(catcoo, + S_NLNGUNITS)) + } + + # Set the output catalog dec units. + i = at_stati (at, FODECUNITS) + if (i <= 0) { + Memc[cfield] = EOS + } else if (at_wrdstr (i, Memc[cfield], SZ_LINE, + AT_DEC_UNITS) <= 0) { + Memc[cfield] = EOS + } + if (Memc[cfield] == EOS) { + call sk_seti (outcoo, S_NLATUNITS, sk_stati (catcoo, + S_NLATUNITS)) + } else { + i = strdic (Memc[cfield], Memc[cfield], SZ_LINE, + SKY_LAT_UNITLIST) + if (i > 0) + call sk_seti (outcoo, S_NLATUNITS, i) + else + call sk_seti (outcoo, S_NLATUNITS, sk_stati(catcoo, + S_NLATUNITS)) + } + + } + } else { + outcoo = NULL + } + + # Open the image coordinate system. + if (im == NULL) { + imcoo = NULL + mwim = NULL + } else { + imstat = sk_decim (im, "logical", mwim, imcoo) + if (imstat == ERR || mwim == NULL) { + if (mwim != NULL) + call mw_close (mwim) + mwim = NULL + call sk_close (outcoo) + outcoo = NULL + } else { + call sk_seti (imcoo, S_NLNGUNITS, SKY_DEGREES) + call sk_seti (imcoo, S_NLATUNITS, SKY_DEGREES) + } + } + + call sfree (sp) +end + + +# AT_COOFIELDS -- Get the sequence number of the coordinate output fields. + +procedure at_coofields (at, res, flist, raname, decname, rafield, decfield, + xpfield, ypfield, xcfield, ycfield) + + +pointer at #I the astrometry package descriptor +pointer res #I the output results descriptor +pointer flist #I the output field list descriptor +char raname[ARB] #O the catalog ra name +char decname[ARB] #O the catalog dec name +int rafield #O the output ra field no +int decfield #O the output dec field no +int xpfield #O the output xp field no +int ypfield #O the output yp field no +int xcfield #O the output xp field no +int ycfield #O the output yp field no + +pointer sp, xpname, ypname, str +int i +int at_flnexpr(), cq_fnumber() +bool streq() + +begin + # Get working space. + call smark (sp) + call salloc (xpname, FL_SZ_EXPR, TY_CHAR) + call salloc (ypname, FL_SZ_EXPR, TY_CHAR) + call salloc (str, FL_SZ_EXPR, TY_CHAR) + + # Initialize. + rafield = 0 + decfield = 0 + xpfield = 0 + ypfield = 0 + xcfield = 0 + ycfield = 0 + + # Get the ra and dec field names. + call at_stats (at, FIRA, Memc[str], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[str], i, raname, FL_SZ_EXPR) == EOF) + raname[1] = EOS + call at_stats (at, FIDEC, Memc[str], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[str], i, decname, FL_SZ_EXPR) == EOF) + decname[1] = EOS + + # Get the predicted x and y field names. + call at_stats (at, FIXP, Memc[str], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[str], i, Memc[xpname], FL_SZ_EXPR) == EOF) + Memc[xpname] = EOS + call at_stats (at, FIYP, Memc[str], FL_SZ_EXPR) + i = 1 + if (at_flnexpr (res, Memc[str], i, Memc[ypname], FL_SZ_EXPR) == EOF) + Memc[ypname] = EOS + + # Get the center x and y field names. Ignore this for now. + + # Check to see whether the field names are in the input catalog + # and whether at least one of them is in the output catalog. + if (cq_fnumber (res, raname) > 0 && cq_fnumber (res, decname) > 0) { + do i = 0, FL_NFIELDS(flist) - 1 { + if (streq (raname, Memc[FL_FNAMES(flist)+i* + (CQ_SZ_QPNAME+1)])) { + rafield = i + 1 + } else if (streq (decname, Memc[FL_FNAMES(flist)+i* + (CQ_SZ_QPNAME+1)])) { + decfield = i + 1 + } else if (streq (Memc[xpname], Memc[FL_FNAMES(flist)+i* + (CQ_SZ_QPNAME+1)])) { + xpfield = i + 1 + } else if (streq (Memc[ypname], Memc[FL_FNAMES(flist)+i* + (CQ_SZ_QPNAME+1)])) { + ypfield = i + 1 + } + #if (rafield > 0 && decfield > 0) + #break + } + } + + call sfree (sp) +end + + +# AT_MKRECORD -- Format the output catalog record. + +int procedure at_mkrecord (flist, res, record, maxch, recno, rafield, decfield, + raval, decval, xpfield, ypfield, xpval, ypval) + +pointer flist #I the output field list descriptor +pointer res #I the output results descriptor +char record[ARB] #O the output record +int maxch #I the maximum size of a record +int recno #I the current record number +int rafield #I the output ra field +int decfield #I the output dec field +double raval #I the input ra value +double decval #I the input dec value +int xpfield #I the output predicted x field +int ypfield #I the output predicted y field +double xpval #I the input predicted x value +double ypval #I the input predicted y value + +pointer sp, newval, eptr, rptr, o +int i, j, k, op, findex, nchars +pointer evvexpr(), locpr() +int gstrcpy(), cq_rstati(), cq_gvalc(), strlen() +extern at_getop() + +begin + call smark (sp) + call salloc (newval, SZ_LINE, TY_CHAR) + + # Initialize. + findex = 0 + op = 1 + record[op] = EOS + o = NULL + eptr = FL_FLIST(flist) + rptr = FL_FRANGES(flist) + + # Loop over the expressions. + do i = 1, FL_NEXPR(flist) { + + # The output field is an expression. + if (IS_INDEFI(Memi[rptr])) { + + # Evaluate the expression. + if (xpfield == (findex + 1)) { + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (xpval) + } else if (ypfield == (findex + 1)) { + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (ypval) + } else { + + if (o != NULL) + call evvfree (o) + o = evvexpr (Memc[eptr], locpr (at_getop), res, 0, res, 0) + + # Encode the expression in a string. + switch (O_TYPE(o)) { + case TY_CHAR: + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargstr (O_VALC(o)) + case TY_INT: + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargi (O_VALI(o)) + case TY_REAL: + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargr (O_VALR(o)) + case TY_DOUBLE: + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (O_VALD(o)) + default: + call sprintf (Memc[newval], maxch, Memc[FL_FFMTS(flist)+ + findex * (CQ_SZ_QPFMTS + 1)]) + call pargstr (O_VALC(o)) + } + } + + # Copy the string to the output record. + if (Memi[FL_FSIZES(flist)+findex] == 0) { + op = op + gstrcpy (" ", record[op], maxch - op + 1) + op = op + gstrcpy (Memc[newval], record[op], maxch - op + 1) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + } else { + nchars = min (Memi[FL_FSIZES(flist)+findex], + strlen (Memc[newval])) + do k = 1, Memi[FL_FSIZES(flist)+findex] - nchars - 1 { + if (op > maxch) + break + record[op] = ' ' + op = op + 1 + } + op = op + gstrcpy (Memc[newval], record[op], min (nchars, + maxch - op + 1)) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + } + + findex = findex + 1 + + # The field expression are input catalog columns. + } else if (Memi[rptr] >= 1 && Memi[rptr+1] <= cq_rstati (res, + CQNFIELDS)) { + + # Loop over the fields in each range. + do j = max (1, Memi[rptr]), min (Memi[rptr+1], cq_rstati(res, + CQNFIELDS)), Memi[rptr+2] { + + # Encode the record values. + if (rafield == (findex + 1)) { + call sprintf (Memc[newval], SZ_LINE, + Memc[FL_FFMTS(flist)+ findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (raval) + nchars = strlen (Memc[newval]) + } else if (decfield == (findex + 1)) { + call sprintf (Memc[newval], SZ_LINE, + Memc[FL_FFMTS(flist)+ findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (decval) + nchars = strlen (Memc[newval]) + } else if (xpfield == (findex + 1)) { + call sprintf (Memc[newval], SZ_LINE, + Memc[FL_FFMTS(flist)+ findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (xpval) + nchars = strlen (Memc[newval]) + } else if (ypfield == (findex + 1)) { + call sprintf (Memc[newval], SZ_LINE, + Memc[FL_FFMTS(flist)+ findex * (CQ_SZ_QPFMTS + 1)]) + call pargd (ypval) + nchars = strlen (Memc[newval]) + } else { + nchars = cq_gvalc (res, recno, Memc[FL_FNAMES(flist)+ + findex*(CQ_SZ_QPNAME+1)], Memc[newval], SZ_LINE) + } + + # Copy the string to the output record. + if (Memi[FL_FSIZES(flist)+findex] == 0) { + if ((j == 1) && (! IS_WHITE(Memc[newval]))) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + else if (rafield == (findex + 1) || decfield == + (findex + 1) || xpfield == (findex + 1) || + ypfield == (findex + 1)) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + op = op + gstrcpy (Memc[newval], record[op], + maxch - op + 1) + if (rafield == (findex + 1) || decfield == + (findex + 1) || xpfield == (findex + 1) || + ypfield == (findex + 1)) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + } else { + nchars = min (Memi[FL_FSIZES(flist)+findex], nchars) + do k = 1, Memi[FL_FSIZES(flist)+findex] - nchars - 1 { + if (op > maxch) + break + record[op] = ' ' + op = op + 1 + } + op = op + gstrcpy (Memc[newval], record[op], + min (nchars, maxch - op + 1)) + op = op + gstrcpy (" ", record[op], maxch - op + 1) + } + + findex = findex + 1 + } + } + + # Increment the expression and ranges pointers. + eptr = eptr + FL_SZ_EXPR + 1 + rptr = rptr + 3 + } + if (o != NULL) + call evvfree (o) + + # Append a newline and EOS to the data. + if (record[1] != EOS) { + record[op] = '\n' + record[op+1] = EOS + } + + call sfree (sp) + + return (op - 1) +end + + +# AT_NFLIST -- Add new fields to the current output field list and optionally +# specify the field names, field types, field units, and field formats. + +procedure at_nflist (at, nfields, nfnames, nftypes, nfunits, nformats, append) + +pointer at #I the astrometry package descriptors +int nfields #I the number of new fields +char nfnames[ARB] #I the new field names list +char nftypes[ARB] #I the new field types list +char nfunits[ARB] #I the new field units list +char nformats[ARB] #I the new field formats list +bool append #I append the new fields + +pointer sp, str1, str2 +int i, op1, op2 +int gstrcpy(), strlen + +begin + call smark (sp) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + + # Set the new output field expressions to INDEF and retrieve the + # original user fields value. + call at_stats (at, FIELDS, Memc[str1], SZ_LINE) + op1 = strlen (Memc[str1]) + op2 = 0 + do i = 1, nfields { + if (i == 1) { + op2 = op2 + gstrcpy ("INDEF", Memc[str2+op2], SZ_LINE - op2) + } else { + op2 = op2 + gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy ("INDEF", Memc[str2+op2], SZ_LINE - op2) + } + } + + # Construct the new output fields string. + if (append) { + op1 = op1 + gstrcpy (",", Memc[str1+op1], SZ_LINE - op1) + op1 = op1 + gstrcpy (Memc[str2], Memc[str1+op1], SZ_LINE - op1) + call at_sets (at, FIELDS, Memc[str1]) + } else { + op2 = op 2+ gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (Memc[str1], Memc[str2+op2], SZ_LINE - op2) + call at_sets (at, FIELDS, Memc[str2]) + } + + # Construct the new field names. + call at_stats (at, FNAMES, Memc[str1], SZ_LINE) + op1 = strlen (Memc[str1]) + op2 = 0 + if (append) { + op1 = op1 + gstrcpy (",", Memc[str1+op1], SZ_LINE - op1) + op1 = op1 + gstrcpy (nfnames, Memc[str1+op1], SZ_LINE - op1) + call at_sets (at, FNAMES, Memc[str1]) + } else { + op2 = op2 + gstrcpy (nfnames, Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (Memc[str1], Memc[str2+op2], SZ_LINE - op2) + call at_sets (at, FNAMES, Memc[str2]) + } + + # Construct the new field types. + call at_stats (at, FNTYPES, Memc[str1], SZ_LINE) + op1 = strlen (Memc[str1]) + op2 = 0 + if (append) { + op1 = op1 + gstrcpy (",", Memc[str1+op1], SZ_LINE - op1) + op1 = op1 + gstrcpy (nftypes, Memc[str1+op1], SZ_LINE - op1) + call at_sets (at, FNTYPES, Memc[str1]) + } else { + op2 = op2 + gstrcpy (nftypes, Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (Memc[str1], Memc[str2+op2], SZ_LINE - op2) + call at_sets (at, FNTYPES, Memc[str2]) + } + + # Construct the new field units. + call at_stats (at, FNUNITS, Memc[str1], SZ_LINE) + op1 = strlen (Memc[str1]) + op2 = 0 + if (append) { + op1 = op1 + gstrcpy (",", Memc[str1+op1], SZ_LINE - op1) + op1 = op1 + gstrcpy (nfunits, Memc[str1+op1], SZ_LINE - op1) + call at_sets (at, FNUNITS, Memc[str1]) + } else { + op2 = op2 + gstrcpy (nfunits, Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (Memc[str1], Memc[str2+op2], SZ_LINE - op2) + call at_sets (at, FNUNITS, Memc[str2]) + } + + # Construct the new field units. + call at_stats (at, FNFORMATS, Memc[str1], SZ_LINE) + op1 = strlen (Memc[str1]) + op2 = 0 + if (append) { + op1 = op1 + gstrcpy (",", Memc[str1+op1], SZ_LINE - op1) + op1 = op1 + gstrcpy (nformats, Memc[str1+op1], SZ_LINE - op1) + call at_sets (at, FNFORMATS, Memc[str1]) + } else { + op2 = op2 + gstrcpy (nformats, Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (",", Memc[str2+op2], SZ_LINE - op2) + op2 = op2 + gstrcpy (Memc[str1], Memc[str2+op2], SZ_LINE - op2) + call at_sets (at, FNFORMATS, Memc[str2]) + } + + call sfree (sp) +end + + +# AT_FLINIT -- Initialize the field list structure. This routines: 1) +# creates a list of fields, field ranges, and field expressions, 2) determines +# whether an output field is an input field or a field expression, and 3) +# assembles the information required to write a catalog header. + + +pointer procedure at_flinit (at, res) + +pointer at #I the astrometry package descriptor +pointer res #I results descriptor + +pointer sp, fields +pointer fl +int nexpr, nfields +int at_flelist(), at_flranges() + +begin + # Get some working space. + call smark (sp) + call salloc (fields, SZ_LINE, TY_CHAR) + + # Get the user field list. + call at_stats (at, FIELDS, Memc[fields], SZ_LINE) + + # Allocate the field list descriptor. + call calloc (fl, FL_FLENGTH, TY_STRUCT) + + # Create the field expression list. + nexpr = at_flelist (res, Memc[fields], fl) + if (nexpr > 0) { + + # Determine which individual fields are to be output. + nfields = at_flranges (res, fl) + + # Compile the new header info. + if (nfields > 0) + call at_flflist (at, res, fl) + } + + call sfree (sp) + + return (fl) +end + + +# AT_FLFREE -- Free the field list structure. + +procedure at_flfree (fl) + +pointer fl #I the field list descriptor + +begin + if (FL_FLIST(fl) != NULL) + call mfree (FL_FLIST(fl), TY_CHAR) + if (FL_FRANGES(fl) != NULL) + call mfree (FL_FRANGES(fl), TY_STRUCT) + + if (FL_FNAMES(fl) != NULL) + call mfree (FL_FNAMES(fl), TY_CHAR) + if (FL_FOFFSETS(fl) != NULL) + call mfree (FL_FOFFSETS(fl), TY_INT) + if (FL_FSIZES(fl) != NULL) + call mfree (FL_FSIZES(fl), TY_INT) + if (FL_FTYPES(fl) != NULL) + call mfree (FL_FTYPES(fl), TY_INT) + if (FL_FUNITS(fl) != NULL) + call mfree (FL_FUNITS(fl), TY_CHAR) + if (FL_FFMTS(fl) != NULL) + call mfree (FL_FFMTS(fl), TY_CHAR) + + call mfree (fl, TY_STRUCT) +end + + +# AT_FLELIST -- Create the expression list from the user field list. + +int procedure at_flelist (res, fields, fl) + +pointer res #I the results descriptor +char fields[ARB] #I the user field list +pointer fl #O the field list descriptor + +int i, ip, fp, nexpr +int at_flnexpr() + +begin + # Allocate space for the expression list. + call malloc (FL_FLIST(fl), FL_MAX_NEXPR * (FL_SZ_EXPR + 1), TY_CHAR) + + # Decode the user field list into a list of comma separated + # expressions. Expressions may be field names (e.g. "ra" or "f2"), + # field ranges (e.g. "f[*]" or "f[1-4]"), or field expressions + # (e.g. "f2 - f3" or "mag2 - mag1"). + + ip = 1 + fp = FL_FLIST(fl) + nexpr = 0 + do i = 1, FL_MAX_NEXPR { + if (at_flnexpr (res, fields, ip, Memc[fp], FL_SZ_EXPR) == EOF) + break + #call strlwr (Memc[fp]) + fp = fp + FL_SZ_EXPR + 1 + nexpr = nexpr + 1 + } + call realloc (FL_FLIST(fl), nexpr * (FL_SZ_EXPR + 1), TY_CHAR) + FL_NEXPR(fl) = nexpr + + return (nexpr) +end + + +# AT_FLNEXPR -- Get the next expression from an expression list. + +int procedure at_flnexpr (res, exprlist, ip, expr, maxch) + +pointer res #I pointer to the results descriptor +char exprlist[ARB] #I the input expression list +int ip #I pointer into the expression list +char expr[ARB] #O the output expression +int maxch #I maximum length of the output expression + +int ep, op, token, fnum +int ctotok(), strlen(), cq_fnumber(), ctoi(), cq_rstati(), cq_fname() + +begin + # Decode the column labels. + op = 1 + while (exprlist[ip] != EOS) { + + token = ctotok (exprlist, ip, expr[op], maxch) + if (expr[op] == EOS) + next + + if ((token == TOK_PUNCTUATION) && (expr[op] == ',')) { + if (op == 1) + next + else + break + } + + # Replace generic identifiers with their catalog equivalents. + if (token == TOK_IDENTIFIER) { + fnum = cq_fnumber (res, expr[op]) + if (fnum <= 0 && expr[op] == 'f') { + ep = 2 + if (ctoi (expr[op], ep, fnum) <= 0) + fnum = 0 + if (fnum < 1 || fnum > cq_rstati (res, CQRNRECS)) + fnum = 0 + if (fnum > 0) { + if (cq_fname (res, fnum, expr[op], maxch) != fnum) + ; + } + } + } + + + op = op + strlen (expr[op]) + } + + expr[op] = EOS + if ((exprlist[ip] == EOS) && (op == 1)) + return (EOF) + else + return (op - 1) +end + + +# AT_FLNITEM -- Get the next expression from an expression list. + +int procedure at_flnitem (itemlist, ip, item, maxch) + +char itemlist[ARB] #I the input item list +int ip #I pointer into the item list +char item[ARB] #O the output item +int maxch #I maximum length of the output item + +int op, token +int ctotok(), strlen() + +begin + # Decode the column labels. + op = 1 + while (itemlist[ip] != EOS) { + + token = ctotok (itemlist, ip, item[op], maxch) + if (item[op] == EOS) + next + + if ((token == TOK_PUNCTUATION) && (item[op] == ',')) { + if (op == 1) + next + else + break + } + + op = op + strlen (item[op]) + } + + item[op] = EOS + if ((itemlist[ip] == EOS) && (op == 1)) + return (EOF) + else + return (op - 1) +end + + +# AT_FLRANGES -- Get the field ranges for each output field. + +int procedure at_flranges (res, fl) + +pointer res #I the results descriptor +pointer fl #I the field list descriptor + +pointer sp, fname, fptr, rptr +int i, nin, nout, lindex, rindex, ip1, ip2, f1, f2 +char lbracket, rbracket +int cq_rstati(), strldx(), ctoi(), ctotok(), cq_fnumber() +bool streq() +data lbracket /'['/, rbracket /']'/ + +begin + # Allocate working space. + call smark (sp) + call salloc (fname, FL_SZ_EXPR, TY_CHAR) + + # Allocate space for the ranges list. + call calloc (FL_FRANGES(fl), 3 * FL_NEXPR(fl) + 1, TY_INT) + + # Initialize. + nin = cq_rstati (res, CQNFIELDS) + nout = 0 + + # Loop over the expressions. Fields which cannot be decoded + # have zero-valued range entries. Expression fields have INDEFI + # valued range entries. + fptr = FL_FLIST(fl) + rptr = FL_FRANGES(fl) + do i = 1, FL_NEXPR (fl) { + + lindex = strldx (lbracket, Memc[fptr]) + rindex = strldx (rbracket, Memc[fptr]) + ip1 = 1 + ip2 = lindex + 1 + + # Decode generic field ranges. + if (Memc[fptr] == 'f' && lindex == 2 && rindex > lindex) { + + # Find the range limits. + if (Memc[fptr+lindex] == '*') { + f1 = 1 + f2 = nin + } else { + if (ctoi (Memc[fptr], ip2, f1) <= 0) + f1 = 0 + else if (f1 < 1 || f1 > nin) + f1 = 0 + if (ctoi (Memc[fptr], ip2, f2) <= 0) + f2 = 0 + else + f2 = -f2 + if (f2 < 1 || f2 > nin) + f2 = 0 + } + + # Valid range. + if (f1 > 0 && f2 > f1) { + Memi[rptr] = f1 + Memi[rptr+1] = f2 + Memi[rptr+2] = 1 + nout = nout + f2 - f1 + 1 + + # Field cannot be decoded. + } else { + Memi[rptr] = 0 + Memi[rptr+1] = 0 + Memi[rptr+2] = 0 + } + + # Decode fields and expressions. + } else if (ctotok (Memc[fptr], ip1, Memc[fname], FL_SZ_EXPR) == + TOK_IDENTIFIER) { + + # Find the field number. + f1 = cq_fnumber (res, Memc[fptr]) + if (f1 <= 0 && streq (Memc[fptr], Memc[fname])) { + if (Memc[fptr] != 'f') + f1 = 0 + else { + f2 = 1 + if (ctoi (Memc[fptr+1], f2, f1) <= 0) + f1 = 0 + else if (f1 < 1 || f1 > nin) + f1 = 0 + } + } + + # Valid single field. + if (f1 > 0) { + Memi[rptr] = f1 + Memi[rptr+1] = f1 + Memi[rptr+2] = 1 + nout = nout + 1 + + # Field is an expression. + } else if (ctotok (Memc[fptr], ip1, Memc[fname], + FL_SZ_EXPR) != TOK_EOS) { + Memi[rptr] = INDEFI + Memi[rptr+1] = INDEFI + Memi[rptr+2] = INDEFI + nout = nout + 1 + + # Field cannot be decoded. + } else { + Memi[rptr] = 0 + Memi[rptr+1] = 0 + Memi[rptr+2] = 0 + } + + # What's left over is an expression field. + } else { + + Memi[rptr] = INDEFI + Memi[rptr+1] = INDEFI + Memi[rptr+2] = INDEFI + nout = nout + 1 + } + + fptr = fptr + FL_SZ_EXPR + 1 + rptr = rptr + 3 + } + + # Store the field counts. + FL_NFIELDS(fl) = nout + + call sfree (sp) + + return (nout) +end + + +# AT_FLFLIST -- Assemble the field header info for the new catalog. + +procedure at_flflist (at, res, fl) + +pointer at #I the astrometry package descriptor +pointer res #I the results descriptor +pointer fl #I the field list descriptor + +pointer sp, fnames, fntypes, fnunits, fnfmts, franame, fdecname, xpname +pointer ypname, xcname, ycname, str, rptr +int i, j, ip, nfields, fnp, ftp, fup, ffp, rtype, foffset, ival +int cq_rstati(), ctotok(), at_dtype(), ctoi(), cq_finfon() +int at_wrdstr(), at_stati(), at_flnitem() +bool streq(), strne() + +begin + # Get the number of output fields. + nfields = FL_NFIELDS(fl) + if (nfields <= 0) + return + + # Get some working space. + call smark (sp) + call salloc (fnames, SZ_LINE, TY_CHAR) + call salloc (fntypes, SZ_LINE, TY_CHAR) + call salloc (fnunits, SZ_LINE, TY_CHAR) + call salloc (fnfmts, SZ_LINE, TY_CHAR) + call salloc (franame, CQ_SZ_QPNAME, TY_CHAR) + call salloc (fdecname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (xpname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (ypname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (xcname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (ycname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the user parameters defining the names, types, units, and + # formats of the new fields. + call at_stats (at, FNAMES, Memc[fnames], SZ_LINE) + call at_stats (at, FNTYPES, Memc[fntypes], SZ_LINE) + call at_stats (at, FNUNITS, Memc[fnunits], SZ_LINE) + call at_stats (at, FNFORMATS, Memc[fnfmts], SZ_LINE) + fnp = 1 + ftp = 1 + fup = 1 + ffp = 1 + + # Get the special coordinate field names. + call at_stats (at, FIRA, Memc[franame], CQ_SZ_FNAME) + call at_stats (at, FIDEC, Memc[fdecname], CQ_SZ_FNAME) + call at_stats (at, FIXP, Memc[xpname], CQ_SZ_FNAME) + call at_stats (at, FIYP, Memc[ypname], CQ_SZ_FNAME) + call at_stats (at, FIXC, Memc[xcname], CQ_SZ_FNAME) + call at_stats (at, FIYC, Memc[ycname], CQ_SZ_FNAME) + + # Allocate space for the header field names, offsets, sizes, data + # types, units, and formats. + call calloc (FL_FNAMES(fl), nfields * (CQ_SZ_QPNAME + 1), TY_CHAR) + call calloc (FL_FOFFSETS(fl), nfields, TY_INT) + call calloc (FL_FSIZES(fl), nfields, TY_INT) + call calloc (FL_FTYPES(fl), nfields, TY_INT) + call calloc (FL_FUNITS(fl), nfields * (CQ_SZ_QPUNITS + 1), TY_CHAR) + call calloc (FL_FFMTS(fl), nfields * (CQ_SZ_QPFMTS + 1), TY_CHAR) + + # Get the output type. This is the same as the input type. + rtype = cq_rstati (res, CQRTYPE) + + # Loop over the ranges list. + nfields = 0 + foffset = 1 + rptr = FL_FRANGES(fl) + do i = 1, FL_NEXPR(fl) { + + # Skip non-decodable fields. + if (Memi[rptr] == 0) { + rptr = rptr + 3 + next + } + + # The field is an input catalog field. + if (! IS_INDEFI (Memi[rptr])) { + + do j = Memi[rptr], Memi[rptr+1] { + + # Get the field name, the field offset, size, data type, + # units and default format. + if (cq_finfon (res, j, Memc[FL_FNAMES(fl)+nfields* + (CQ_SZ_QPNAME+1)], CQ_SZ_QPNAME, Memi[FL_FOFFSETS(fl)+ + nfields], Memi[FL_FSIZES(fl)+nfields], + Memi[FL_FTYPES(fl)+nfields], + Memc[FL_FUNITS(fl)+nfields*(CQ_SZ_QPUNITS+1)], + CQ_SZ_QPUNITS, Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) != j) + ; + + # Correct the field offset and field size. + switch (rtype) { + case CQ_STEXT: + Memi[FL_FOFFSETS(fl)+nfields] = nfields + 1 + Memi[FL_FSIZES(fl)+nfields] = 0 + case CQ_BTEXT: + Memi[FL_FOFFSETS(fl)+nfields] = foffset + foffset = foffset + Memi[FL_FSIZES(fl)+nfields] + default: + call error (0, "Unknown output catalog type") + ; + } + + # Correct for coordinate units and format transform here. + if (streq (Memc[franame], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)])) { + ival = at_stati (at, FORAUNITS) + if (ival <= 0) + Memc[str] = EOS + else if (at_wrdstr (ival, Memc[str], SZ_FNAME, + AT_RA_UNITS) <= 0) + Memc[str] = EOS + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FUNITS(fl)+nfields*(CQ_SZ_QPUNITS+1)])) + call strcpy (Memc[str], Memc[FL_FUNITS(fl)+ + nfields*(CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) + call at_stats (at, FORAFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + if (streq (Memc[fdecname], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)])) { + ival = at_stati (at, FODECUNITS) + if (ival <= 0) + Memc[str] = EOS + else if (at_wrdstr (ival, Memc[str], SZ_FNAME, + AT_DEC_UNITS) <= 0) + Memc[str] = EOS + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FUNITS(fl)+nfields*(CQ_SZ_QPUNITS+1)])) + call strcpy (Memc[str], Memc[FL_FUNITS(fl)+ + nfields*(CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) + call at_stats (at, FODECFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + + # Correct for pixel coordinate formats here. + if (streq (Memc[xpname], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)]) || streq (Memc[xcname], + Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)])) { + call at_stats (at, FOXFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + if (streq (Memc[ypname], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)]) || streq (Memc[ycname], + Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)])) { + call at_stats (at, FOYFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + + + nfields = nfields + 1 + } + + # This field is a new field. + } else { + + # Get the field names. The default is f#. + ip = 1 + if (at_flnitem (Memc[fnames], fnp, Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)], CQ_SZ_QPNAME) == EOF) { + call sprintf (Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)], + CQ_SZ_QPNAME, "f%d") + call pargi (nfields + 1) + } else if (ctotok (Memc[FL_FNAMES(fl)+nfields* + (CQ_SZ_QPNAME+1)], ip, Memc[FL_FNAMES(fl)+nfields* + (CQ_SZ_QPNAME+1)], CQ_SZ_QPNAME) != TOK_IDENTIFIER) { + call sprintf (Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)], + CQ_SZ_QPNAME, "f%d") + call pargi (nfields + 1) + } + + # Get the data types. The default for now is type real. + ip = 1 + if (at_flnitem (Memc[fntypes], ftp, Memc[str], + SZ_FNAME) == EOF) { + Memi[FL_FTYPES(fl)+nfields] = TY_REAL + } else if (ctotok (Memc[str], ip, Memc[str], SZ_FNAME) != + TOK_IDENTIFIER) { + Memi[FL_FTYPES(fl)+nfields] = TY_REAL + } else { + #call strlwr (Memc[str]) + #Memi[FL_FTYPES(fl)+nfields] = cq_dtype(Memc[str]) + Memi[FL_FTYPES(fl)+nfields] = at_dtype(Memc[str]) + } + + # Get the data units. The default is INDEF. + ip = 1 + if (at_flnitem (Memc[fnunits], fup, Memc[FL_FUNITS(fl)+ + nfields*(CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) == EOF) { + call strcpy ("INDEF", Memc[FL_FUNITS(fl)+nfields* + (CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) + } else if (ctotok (Memc[FL_FUNITS(fl)+nfields* + (CQ_SZ_QPUNITS+1)], ip, Memc[FL_FUNITS(fl)+nfields* + (CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) != TOK_IDENTIFIER) { + call strcpy ("INDEF", Memc[FL_FUNITS(fl)+nfields* + (CQ_SZ_QPUNITS+1)], CQ_SZ_QPUNITS) + } + + # Get the data formats. The default is %10s, %10d, and %10g + # for character, integer, and floating data points respectively. + ip = 1 + if (at_flnitem (Memc[fnfmts], ffp, Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) == EOF) { + switch (Memi[FL_FTYPES(fl)+nfields]) { + case TY_CHAR: + call strcpy ("%10s", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + case TY_SHORT, TY_INT, TY_LONG: + call strcpy ("%10d", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + case TY_REAL, TY_DOUBLE: + call strcpy ("%10g", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + } else if (Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)] != '%') { + switch (Memi[FL_FTYPES(fl)+nfields*(CQ_SZ_QPFMTS+1)]) { + case TY_CHAR: + call strcpy ("%10s", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + case TY_SHORT, TY_INT, TY_LONG: + call strcpy ("%10d", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + case TY_REAL, TY_DOUBLE: + call strcpy ("%10g", Memc[FL_FFMTS(fl)+nfields* + (CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + } + + # Correct for pixel coordinate formats here. + if (streq (Memc[xpname], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)]) || streq (Memc[xcname], + Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)])) { + call at_stats (at, FOXFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + if (streq (Memc[ypname], Memc[FL_FNAMES(fl)+ + nfields*(CQ_SZ_QPNAME+1)]) || streq (Memc[ycname], + Memc[FL_FNAMES(fl)+nfields*(CQ_SZ_QPNAME+1)])) { + call at_stats (at, FOYFORMAT, Memc[str], SZ_FNAME) + if (Memc[str] != EOS && strne (Memc[str], + Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)])) + call strcpy (Memc[str], Memc[FL_FFMTS(fl)+ + nfields*(CQ_SZ_QPFMTS+1)], CQ_SZ_QPFMTS) + } + + # Get the field width. + ip = 2 + if (ctoi (Memc[FL_FFMTS(fl)+nfields*(CQ_SZ_QPFMTS+1)], ip, + ival) <= 0) + ival = 10 + else if (ival <= 0 || IS_INDEFI(ival)) + ival = 10 + + # Get the field offset and field size. Note the extra + # character added to the width ... + switch (rtype) { + case CQ_STEXT: + Memi[FL_FOFFSETS(fl)+nfields] = nfields + 1 + Memi[FL_FSIZES(fl)+nfields] = 0 + case CQ_BTEXT: + Memi[FL_FOFFSETS(fl)+nfields] = foffset + Memi[FL_FSIZES(fl)+nfields] = ival + 1 + foffset = foffset + Memi[FL_FSIZES(fl)+nfields] + default: + call error (0, "Unknown output catalog type") + } + + nfields = nfields + 1 + + } + + rptr = rptr + 3 + } + + call sfree (sp) +end + + +# AT_GETOP -- Fetch an operand from the data structure. + +procedure at_getop (res, operand, o) + +pointer res #I pointer to the data structure +char operand[ARB] #I name of operand to be returned +pointer o #I pointer to output operand + +pointer sp, fvalue +int fieldno, nchars +int cq_fnumber(), cq_ftype(), cq_gvald(), cq_gvali(), cq_gvalc() +int cq_rstati() + +begin + fieldno = cq_fnumber (res, operand) + if (fieldno <= 0) + call error (0, "Illegal operand in expression") + + switch (cq_ftype (res, operand)) { + + case TY_CHAR: + call smark (sp) + call salloc (fvalue, SZ_LINE, TY_CHAR) + nchars = cq_gvalc (res, cq_rstati(res, CQRECPTR), operand, + Memc[fvalue], SZ_LINE) + if (nchars <= 0) { + call strcpy ("INDEF", Memc[fvalue], 5) + nchars = 5 + } + O_TYPE(o) = TY_CHAR + O_LEN(o) = nchars + O_FLAGS(o) = O_FREEVAL + call malloc (O_VALP(o), nchars, TY_CHAR) + call strcpy (Memc[fvalue], O_VALC(o), nchars) + call sfree (sp) + + case TY_SHORT, TY_INT, TY_LONG: + O_TYPE(o) = TY_INT + O_LEN(o) = 0 + O_FLAGS(o) = 0 + nchars = cq_gvali (res, cq_rstati (res, CQRECPTR), operand, + O_VALI(o)) + + case TY_REAL, TY_DOUBLE: + O_TYPE(o) = TY_DOUBLE + O_LEN(o) = 0 + O_FLAGS(o) = 0 + nchars = cq_gvald (res, cq_rstati(res, CQRECPTR), operand, + O_VALD(o)) + + default: + call smark (sp) + call salloc (fvalue, SZ_LINE, TY_CHAR) + nchars = cq_gvalc (res, cq_rstati(res, CQRECPTR), operand, + Memc[fvalue], SZ_LINE) + if (nchars <= 0) { + call strcpy ("INDEF", Memc[fvalue], 5) + nchars = 5 + } + O_TYPE(o) = TY_CHAR + O_LEN(o) = nchars + O_FLAGS(o) = O_FREEVAL + call malloc (O_VALP(o), nchars, TY_CHAR) + call strcpy (Memc[fvalue], O_VALC(o), nchars) + call sfree (sp) + } +end diff --git a/noao/astcat/src/agetcat/athedit.x b/noao/astcat/src/agetcat/athedit.x new file mode 100644 index 00000000..73815db0 --- /dev/null +++ b/noao/astcat/src/agetcat/athedit.x @@ -0,0 +1,614 @@ +include "../../lib/astrom.h" +include "../../lib/aimpars.h" +include <pkg/cq.h> + + +# AT_HEDIT -- Add a set of standard keywords to the image header. + +procedure at_hedit (im, res, at, update, verbose) + +pointer im #I the input image descriptor +pointer res #I the image results descriptor +pointer at #I the astrometry package descriptor +bool update #I update the header ? +bool verbose #I verbose mode ? + +begin + if (res != NULL) + call at_dbkey (im, res, update, verbose) + if (at != NULL) + call at_parkey (im, at, update, verbose) +end + + +# AT_DBKEY -- Add a set of standard keywords required by astrometric +# reductions to the image header. New keywords will only be added if +# the keyword name is defined in the in the image survey database and the +# standard keyword does not already exist in the image header, or if the +# keyword has a default value in the image survey database. + +procedure at_dbkey (im, res, update, verbose) + +pointer im #I the input image descriptor +pointer res #I the image results descriptor +bool update #I update the header ? +bool verbose #I verbose mode ? + +pointer sp, kfield, kname, kvalue, kunits +int i, nkey, ktype +int cq_istati(), cq_kinfon(), imaccf(), at_akeyword() +bool streq() + +begin + call smark (sp) + call salloc (kfield, CQ_SZ_QPNAME, TY_CHAR) + call salloc (kname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (kvalue, CQ_SZ_QPVALUE, TY_CHAR) + call salloc (kunits, CQ_SZ_QPUNITS, TY_CHAR) + + # Loop over the keywords. + nkey = cq_istati (res, CQNIMPARS) + do i = 1, nkey { + + # Get the keyword information. + if (cq_kinfon (res, i, Memc[kfield], CQ_SZ_QPNAME, Memc[kname], + CQ_SZ_QPNAME, Memc[kvalue], CQ_SZ_QPVALUE, ktype, Memc[kunits], + CQ_SZ_QPUNITS) != i) + next + + # The keyword names is INDEF. + if (streq (Memc[kname], "INDEF")) { + + # Go to next keyword if the keyword value is also INDEF. + if (streq (Memc[kvalue], "INDEF")) + next + + # Add keyword with its default value if it does not exist. + if (imaccf (im, Memc[kfield]) == NO) { + if (at_akeyword (im, Memc[kfield], Memc[kvalue], ktype, + Memc[kunits], update) == OK) { + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding survey keyword %s = %s to header\n") + call pargstr (Memc[kfield]) + call pargstr (Memc[kvalue]) + } + #} else if (update || verbose) { + } else if (verbose) { + call printf ( + " Error adding survey keyword %s to header\n") + call pargstr (Memc[kfield]) + } + #} else if (update || verbose) { + } else if (verbose) { + call printf ( + " Warning survey keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + + # The keyword name is defined and it exists in the image. + } else if (imaccf (im, Memc[kname]) == YES) { + + # Add the new keyword with the old keyword value to the image + call imgstr (im, Memc[kname], Memc[kvalue], CQ_SZ_QPVALUE) + if (imaccf (im, Memc[kfield]) == NO) { + if (at_akeyword (im, Memc[kfield], Memc[kvalue], ktype, + Memc[kunits], update) == OK) { + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding survey keyword %s = %s to header\n") + call pargstr (Memc[kfield]) + call pargstr (Memc[kvalue]) + } + #} else if (update || verbose) { + } else if (verbose) { + call printf ( + " Error adding survey keyword %s to header\n") + call pargstr (Memc[kfield]) + } + #} else if (update || verbose) { + } else if (verbose) { + call printf ( + " Warning survey keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + + # The keywords names does not exist in the image. + #} else if (update || verbose) { + } else if (verbose) { + + call printf ( + " Warning survey keyword %s value %s does not exist\n") + call pargstr (Memc[kfield]) + call pargstr (Memc[kname]) + } + } + + call sfree (sp) +end + + +# AT_PARKEY -- Add a set of standard keywords required by astrometric +# reductions to the image header. New keywords will only be added if +# the keyword name is defined in default AIMPARS parameter file and the +# standard keyword does not already exist in the image header, or if the +# keyword has a default value in the AIMPARS parameter file. + +procedure at_parkey (im, at, update, verbose) + +pointer im #I the input image descriptor +pointer at #I the astrometry package descriptor +bool update #I update the header ? +bool verbose #I verbose mode ? + +double dval +real rval +pointer imst, sym, sp, kfield, kvalue +int i, key +double at_statd(), imgetd() +real at_statr(), imgetr() +pointer at_statp(), stfind() +int at_wrdstr(), imaccf() +bool streq() +errchk imgetd(), imgetr(), imgstr() + +begin + if (at_statp (at, PIMPARS) == NULL) + return + imst = at_statp (at, IMST) + if (imst == NULL) + return + + call smark (sp) + call salloc (kfield, SZ_FNAME, TY_CHAR) + call salloc (kvalue, SZ_FNAME, TY_CHAR) + + # Loop over the keywords. + do i = 1, AT_NIMFIELDS { + + # Get the parameter name. + key = at_wrdstr (i, Memc[kfield], SZ_FNAME, AT_IMFIELDS) + switch (key) { + case HDR_OBSERVAT: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + call at_stats (at, OBSERVAT, Memc[kvalue], SZ_FNAME) + else iferr (call imgstr (im, AT_IMSTKVAL(sym), + Memc[kvalue], SZ_FNAME)) + call at_stats (at, OBSERVAT, Memc[kvalue], SZ_FNAME) + } else + call at_stats (at, OBSERVAT, Memc[kvalue], SZ_FNAME) + + if (! streq (Memc[kvalue], "INDEF")) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imastr (im, Memc[kfield], Memc[kvalue]) + if (verbose) { + call printf ( + " Adding default keyword %s = %s to header\n") + call pargstr (Memc[kfield]) + call pargstr (Memc[kvalue]) + } + } + } + + case HDR_ESITELNG: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + dval = at_statd (at, ESITELNG) + else iferr (dval = imgetd (im, AT_IMSTKVAL(sym))) + dval = at_statd (at, ESITELNG) + } else + dval = at_statd (at, ESITELNG) + if (! IS_INDEFD(dval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddd (im, Memc[kfield], dval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %h to header\n") + call pargstr (Memc[kfield]) + call pargd (dval) + } + } + } + + case HDR_ESITELAT: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + dval = at_statd (at, ESITELAT) + else iferr (dval = imgetd (im, AT_IMSTKVAL(sym))) + dval = at_statd (at, ESITELAT) + } else + dval = at_statd (at, ESITELAT) + if (! IS_INDEFD(dval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddd (im, Memc[kfield], dval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %h to header\n") + call pargstr (Memc[kfield]) + call pargd (dval) + } + } + } + + case HDR_ESITEALT: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, ESITEALT) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, ESITEALT) + } else + rval = at_statr (at, ESITEALT) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %0.1f to header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_ESITETZ: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, ESITETZ) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, ESITETZ) + } else + rval = at_statr (at, ESITETZ) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + if (verbose) { + call printf ( + " Adding default keyword %s = %0.1f to header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_EMJDOBS: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + dval = at_statd (at, EMJDOBS) + else iferr (dval = imgetd (im, AT_IMSTKVAL(sym))) + dval = at_statd (at, EMJDOBS) + } else + dval = at_statd (at, EMJDOBS) + if (! IS_INDEFD(dval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddd (im, Memc[kfield], dval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %0.5f to header\n") + call pargstr (Memc[kfield]) + call pargd (dval) + } + } + } + + case HDR_EDATAMIN: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, EDATAMIN) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, EDATAMIN) + } else + rval = at_statr (at, EDATAMIN) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %g to header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_EDATAMAX: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, EDATAMAX) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, EDATAMAX) + } else + rval = at_statr (at, EDATAMAX) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %g to image\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_EGAIN: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, EGAIN) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, EGAIN) + } else + rval = at_statr (at, EGAIN) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %0.1f to image\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_ERDNOISE: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, ERDNOISE) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, ERDNOISE) + } else + rval = at_statr (at, ERDNOISE) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %0.1f to image\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_EWAVLEN: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, EWAVLEN) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, EWAVLEN) + } else + rval = at_statr (at, EWAVLEN) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %0.1f to header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_ETEMP: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, ETEMP) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, ETEMP) + } else + rval = at_statr (at, ETEMP) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s to image header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + case HDR_EPRESS: + sym = stfind (imst, Memc[kfield]) + if (sym != NULL) { + if (streq (AT_IMSTKVAL(sym), "INDEF")) + rval = at_statr (at, EPRESS) + else iferr (rval = imgetr (im, AT_IMSTKVAL(sym))) + rval = at_statr (at, EPRESS) + } else + rval = at_statr (at, EPRESS) + if (! IS_INDEFR(rval)) { + if (imaccf (im, Memc[kfield]) == YES) { + #if (update || verbose) { + if (verbose) { + call printf (" Keyword %s already exists\n") + call pargstr (Memc[kfield]) + } + } else { + if (update) + call imaddr (im, Memc[kfield], rval) + #if (update || verbose) { + if (verbose) { + call printf ( + " Adding default keyword %s = %g to header\n") + call pargstr (Memc[kfield]) + call pargr (rval) + } + } + } + + default: + ; + } + } + + call sfree (sp) +end + + +# AT_AKEYWORD -- Add a new keyword to the image header. Note that at present +# nothing is done with the units information although this may be used in the +# future. + +int procedure at_akeyword (im, kname, kvalue, ktype, kunits, update) + +pointer im #I the image descriptor +char kname[ARB] #I the image keyword name +char kvalue[ARB] #I the image keyword value +int ktype #I the image keyword data type +char kunits[ARB] #I the image keyword units (not used) +bool update #I actually update the header ? + +double dval +real rval +long lval +int ip, stat +int ctod(), ctor(), ctol() + +begin + stat = OK + + switch (ktype) { + + case TY_DOUBLE: + ip = 1 + if (ctod (kvalue, ip, dval) > 0) { + if (update) + call imaddd (im, kname, dval) + } else + stat = ERR + + case TY_REAL: + ip = 1 + if (ctor (kvalue, ip, rval) > 0) { + if (update) + call imaddr (im, kname, rval) + } else + stat = ERR + + case TY_LONG, TY_INT, TY_SHORT: + ip = 1 + if (ctol (kvalue, ip, lval) > 0) { + if (update) + call imaddl (im, kname, lval) + } else + stat = ERR + + default: + if (update) + call imastr (im, kname, kvalue) + } + + return (stat) +end diff --git a/noao/astcat/src/agetcat/atincat.x b/noao/astcat/src/agetcat/atincat.x new file mode 100644 index 00000000..1d0ce61e --- /dev/null +++ b/noao/astcat/src/agetcat/atincat.x @@ -0,0 +1,70 @@ +# AT_GAPARS -- Read in the algorithm parameters for the AGETCAT task. + +procedure at_gapars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Initialize the region parameters. + call at_grcpset ("aregpars", at) + + # Initialize the catalog filter / selection parameters. + call at_gfspset ("afiltpars", at) +end + + +# AT_FAPARS -- Read in the algorithm parameters for the AFILTCAT task. + +procedure at_fapars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Initialize the catalog filter / selection parameters. + call at_gfspset ("afiltpars", at) +end + + +# AT_GIAPARS -- Read in the algorithm parameters for the AGETIM task. + +procedure at_giapars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Initialize the region parameters. + call at_grcpset ("aregpars", at) + + # Initialize the default wcs parameters. + #call at_gwcpset ("awcspars", at) + + # Initialize the default image data parameters. + #call at_gimpset ("aimpars", at) +end + + +# AT_HAPARS -- Read in the algorithm parameters for the AHEDIT task. + +procedure at_hapars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Initialize the default wcs parameters. + call at_gwcpset ("awcspars", at) + + # Initialize the default image data parameters. + call at_gimpset ("aimpars", at) +end + + +# AT_IAPARS -- Read in the algorithm parameters for the AIMFIND task. + +procedure at_iapars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Initialize the catalog filter / selection parameters. + call at_gfspset ("afiltpars", at) +end diff --git a/noao/astcat/src/agetcat/atoutcat.x b/noao/astcat/src/agetcat/atoutcat.x new file mode 100644 index 00000000..a3215003 --- /dev/null +++ b/noao/astcat/src/agetcat/atoutcat.x @@ -0,0 +1,72 @@ + + +# AT_GPPARS -- Update the AGETCAT task algorithm parameter sets. + +procedure at_gppars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Update the region definition parameters. + call at_prcpset ("aregpars", at) + + # Update the catalog filtering parameters. + call at_pfspset ("afiltpars", at) +end + + +# AT_FPPARS -- Update the AFILTCAT task algorithm parameter sets. + +procedure at_fppars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Update the catalog filtering parameters. + call at_pfspset ("afiltpars", at) +end + + +# AT_GIPPARS -- Update the AGETIM task algorithm parameter sets. + +procedure at_gippars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Update the region definition parameters. + call at_prcpset ("aregpars", at) + + # Update the default wcs parameters. + #call at_pwcpset ("awcspars", at) + + # Update the default image data parameters. + #call at_pimpset ("aimpars", at) +end + + +# AT_HPPARS -- Update the AHEDIT task algorithm parameter sets. + +procedure at_hppars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Update the default wcs parameters. + call at_pwcpset ("awcspars", at) + + # Update the default image data parameters. + call at_pimpset ("aimpars", at) +end + + +# AT_IPPARS -- Update the AIMFIND task algorithm parameter sets. + +procedure at_ippars (at) + +pointer at #I the pointer to the main astrom structure + +begin + # Update the catalog filtering parameters. + call at_pfspset ("afiltpars", at) +end diff --git a/noao/astcat/src/agetcat/atrcquery.x b/noao/astcat/src/agetcat/atrcquery.x new file mode 100644 index 00000000..393f7b4c --- /dev/null +++ b/noao/astcat/src/agetcat/atrcquery.x @@ -0,0 +1,522 @@ +include <math.h> +include <pkg/cq.h> +include <pkg/skywcs.h> +include "../../lib/astrom.h" +include "../../lib/acatalog.h" + +# AT_RCQUERY -- Format the catalog query for the specified field using +# field data supplied by the user and stored in a symbol table and query +# information stored in the catalog database. + +int procedure at_rcquery (at, cq, fieldno) + +pointer at #I the astrometry pacakge descriptor +pointer cq #I the database descriptor +int fieldno #I the field number descriptor + +double ra, dec, width +pointer sp, qsystem, fsystem, qpname, qpvalue, qpunits, qpformats, raformats +pointer decformats, symbol, qcoo, fcoo, mw +int i, stat, parno, units, nqpars + +pointer stfind(), at_statp() +int sk_decwcs(), cq_nqpars(), cq_gqparn(), cq_sqpar(), strdic(), at_wrdstr() +int sk_stati() +bool streq() +errchk cq_fgwrd() + +begin + call smark (sp) + + # Allocate space for the coordinate system descriptions. + call salloc (qsystem, SZ_FNAME, TY_CHAR) + call salloc (fsystem, SZ_FNAME, TY_CHAR) + + # Fetch the field center symbol. + call sprintf (Memc[qsystem], SZ_FNAME, "%s%d") + call pargstr (DEF_RCST_ROOTNAME) + call pargi (fieldno) + symbol = stfind (at_statp(at, RCST), Memc[qsystem]) + if (symbol == NULL) { + call sfree (sp) + return (ERR) + } + + # Determine the query coordinate system. If the query coordinate system + # is undefined, set it to the current catalog coordinate system. If + # the catalog system is undefined set it to the global default. + iferr (call cq_fgwrd (cq, "qsystem", Memc[qsystem], SZ_FNAME)) { + iferr (call cq_fgwrd (cq, "csystem", Memc[qsystem], SZ_FNAME)) + call strcpy ("DEF_CATSYSTEM", Memc[qsystem], SZ_FNAME) + } + if (Memc[qsystem] == EOS || streq (Memc[qsystem], "INDEF")) + call strcpy ("DEF_CATSYSTEM", Memc[qsystem], SZ_FNAME) + + # Open the query coordinate system data structure. + stat = sk_decwcs (Memc[qsystem], mw, qcoo, NULL) + if (stat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + call sk_close (qcoo) + call sfree (sp) + return (ERR) + } + + # Determine the field center coordinate system. If the field center + # coordinate system is undefined, set it to the query coordinate + # system. + if (AT_RCSTSYSTEM(symbol) == EOS || streq (AT_RCSTSYSTEM(symbol), + "INDEF")) + call strcpy (Memc[qsystem], Memc[fsystem], SZ_FNAME) + else + call strcpy (AT_RCSTSYSTEM(symbol), Memc[fsystem], SZ_FNAME) + + # Open the field center coordinate system data structure. + stat = sk_decwcs (Memc[fsystem], mw, fcoo, NULL) + if (stat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + call sk_close (fcoo) + call sk_close (qcoo) + call sfree (sp) + return (ERR) + } + + # Allocate space for the query parameter description. + call salloc (qpname, CQ_SZ_QPNAME, TY_CHAR) + call salloc (qpvalue, CQ_SZ_QPVALUE, TY_CHAR) + call salloc (qpunits, CQ_SZ_QPUNITS, TY_CHAR) + call salloc (qpformats, CQ_SZ_QPFMTS, TY_CHAR) + call salloc (raformats, CQ_SZ_QPFMTS, TY_CHAR) + call salloc (decformats, CQ_SZ_QPFMTS, TY_CHAR) + + # Loop through the query parameter list encoding the non-coordinate + # system parameters. + nqpars = cq_nqpars (cq) + do i = 1, nqpars { + + # Get the query parameter description. + if (cq_gqparn (cq, i, Memc[qpname], CQ_SZ_QPNAME, Memc[qpvalue], + CQ_SZ_QPVALUE, Memc[qpunits], CQ_SZ_QPUNITS, Memc[qpformats], + CQ_SZ_QPFMTS) != i) + next + + parno = strdic (Memc[qpname], Memc[qpname], CQ_SZ_QPNAME, + AT_QRCFIELDS) + if (parno <= 0) + next + + # Field center right ascension. Set the units and save the format + # for later use since we cannot perform the coordinate + # transformation until both ra and dec units are decoded. + switch (parno) { + case AT_QRCRA: + units = strdic (Memc[qpunits], Memc[qpunits], CQ_SZ_QPUNITS, + SKY_LNG_UNITLIST) + if (units > 0) + call sk_seti (qcoo, S_NLNGUNITS, units) + switch (AT_RCSTRAUNITS(symbol)) { + case AT_DEGREES: + units = SKY_DEGREES + case AT_RADIANS: + units = SKY_RADIANS + case AT_HOURS: + units = SKY_HOURS + default: + units = sk_stati (fcoo, S_NLNGUNITS) + } + call sk_seti (fcoo, S_NLNGUNITS, units) + + call strcpy (Memc[qpformats], Memc[raformats], CQ_SZ_QPFMTS) + + # Field center declination. Set the units and save the format + # for later use since we cannot perform the coordinate + # transformation until both ra and dec units are decoded. + case AT_QRCDEC: + units = strdic (Memc[qpunits], Memc[qpunits], CQ_SZ_QPUNITS, + SKY_LAT_UNITLIST) + if (units > 0) + call sk_seti (qcoo, S_NLATUNITS, units) + switch (AT_RCSTDECUNITS(symbol)) { + case AT_DEGREES: + units = SKY_DEGREES + case AT_RADIANS: + units = SKY_RADIANS + default: + units = sk_stati (fcoo, S_NLATUNITS) + } + call sk_seti (fcoo, S_NLATUNITS, units) + + call strcpy (Memc[qpformats], Memc[decformats], CQ_SZ_QPFMTS) + + # Width. Input units are minutes. Output units are minutes or + # degrees. + case AT_QRCWIDTH: + width = max (AT_RCSTRAWIDTH(symbol), AT_RCSTDECWIDTH(symbol)) + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Radius. Input units are minutes. Output units are minutes or + # degrees. + case AT_QRCRADIUS: + width = max (AT_RCSTRAWIDTH(symbol), + AT_RCSTDECWIDTH(symbol)) / 2.0d0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Half width. Input units are minutes. Output units are minutes or + # degrees. + case AT_QRCHWIDTH: + width = max (AT_RCSTRAWIDTH(symbol), + AT_RCSTDECWIDTH(symbol)) / 2.0d0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Ra width. Input units are minutes. Output units are minutes or + # degrees. + case AT_QRCRAWIDTH: + width = AT_RCSTRAWIDTH(symbol) + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Dec width. Input units are minutes. Output units are minutes or + # degrees. + case AT_QRCDECWIDTH: + width = AT_RCSTDECWIDTH(symbol) + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Ra half width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCRAHWIDTH: + width = AT_RCSTRAWIDTH(symbol) / 2.0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Dec half width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCDECHWIDTH: + width = AT_RCSTDECWIDTH(symbol) / 2.0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # X width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCXWIDTH: + width = AT_RCSTRAWIDTH(symbol) + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Y width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCYWIDTH: + width = AT_RCSTDECWIDTH(symbol) + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # X half width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCXHWIDTH: + width = AT_RCSTRAWIDTH(symbol) / 2.0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + # Y half width. Input units are minutes. Output units are minutes + # or degrees. + case AT_QRCYHWIDTH: + width = AT_RCSTDECWIDTH(symbol) / 2.0 + if (streq (Memc[qpunits], "degrees")) + width = width / 60.0d0 + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[qpformats]) + call pargd (width) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) != i) + ; + + + } + + } + + # Transform the ra and dec from the field center coordinate system to + # the query coordinate system and reformat the query. + call sk_ultran (fcoo, qcoo, AT_RCSTRA(symbol), AT_RCSTDEC(symbol), ra, + dec, 1) + if (at_wrdstr (AT_QRCRA, Memc[qpname], CQ_SZ_QPNAME, + AT_QRCFIELDS) > 0) { + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[raformats]) + call pargd (ra) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) == 0) + ; + } + if (at_wrdstr (AT_QRCDEC, Memc[qpname], CQ_SZ_QPNAME, + AT_QRCFIELDS) > 0) { + call sprintf (Memc[qpvalue], CQ_SZ_QPVALUE, Memc[decformats]) + call pargd (dec) + if (cq_sqpar (cq, Memc[qpname], Memc[qpvalue]) == 0) + ; + } + + # Cleanup. + call sk_close (fcoo) + call sk_close (qcoo) + call sfree (sp) + + return (OK) +end + + +# AT_RCREGION -- Determine the region extraction parameters for the specified +# field using field data supplied by the user and stored in a symbol table and +# information stored in the catalog . + +int procedure at_rcregion (at, cres, fieldno, ra, dec, rawidth, decwidth) + +pointer at #I the astrometry pacakge descriptor +pointer cres #I the catalog results descriptor +int fieldno #I the field number descriptor +double ra #O the field center ra in degrees +double dec #O the field center dec in degrees +double rawidth #O the field ra width in degrees +double decwidth #O the field dec width in degrees + +pointer sp, qsystem, fsystem, raname, decname, raunits, decunits +pointer symbol, mw, qcoo, fcoo +int stat, units +pointer stfind(), at_statp() +int sk_decwcs(), cq_hinfo(), strdic(), sk_stati() +bool streq() +errchk at_stats() + +begin + call smark (sp) + + # Allocate space for the coordinate system descriptions. + call salloc (qsystem, SZ_FNAME, TY_CHAR) + call salloc (fsystem, SZ_FNAME, TY_CHAR) + call salloc (raname, CQ_SZ_FNAME, TY_CHAR) + call salloc (decname, CQ_SZ_FNAME, TY_CHAR) + call salloc (raunits, CQ_SZ_FNAME, TY_CHAR) + call salloc (decunits, CQ_SZ_FNAME, TY_CHAR) + + # Fetch the field center symbol. + call sprintf (Memc[qsystem], SZ_FNAME, "%s%d") + call pargstr (DEF_RCST_ROOTNAME) + call pargi (fieldno) + symbol = stfind (at_statp(at, RCST), Memc[qsystem]) + if (symbol == NULL) { + call sfree (sp) + return (ERR) + } + + # Set the query coordinate system to the catalog coordinate system. + # the catalog system is undefined set it to the global default. + if (cq_hinfo (cres, "csystem", Memc[qsystem], SZ_FNAME) <= 0) + Memc[qsystem] = EOS + if (Memc[qsystem] == EOS || streq (Memc[qsystem], "INDEF")) + call strcpy ("DEF_CATSYSTEM", Memc[qsystem], SZ_FNAME) + + # Open the query coordinate system data structure. + stat = sk_decwcs (Memc[qsystem], mw, qcoo, NULL) + if (stat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + call sk_close (qcoo) + call sfree (sp) + return (ERR) + } + + # Determine the field center coordinate system. If the field center + # coordinate system is undefined, set it to the query coordinate + # system. + if (AT_RCSTSYSTEM(symbol) == EOS || streq (AT_RCSTSYSTEM(symbol), + "INDEF")) + call strcpy (Memc[qsystem], Memc[fsystem], SZ_FNAME) + else + call strcpy (AT_RCSTSYSTEM(symbol), Memc[fsystem], SZ_FNAME) + + # Open the field center coordinate system data structure. + stat = sk_decwcs (Memc[fsystem], mw, fcoo, NULL) + if (stat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + call sk_close (fcoo) + call sk_close (qcoo) + call sfree (sp) + return (ERR) + } + + # Get the names of the columns containing ra and dec. + iferr (call at_stats (at, FIRA, Memc[raname], CQ_SZ_FNAME)) + call strcpy ("ra", Memc[raname], SZ_FNAME) + iferr (call at_stats (at, FIDEC, Memc[decname], CQ_SZ_FNAME)) + call strcpy ("dec", Memc[decname], SZ_FNAME) + + # Get the query ra units. + call cq_funits (cres, Memc[raname], Memc[raunits], CQ_SZ_FUNITS) + units = strdic (Memc[raunits], Memc[raunits], CQ_SZ_FUNITS, + SKY_LNG_UNITLIST) + if (units > 0) + call sk_seti (qcoo, S_NLNGUNITS, units) + else + units = sk_stati (qcoo, S_NLNGUNITS) + switch (AT_RCSTRAUNITS(symbol)) { + case AT_DEGREES: + units = SKY_DEGREES + case AT_RADIANS: + units = SKY_RADIANS + case AT_HOURS: + units = SKY_HOURS + default: + ; + } + call sk_seti (fcoo, S_NLNGUNITS, units) + + # Get the query dec units. + call cq_funits (cres, Memc[decname], Memc[decunits], CQ_SZ_FUNITS) + units = strdic (Memc[decunits], Memc[decunits], CQ_SZ_FUNITS, + SKY_LAT_UNITLIST) + if (units > 0) + call sk_seti (qcoo, S_NLATUNITS, units) + else + units = sk_stati (qcoo, S_NLATUNITS) + switch (AT_RCSTDECUNITS(symbol)) { + case AT_DEGREES: + units = SKY_DEGREES + case AT_RADIANS: + units = SKY_RADIANS + case AT_HOURS: + units = SKY_HOURS + default: + ; + } + call sk_seti (fcoo, S_NLATUNITS, units) + + # Transform the ra and dec from the field center coordinate system to + # the query coordinate system and convert the units to degrees. + call sk_ultran (fcoo, qcoo, AT_RCSTRA(symbol), AT_RCSTDEC(symbol), ra, + dec, 1) + + # Transform the ra, dec, and width parameters to degrees. + switch (sk_stati(qcoo, S_NLNGUNITS)) { + case SKY_HOURS: + ra = 15.0d0 * ra + case SKY_DEGREES: + ; + case SKY_RADIANS: + ra = DRADTODEG (ra) + default: + ; + } + switch (sk_stati(qcoo, S_NLATUNITS)) { + case SKY_HOURS: + dec = 15.0d0 * dec + case SKY_DEGREES: + ; + case SKY_RADIANS: + dec = DRADTODEG (dec) + default: + ; + } + rawidth = AT_RCSTRAWIDTH(symbol) / 60.0d0 + decwidth = AT_RCSTDECWIDTH(symbol) / 60.0d0 + + # Cleanup. + call sk_close (fcoo) + call sk_close (qcoo) + call sfree (sp) + + return (OK) +end + + +# AT_RCLIMITS -- Given the ra, dec, ra width, and dec width of the field +# compute the field corners. + +procedure at_rclimits (ra, dec, rawidth, decwidth, ra1, ra2, dec1, dec2) + +double ra #I the field center ra in degrees +double dec #I the field center dec in degrees +double rawidth #I the field ra width in degrees +double decwidth #I the field dec width in degrees +double ra1 #O lower ra limit in degrees +double ra2 #O upper ra limit in degrees +double dec1 #O lower dec limit in degrees +double dec2 #O upper dec limit in degrees + +double cosdec, dra + +begin + # Find the field corners. + dec1 = dec - 0.5d0 * decwidth + dec2 = dec + 0.5d0 * decwidth + if (dec1 <= -90.0d0) { + dec1 = -90.0d0 + dec2 = min (dec + 0.5d0 * decwidth, 90.0d0) + ra1 = 0.0d0 + ra2 = 360.0d0 + return + } else if (dec2 >= 90.0d0) { + dec2 = 90.0d0 + dec1 = max (dec - 0.5d0 * decwidth, -90.0d0) + ra1 = 0.0d0 + ra2 = 360.0d0 + } else { + if (dec > 0.0d0) + cosdec = cos (DEGTORAD(dec2)) + else + cosdec = cos (DEGTORAD(dec1)) + dra = 0.50d0 * rawidth / cosdec + if (dra >= 180.0d0) { + ra1 = 0.0d0 + ra2 = 360.0d0 + } else { + ra1 = ra - dra + if (ra1 < 0.0d0) + ra1 = ra1 + 360.0d0 + ra2 = ra + dra + if (ra2 > 360.0d0) + ra2 = ra2 - 360.0d0 + } + } + +end diff --git a/noao/astcat/src/agetcat/atrcrd.x b/noao/astcat/src/agetcat/atrcrd.x new file mode 100644 index 00000000..7a770c4f --- /dev/null +++ b/noao/astcat/src/agetcat/atrcrd.x @@ -0,0 +1,314 @@ +include <fset.h> +include <imhdr.h> +include <mwset.h> +include <pkg/skywcs.h> +include "../../lib/astrom.h" +include "../../lib/acatalog.h" + +# AT_RCLIST -- Create a list of field centers. + +int procedure at_rclist (at, rcsource) + +pointer at #I the astrometry descriptor +char rcsource[ARB] #I the source of the regions list + +pointer sp, symname, st, sym +int nfields, fd, imlist +double at_statd() +pointer at_statp(), stopen(), stenter() +int at_stati(), access(), open(), imtopen(), at_rcread(), at_rcwcsim() +bool streq() + +begin + # Store the rcsource name in the data structure. + call at_sets (at, RCSOURCE, rcsource) + + # Check that the field center pointer is defined. + if (at_statp (at, PRCENTER) == NULL) + return (0) + + call smark (sp) + call salloc (symname, SZ_FNAME, TY_CHAR) + + # Open the symbol table. + if (at_statp (at, RCST) != NULL) + call stclose (at_statp(at, RCST)) + st = stopen ("fclist", 2 * DEF_LEN_RCST, DEF_LEN_RCST, + 10 * DEF_LEN_RCST) + call at_setp (at, RCST, st) + + # The source is the fcpars parameter set. + if (streq (rcsource, "pars")) { + + if (at_statd (at, RCRA) < 0.0d0 || at_statd(at, RCRA) > 360.0d0) + nfields = 0 + else if (at_statd (at, RCDEC) < -90.0d0 || at_statd(at, RCDEC) > + 90.0d0) + nfields = 0 + else if (at_statd (at, RCRAWIDTH) / 60.0d0 <= 0.0d0 || + at_statd(at, RCRAWIDTH) / 60.0d0 > 360.0d0) + nfields = 0 + else if (at_statd (at, RCDECWIDTH) / 60.0d0 <= 0.0d0 || + at_statd(at, RCDECWIDTH) / 60.0d0 > 180.0d0) + nfields = 0 + else { + call sprintf (Memc[symname], SZ_FNAME, "%s1") + call pargstr (DEF_RCST_ROOTNAME) + sym = stenter (st, Memc[symname], LEN_RCST_STRUCT) + AT_RCSTRA(sym) = at_statd (at, RCRA) + AT_RCSTDEC(sym) = at_statd (at, RCDEC) + AT_RCSTRAWIDTH(sym) = at_statd (at, RCRAWIDTH) + AT_RCSTDECWIDTH(sym) = at_statd (at, RCDECWIDTH) + AT_RCSTRAUNITS(sym) = at_stati (at, RCRAUNITS) + AT_RCSTDECUNITS(sym) = at_stati (at, RCDECUNITS) + call at_stats (at, RCSYSTEM, Memc[symname], SZ_FNAME) + call strcpy (Memc[symname], AT_RCSTSYSTEM(sym), SZ_FNAME) + call strcpy ("pars", AT_RCSTSOURCE(sym), SZ_FNAME) + call strcpy ("", AT_RCSTNAME(sym), SZ_FNAME) + nfields = 1 + } + + # The source is a text file. + } else if (access (rcsource, READ_ONLY, TEXT_FILE) == YES) { + + fd = open (rcsource, READ_ONLY, TEXT_FILE) + nfields = at_rcread (fd, at, st) + call close (fd) + + # The field center source is a list of images. Assume for now that + # images with celestial coordinate systems have a wcs system name + # of "image". This is true of images with a standard FITS wcs and + # for images with a wcs created by the core IRAF tasks. + } else { + imlist = imtopen (rcsource) + nfields = at_rcwcsim (imlist, at, st) + call imtclose (imlist) + } + + call sfree (sp) + + return (nfields) +end + + +# AT_RCREAD -- Read in the field center information from a text file. + +int procedure at_rcread (fd, at, st) + +int fd #I the field center file descriptor +pointer at #I the astrometry descriptor +pointer st #I the field center symbol table descriptor. + +double ra, dec, rawidth, decwidth +pointer sp, symname, sym +int nfields +pointer stenter() +int fscan(), nscan(), at_stati(), strdic() + +begin + call smark (sp) + call salloc (symname, SZ_FNAME, TY_CHAR) + + nfields = 0 + while (fscan(fd) != EOF) { + + # Get the minimum number of fields. + call gargd (ra) + call gargd (dec) + call gargd (rawidth) + call gargd (decwidth) + if (nscan() < 4) + next + if (ra < 0.0d0 || ra > 360.0d0) + next + if (dec < -90.0d0 || dec > 90.0d0) + next + if (rawidth / 60.0d0 <= 0.0d0 || rawidth / 60.0d0 > 360.0d0) + next + if (decwidth / 60.0d0 <= 0.0d0 || decwidth / 60.0d0 > 180.0d0) + next + + # Get the next symbols. + nfields = nfields + 1 + call sprintf (Memc[symname], SZ_FNAME, "%s%d") + call pargstr (DEF_RCST_ROOTNAME) + call pargi (nfields) + sym = stenter (st, Memc[symname], LEN_RCST_STRUCT) + + AT_RCSTRA(sym) = ra + AT_RCSTDEC(sym) = dec + AT_RCSTRAWIDTH(sym) = rawidth + AT_RCSTDECWIDTH(sym) = decwidth + + # Set the source and source name. + call strcpy ("file", AT_RCSTSOURCE(sym), SZ_FNAME) + call fstats (fd, F_FILENAME, Memc[symname], SZ_FNAME) + call strcpy (Memc[symname], AT_RCSTNAME(sym), SZ_FNAME) + + # Decode the units. + call gargwrd (Memc[symname], SZ_FNAME) + if (nscan() < 5) { + AT_RCSTRAUNITS(sym) = at_stati (at, RCRAUNITS) + AT_RCSTDECUNITS(sym) = at_stati (at, RCDECUNITS) + call at_stats (at, RCSYSTEM, Memc[symname], SZ_FNAME) + call strcpy (Memc[symname], AT_RCSTSYSTEM(sym), SZ_FNAME) + next + } else + AT_RCSTRAUNITS(sym) = strdic (Memc[symname], Memc[symname], + SZ_FNAME, AT_RA_UNITS) + call gargwrd (Memc[symname], SZ_FNAME) + if (nscan() < 6) { + AT_RCSTDECUNITS(sym) = at_stati (at, RCDECUNITS) + call at_stats (at, RCSYSTEM, Memc[symname], SZ_FNAME) + call strcpy (Memc[symname], AT_RCSTSYSTEM(sym), SZ_FNAME) + next + } else + AT_RCSTDECUNITS(sym) = strdic (Memc[symname], Memc[symname], + SZ_FNAME, AT_DEC_UNITS) + + # Decode the coordinate system. + call gargstr (Memc[symname], SZ_FNAME) + if (Memc[symname] == EOS || nscan() < 7) { + call at_stats (at, RCSYSTEM, Memc[symname], SZ_FNAME) + call strcpy (Memc[symname], AT_RCSTSYSTEM(sym), SZ_FNAME) + } else + call strcpy (Memc[symname], AT_RCSTSYSTEM(sym), SZ_FNAME) + + } + + call sfree (sp) + + return (nfields) +end + + +# AT_RCWCSIM -- Read in the field center information from a list of images. + +int procedure at_rcwcsim (imlist, at, st) + +int imlist #I the image list descriptor +pointer at #I the astrometry descriptor +pointer st #I the field center symbol table descriptor. + +double ra, dec, width +pointer sp, image, symname, im, mw, coo, sym, ct +int nfields +pointer immap(), mw_sctran(), stenter() +int imtgetim(), sk_decim(), sk_stati() + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (symname, SZ_FNAME, TY_CHAR) + + nfields = 0 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # The image must be 2D. + im = immap (Memc[image], READ_ONLY, 0) + if (IM_NDIM(im) != 2) { + call imunmap (im) + next + } + + # The image must have a FITS celestial coordinate system. + if (sk_decim (im, "world", mw, coo) == ERR) { + if (mw != NULL) + call mw_close (mw) + call sk_close (coo) + call imunmap (im) + next + } + + # Find the center of the image. + ct = mw_sctran (mw, "logical", "world", 03B) + call mw_c2trand (ct, double((1.0d0 + IM_LEN(im,1)) / 2.0d0), + double((1.0d0 + IM_LEN(im,2)) / 2.0d0), ra, dec) + if (ra < 0.0d0 || ra > 360.0d0) + next + if (dec < -90.0d0 || dec > 90.0d0) + next + + # Find the width of the field. + call at_gfwidth (im, mw, sk_stati(coo, S_PLNGAX), + sk_stati(coo, S_PLATAX), width) + + # Get the next symol. + nfields = nfields + 1 + call sprintf (Memc[symname], SZ_FNAME, "%s%d") + call pargstr (DEF_RCST_ROOTNAME) + call pargi (nfields) + sym = stenter (st, Memc[symname], LEN_RCST_STRUCT) + + AT_RCSTRA(sym) = ra + AT_RCSTDEC(sym) = dec + AT_RCSTRAWIDTH(sym) = width + AT_RCSTDECWIDTH(sym) = width + AT_RCSTRAUNITS(sym) = AT_DEGREES + AT_RCSTDECUNITS(sym) = AT_DEGREES + call sk_enwcs (coo, AT_RCSTSYSTEM(sym), SZ_FNAME) + + call strcpy ("image", AT_RCSTSOURCE(sym), SZ_FNAME) + call strcpy (Memc[image], AT_RCSTNAME(sym), SZ_FNAME) + + # Cleanup. + call sk_close (coo) + call mw_close (mw) + call imunmap (im) + } + + call sfree (sp) + + return (nfields) +end + + +define NEWCD Memd[ncd+(($2)-1)*ndim+($1)-1] + +# AT_GFWIDTH -- Estimate the field width in arcminutes from the size of the +# image and the image wcs. + +procedure at_gfwidth (im, mw, lngax, latax, width) + +pointer im #I the input image desciptor +pointer mw #I the input wcs descriptor +int lngax #I the longitude axis +int latax #I the latitude axis +double width #O the output field width in minutes of arc + +double scale +pointer r, cd, ltm, iltm, ncd +int ndim +int mw_stati() + +begin + # Get the dimension of the wcs. + ndim = mw_stati (mw, MW_NPHYSDIM) + + # Allocate working memory. + call malloc (r, ndim * ndim, TY_DOUBLE) + call malloc (cd, ndim * ndim, TY_DOUBLE) + call malloc (ltm, ndim * ndim, TY_DOUBLE) + call malloc (iltm, ndim * ndim, TY_DOUBLE) + call malloc (ncd, ndim * ndim, TY_DOUBLE) + + # Compute the original world to logical transformation. + call mw_gwtermd (mw, Memd[r], Memd[r], Memd[cd], ndim) + call mw_gltermd (mw, Memd[ltm], Memd[r], ndim) + call mwinvertd (Memd[ltm], Memd[iltm], ndim) + call mwmmuld (Memd[cd], Memd[iltm], Memd[ncd], ndim) + + # Estimate the scale. + scale = max (sqrt (NEWCD(lngax,lngax)**2 + NEWCD(lngax,latax)**2), + sqrt (NEWCD(latax,lngax)**2 + NEWCD(latax,latax)**2)) + + # Compute the width + width = 60.0d0 * scale * max (IM_LEN(im,1), IM_LEN(im,2)) + + # Free the space. + call mfree (r, TY_DOUBLE) + call mfree (cd, TY_DOUBLE) + call mfree (ncd, TY_DOUBLE) + call mfree (ltm, TY_DOUBLE) + call mfree (iltm, TY_DOUBLE) +end diff --git a/noao/astcat/src/agetcat/atrcsym.x b/noao/astcat/src/agetcat/atrcsym.x new file mode 100644 index 00000000..268ea86b --- /dev/null +++ b/noao/astcat/src/agetcat/atrcsym.x @@ -0,0 +1,29 @@ +include "../../lib/astrom.h" + +# AT_RCSYM -- Return the symbol for the specified field nymber. + +pointer procedure at_rcsym (at, fieldno) + +pointer at #I the astrometry package descriptor +int fieldno #I the region whose symbol is to be locate + +pointer sp, symname, st, sym +pointer at_statp(), stfind() + +begin + st = at_statp (at, RCST) + if (st == NULL) + return (NULL) + + call smark (sp) + call salloc (symname, SZ_FNAME, TY_CHAR) + + call sprintf (Memc[symname], SZ_FNAME, "%s%d") + call pargstr (DEF_RCST_ROOTNAME) + call pargi (fieldno) + sym = stfind (st, Memc[symname]) + + call sfree (sp) + + return (sym) +end diff --git a/noao/astcat/src/agetcat/attquery.x b/noao/astcat/src/agetcat/attquery.x new file mode 100644 index 00000000..3f8ff6d4 --- /dev/null +++ b/noao/astcat/src/agetcat/attquery.x @@ -0,0 +1,183 @@ +include <math.h> +include <pkg/cq.h> +include <pkg/skywcs.h> +include "../../lib/astrom.h" + +# AT_TQUERY -- Extract catalog objects from a text file stored in a results +# query structure. + +pointer procedure at_tquery (at, cq, cres, hdrtext, nlines, fieldno) + +pointer at #I the astrometry package descriptor +pointer cq #I the astrometric catalog descriptor +pointer cres #I the input catalog results descriptor +char hdrtext[ARB] #I the catalog header test +int nlines #I the number of lines in the header text +int fieldno #I the region number + +double rac, decc, ra, dec, rawidth, decwidth, ra1, ra2, dec1, dec2, dist +double tra, trac +pointer sp, csystem, raname, decname, funits, tmpname, line +pointer res, ccoo, mw +int i, fd, strfd, stat, units +pointer cq_fquery() +int cq_rstati(), at_rcregion(), open(), strlen(), cq_hinfo(), sk_decwcs() +int stropen(), getline(), strdic(), cq_gvald(), cq_grecord(), sk_stati() +bool streq() + +begin + # Return if the input catalog is undefined or contains no records. + if (cres == NULL) + return (NULL) + if (cq_rstati (cres, CQRNRECS) <= 0) + return (NULL) + + # Return if the header is undefined. + if (nlines <= 0 || hdrtext[1] == EOS) + return (NULL) + + # Get the region to be extracted. + if (at_rcregion (at, cres, fieldno, rac, decc, rawidth, + decwidth) == ERR) + return (NULL) + + # Compute the ra and dec limits. + call at_rclimits (rac, decc, rawidth, decwidth, ra1, ra2, dec1, dec2) + + # Get some working space. + call smark (sp) + call salloc (csystem, CQ_SZ_FNAME, TY_CHAR) + call salloc (tmpname, SZ_FNAME, TY_CHAR) + call salloc (raname, CQ_SZ_FNAME, TY_CHAR) + call salloc (decname, CQ_SZ_FNAME, TY_CHAR) + call salloc (funits, CQ_SZ_FUNITS, TY_CHAR) + call salloc (line, SZ_LINE, TY_CHAR) + + # Open the catalog coordinate system. + if (cq_hinfo (cres, "csystem", Memc[csystem], SZ_FNAME) <= 0) + Memc[csystem] = EOS + if (Memc[csystem] == EOS || streq (Memc[csystem], "INDEF")) + call strcpy ("DEF_CATSYSTEM", Memc[csystem], SZ_FNAME) + + # Open the query coordinate system data structure. + stat = sk_decwcs (Memc[csystem], mw, ccoo, NULL) + if (stat == ERR || mw != NULL) { + if (mw != NULL) + call mw_close (mw) + call sk_close (ccoo) + call sfree (sp) + return (NULL) + } + + # Open the temporary results file. + call mktemp ("res", Memc[tmpname], SZ_FNAME) + fd = open (Memc[tmpname], NEW_FILE, TEXT_FILE) + + # Write the file header to the temporary results file. + strfd = stropen (hdrtext, strlen(hdrtext), READ_ONLY) + call fprintf (fd, "# BEGIN CATALOG HEADER\n") + while (getline (strfd, Memc[line]) != EOF) { + call fprintf (fd, "# %s") + call pargstr (Memc[line]) + } + call fprintf (fd, "# END CATALOG HEADER\n#\n") + call strclose (strfd) + + # Determine the names of the ra and dec columns. + iferr (call at_stats (at, FIRA, Memc[raname], CQ_SZ_FNAME)) + call strcpy ("ra", Memc[raname], CQ_SZ_FNAME) + iferr (call at_stats (at, FIDEC, Memc[decname], CQ_SZ_FNAME)) + call strcpy ("dec", Memc[decname], CQ_SZ_FNAME) + + # Determine the units of the ra and dec keywords. + call cq_funits (cres, Memc[raname], Memc[funits], CQ_SZ_QPUNITS) + units = strdic (Memc[funits], Memc[funits], CQ_SZ_FUNITS, + SKY_LNG_UNITLIST) + if (units > 0) + call sk_seti (ccoo, S_NLNGUNITS, units) + call cq_funits (cres, Memc[decname], Memc[funits], CQ_SZ_QPUNITS) + units = strdic (Memc[funits], Memc[funits], CQ_SZ_FUNITS, + SKY_LAT_UNITLIST) + if (units > 0) + call sk_seti (ccoo, S_NLATUNITS, units) + + # Loop over the catalog records selecting those that match + # the region description. + do i = 1, cq_rstati (cres, CQRNRECS) { + + # Decode the coordinates. + if (cq_gvald (cres, i, Memc[raname], ra) <= 0) + next + if (cq_gvald (cres, i, Memc[decname], dec) <= 0) + next + + # Determine the coordinate units. + switch (sk_stati(ccoo, S_NLNGUNITS)) { + case SKY_HOURS: + ra = 15.0d0 * ra + case SKY_DEGREES: + ; + case SKY_RADIANS: + ra = DRADTODEG(ra) + default: + ; + } + switch (sk_stati(ccoo, S_NLATUNITS)) { + case SKY_HOURS: + dec = 15.0d0 * dec + case SKY_DEGREES: + ; + case SKY_RADIANS: + dec = DRADTODEG(dec) + default: + ; + } + + # Test the limits + if (dec < dec1 || dec > dec2) + next + if (ra1 < ra2) { + if (ra < ra1 || ra > ra2) + next + } else { + if (ra > ra2 && ra < ra1) + next + } + + # Check the longitude coordinate distance to remove pathologies + # in longitude or latitude strips involving the pole. This is + # an extra test of my own. + if (ra1 < ra2) { + dist = abs (ra - rac) + } else { + if (ra > ra1) + tra = ra - 360.0d0 + else + tra = ra + if (rac > ra1) + trac = rac - 360.0d0 + else + trac = rac + dist = abs (tra - trac) + } + if (abs (2.0d0 * dist *cos(DEGTORAD(dec))) > rawidth) + next + + # Record has been selected. + if (cq_grecord (cres, Memc[line], SZ_LINE, i) <= 0) + next + call putline (fd, Memc[line]) + } + + # Close the tmeporary file. + call close (fd) + + # Query the temporary file and then delete it. + res = cq_fquery (cq, Memc[tmpname], hdrtext) + call delete (Memc[tmpname]) + + # Clean up. + call sfree (sp) + + return (res) +end diff --git a/noao/astcat/src/agetcat/atwcat.x b/noao/astcat/src/agetcat/atwcat.x new file mode 100644 index 00000000..01d9ef6b --- /dev/null +++ b/noao/astcat/src/agetcat/atwcat.x @@ -0,0 +1,197 @@ +include <pkg/cq.h> + +# AT_WNOFILRECS -- Write out the catalog header and records without filtering. + +procedure at_wnofilrecs (fd, res, standard) + +int fd #I the output file descriptor +pointer res #I the results descriptor +bool standard #I write a standard catalog header. + +int nlines, nrecs +int at_wcathdr(), at_wcatrecs() + +begin + # Write out the catalog header. + if (standard) + nlines = at_wcathdr (fd, res) + + # Write out the records. + nrecs = at_wcatrecs (fd, res) +end + + +# AT_WCATHDR -- Write out a catalog header. + +int procedure at_wcathdr (fd, res) + +int fd #I the output file descriptor +pointer res #I the results descriptor + +pointer sp, catname, qpnames, qpvalues, qpunits, fname, fvalue, funits, ffmts +int i, nlines, nfields, fsize, foffset, ftype +int at_wrdstr(), cq_rstati(), cq_hinfon(), cq_finfon() +char cq_itype() + +begin + nlines = 0 + + # Allocate working space. + call smark (sp) + call salloc (catname, SZ_FNAME, TY_CHAR) + call salloc (fname, max (CQ_SZ_QPNAME, CQ_SZ_FNAME), TY_CHAR) + call salloc (fvalue, CQ_SZ_QPVALUE, TY_CHAR) + call salloc (funits, max (CQ_SZ_QPUNITS, CQ_SZ_FUNITS), TY_CHAR) + call salloc (ffmts, CQ_SZ_FFMTS, TY_CHAR) + call salloc (qpnames, SZ_LINE, TY_CHAR) + call salloc (qpvalues, SZ_LINE, TY_CHAR) + call salloc (qpunits, SZ_LINE, TY_CHAR) + + # Write the header banner. + call fprintf (fd, "# BEGIN CATALOG HEADER\n") + nlines = nlines + 1 + + # Write the catalog database and id. + call cq_rstats (res, CQRCATDB, Memc[catname], SZ_FNAME) + call fprintf (fd, "# catdb %s\n") + call pargstr (Memc[catname]) + nlines = nlines + 1 + call cq_rstats (res, CQRCATNAME, Memc[catname], SZ_FNAME) + call fprintf (fd, "# catname %s\n") + call pargstr (Memc[catname]) + nlines = nlines + 1 + + # Write out the query parameter names, values, and units used + # to generate the catalog. + call cq_rstats (res, CQRQPNAMES, Memc[qpnames], SZ_LINE) + call cq_rstats (res, CQRQPVALUES, Memc[qpvalues], SZ_LINE) + call cq_rstats (res, CQRQPUNITS, Memc[qpunits], SZ_LINE) + nfields = cq_rstati (res, CQRNQPARS) + call fprintf (fd, "# nquery %d\n") + call pargi (nfields) + nlines = nlines + 1 + do i = 1, nfields { + if (at_wrdstr (i, Memc[fname], CQ_SZ_QPNAME, Memc[qpnames]) != i) + ; + if (at_wrdstr (i, Memc[fvalue], CQ_SZ_QPVALUE, Memc[qpvalues]) != i) + ; + if (at_wrdstr (i, Memc[funits], CQ_SZ_QPUNITS, Memc[qpunits]) != i) + ; + call fprintf (fd, "# %s %s %s\n") + call pargstr (Memc[fname]) + call pargstr (Memc[fvalue]) + call pargstr (Memc[funits]) + nlines = nlines + 1 + } + + # Write out the results format type. + if (at_wrdstr (cq_rstati(res, CQRTYPE), Memc[fvalue], CQ_SZ_QPVALUE, + CQ_RTYPESTR) <= 0) + call strcpy ("stext", Memc[fvalue], CQ_SZ_QPVALUE) + call fprintf (fd, "# type %s\n") + call pargstr (Memc[fvalue]) + nlines = nlines + 1 + + # Write out the header parameters, + nfields = cq_rstati (res, CQNHEADER) + call fprintf (fd, "# nheader %d\n") + call pargi (nfields) + nlines = nlines + 1 + do i = 1, nfields { + if (cq_hinfon (res, i, Memc[fname], CQ_SZ_QPNAME, Memc[fvalue], + CQ_SZ_QPVALUE) != i) + next + call fprintf (fd, "# %s %s\n") + call pargstr (Memc[fname]) + call pargstr (Memc[fvalue]) + nlines = nlines + 1 + } + + # Write out the field parameters. + nfields = cq_rstati (res, CQNFIELDS) + call fprintf (fd, "# nfields %d\n") + call pargi (nfields) + do i = 1, nfields { + if (cq_finfon (res, i, Memc[fname], CQ_SZ_FNAME, foffset, fsize, + ftype, Memc[funits], CQ_SZ_FUNITS, Memc[ffmts], + CQ_SZ_FFMTS) != i) + next + call fprintf (fd, "# %s %d %d %c %s %s\n") + call pargstr (Memc[fname]) + call pargi (foffset) + call pargi (fsize) + call pargc (cq_itype (ftype)) + call pargstr (Memc[funits]) + call pargstr (Memc[ffmts]) + nlines = nlines + 1 + } + + # Write the header trailer. + call fprintf (fd, "# END CATALOG HEADER\n#\n") + nlines = nlines + 1 + + call sfree (sp) + + return (nlines) +end + + +# AT_WCATRECS -- Write out the catalog records without modification, except +# for the builtin trim parameters. + +int procedure at_wcatrecs (fd, res) + +int fd #I the output file descriptor +pointer res #I the results descriptor + +pointer sp, record +int sz_rec, nrec, recptr, nchars +int cq_rstati(), cq_gnrecord() + +begin + # Allocate space for the record. For now SZ_LINE is the default. + if (cq_rstati(res, CQRECSIZE) > 0) + sz_rec = max (SZ_LINE, cq_rstati (res, CQRECSIZE)) + else + sz_rec = SZ_LINE + nrec = cq_rstati (res, CQRNRECS) + + # Allocate working space. + call smark (sp) + call salloc (record, sz_rec, TY_CHAR) + + # For the moment assume that the simple and blocked text file records + # are newline delimited, and that the simple text file fields are + # whitespace delimited. + + # Write the records. + switch (cq_rstati (res, CQRTYPE)) { + + case CQ_STEXT: + recptr = 0 + while (recptr < nrec) { + nchars = cq_gnrecord (res, Memc[record], sz_rec, recptr) + if (nchars == EOF) + break + call fprintf (fd, "%s") + call pargstr (Memc[record]) + } + + case CQ_BTEXT: + recptr = 0 + while (recptr < nrec) { + nchars = cq_gnrecord (res, Memc[record], sz_rec, recptr) + if (nchars == EOF) + break + call fprintf (fd, "%s") + call pargstr (Memc[record]) + } + + default: + ; + } + + call sfree (sp) + + return (recptr) +end diff --git a/noao/astcat/src/agetcat/atwedit.x b/noao/astcat/src/agetcat/atwedit.x new file mode 100644 index 00000000..7678eb35 --- /dev/null +++ b/noao/astcat/src/agetcat/atwedit.x @@ -0,0 +1,83 @@ +include <imhdr.h> +include <pkg/cq.h> + +# Add a valid WCS to the image header if it does not already have one using +# the image WCS status specified by the wcs keyword in the image survey +# database. If the wcs keyord is "fits", the image is assumed to have a +# valid FITS WCS and no new wcs is computed, if it is "dss" the image is assumed +# to have a valid DSS image header which will be transformed to a valid FITS +# WCS if a FITS WCS is not already present, if it is "none" the image is +# assumed to have no valid WCS and the code will attempt to insert one using +# information in the image results structure. An error status is returned +# only if there is no valid wcs code. + +procedure at_wedit (im, res, at, wcstype, update, verbose) + +pointer im #I the input image descriptor +pointer res #I the image query results descriptor +pointer at #I the astrometry package descriptor +int wcstype #I the default wcs type +bool update #I actually update the header ? +bool verbose #I verbose mode ? + +int cq_istati(), at_mkdss(), at_dbwcs(), at_parwcs() + +begin + # Update WCS from database + if (res != NULL) { + + switch (cq_istati (res, CQWCS)) { + + # Image surveys database indicates image already has a FITS WCS. + case CQ_WFITS: + ; + + # Image surveys database indicates image has a DSS WCS. + case CQ_WDSS: + if (at_mkdss (im, update, verbose) == ERR) { + #if (update || verbose) + if (verbose) + call printf ( + " Error converting DSS wcs to FITS wcs\n") + } + + # Image surveys database indicates image has no WCS. If the proper + # information is not in the image survey then default to awcspars. + default: + if (at_dbwcs (im, res, update, verbose) == ERR) { + #if (update || verbose) + if (verbose) + call printf ( + " Error creating FITS wcs using image survey db\n") + } + } + + } else { + + switch (wcstype) { + + # User parameter indicates image already has a FITS WCS. + case CQ_WFITS: + ; + + # User parameter indicates image has a DSS WCS. + case CQ_WDSS: + if (at_mkdss (im, update, verbose) == ERR) { + #if (update || verbose) + if (verbose) + call printf ( + " Error converting DSS wcs to FITS wcs\n") + } + + default: + if (at == NULL) + ; + else if (at_parwcs (im, at, update, verbose) == ERR) { + #if (update || verbose) + if (verbose) + call printf ( + " Error creating FITS wcs using default parameters\n") + } + } + } +end diff --git a/noao/astcat/src/agetcat/mkpkg b/noao/astcat/src/agetcat/mkpkg new file mode 100644 index 00000000..c8b025b5 --- /dev/null +++ b/noao/astcat/src/agetcat/mkpkg @@ -0,0 +1,31 @@ +# AGETCAT task subdirectory + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + t_aclist.x + t_aslist.x + t_agetcat.x "../../lib/astrom.h" + t_afiltcat.x "../../lib/astrom.h" + t_agetim.x "../../lib/astrom.h" <pkg/cq.h> + t_ahedit.x "../../lib/astrom.h" <pkg/cq.h> + t_aimfind.x "../../lib/astrom.h" <pkg/cq.h> + atrcquery.x <math.h> <pkg/cq.h> <pkg/skywcs.h> \ + "../../lib/astrom.h" "../../lib/acatalog.h" + atrcrd.x <fset.h> <imhdr.h> <mwset.h> <pkg/skywcs.h> \ + "../../lib/astrom.h" "../../lib/acatalog.h" + attquery.x <math.h> <pkg/cq.h> <pkg/skywcs.h> "../../lib/astrom.h" + atwcat.x <pkg/cq.h> + atfcat.x <imhdr.h> <ctotok.h> <evvexpr.h> "../../lib/astrom.h" \ + <ctype.h> "../../lib/acatalog.h" <pkg/cq.h> \ + <pkg/skywcs.h> + athedit.x "../../lib/astrom.h" "../../lib/aimpars.h" <pkg/cq.h> + atwedit.x <imhdr.h> <pkg/cq.h> + atcatinit.x + atincat.x + atoutcat.x + atrcsym.x "../../lib/astrom.h" + ; diff --git a/noao/astcat/src/agetcat/t_aclist.x b/noao/astcat/src/agetcat/t_aclist.x new file mode 100644 index 00000000..4c13e1ed --- /dev/null +++ b/noao/astcat/src/agetcat/t_aclist.x @@ -0,0 +1,112 @@ +# T_ACLIST -- List the supported catalogs. + +procedure t_aclist() + +pointer sp, str1, str2, line, cq +int i, j, catlist, nquery, nheader, nfields +bool verbose +pointer cq_map() +int at_catlist(), fntlenb(), fntrfnb(), cq_setcat(), cq_fgeti(), cq_scan() +bool clgetb() +errchk cq_fgeti() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (str1, SZ_FNAME, TY_CHAR) + call salloc (str2, SZ_FNAME, TY_CHAR) + call salloc (line, SZ_LINE, TY_CHAR) + + # Get the parameters. + call clgstr ("catalogs", Memc[str1], SZ_FNAME) + call clgstr ("catdb", Memc[str2], SZ_FNAME) + verbose = clgetb ("verbose") + + # Get the catalog list. + catlist = at_catlist (Memc[str1], Memc[str2]) + if (fntlenb (catlist) <= 0) { + if (verbose) + call printf ("The catalog list is empty\n") + call fntclsb (catlist) + call sfree (sp) + return + } + + # Open the catalog database. + cq = cq_map (Memc[str2], READ_ONLY) + if (verbose) { + call printf ("\nScanning catalog database %s\n") + call pargstr (Memc[str2]) + } + + # Loop over the catalogs. + if (verbose) + call printf ("Listing the supported catalogs\n") + do i = 1, fntlenb (catlist) { + + # Get the catalog name and set the current catalog. + if (fntrfnb (catlist, i, Memc[str1], SZ_FNAME) == EOF) + break + if (cq_setcat (cq, Memc[str1]) <= 0) { + next + } else { + call printf ("%s\n") + call pargstr (Memc[str1]) + } + + # Do a detailed listing. + if (verbose) { + iferr (nquery = cq_fgeti (cq, "nquery")) + nquery = 0 + call printf ("nquery %d\n") + call pargi (nquery) + if (nquery > 0) { + do j = 1, nquery { + if (cq_scan (cq) == EOF) + break + call gargstr (Memc[line], SZ_LINE) + call printf ("%s\n") + call pargstr (Memc[line]) + } + } + iferr (nheader = cq_fgeti (cq, "nheader")) + nheader = 0 + call printf ("nheader %d\n") + call pargi (nheader) + if (nheader > 0) { + do j = 1, nheader { + if (cq_scan (cq) == EOF) + break + call gargstr (Memc[line], SZ_LINE) + call printf ("%s\n") + call pargstr (Memc[line]) + } + } + iferr (nfields = cq_fgeti (cq, "nfields")) + nfields = 0 + call printf ("nfields %d\n") + call pargi (nfields) + if (nfields > 0) { + do j = 1, nfields { + if (cq_scan (cq) == EOF) + break + call gargstr (Memc[line], SZ_LINE) + call printf ("%s\n") + call pargstr (Memc[line]) + } + } + if (nquery > 0 || nheader > 0 || nfields > 0) + call printf ("\n") + } + } + + + # Close the catalog database. + call cq_unmap (cq) + + # Close the catalog list. + call fntclsb (catlist) + + # Free working memory. + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_afiltcat.x b/noao/astcat/src/agetcat/t_afiltcat.x new file mode 100644 index 00000000..72dfe8b1 --- /dev/null +++ b/noao/astcat/src/agetcat/t_afiltcat.x @@ -0,0 +1,211 @@ +include "../../lib/astrom.h" + +define SZ_HDRTEXT 5 * SZ_LINE + +# T_AFILTCAT -- Filter existing astrometry catalogs. + +procedure t_afiltcat() + +pointer sp, input, output, catdb, catname, infname, outfname, tmpfname, hdrtext +pointer at, cq, res +int icatlist, ocatlist, catno, infd, outfd, nlines +bool standard, filter, update, verbose +pointer cq_map(), cq_fquery() +int fntopnb(), fntlenb(), fntgfnb(), cq_setcat(), open(), at_gcathdr() +int at_pcathdr() +bool streq(), clgetb() +errchk open() + +begin + # Allocate some working space. + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (catdb, SZ_FNAME, TY_CHAR) + call salloc (catname, SZ_FNAME, TY_CHAR) + call salloc (infname, SZ_FNAME, TY_CHAR) + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (tmpfname, SZ_FNAME, TY_CHAR) + call salloc (hdrtext, SZ_HDRTEXT, TY_CHAR) + + # Get the important query parameters. + call clgstr ("input", Memc[input], SZ_FNAME) + call clgstr ("output", Memc[output], SZ_FNAME) + call clgstr ("catdb", Memc[catdb], SZ_FNAME) + call clgstr ("catalogs", Memc[catname], SZ_FNAME) + + standard = clgetb ("standard") + filter = clgetb ("filter") + update = clgetb ("update") + verbose = clgetb ("verbose") + + # Open the input catalog list. + icatlist = fntopnb (Memc[input], NO) + ocatlist = fntopnb (Memc[output], NO) + + # Check that the input and output catalogs are the same size. + if (fntlenb (icatlist) != fntlenb (ocatlist)) { + if (verbose) { + call printf ( + "Input and output file lists lengths are different\n") + call flush (STDOUT) + } + call fntclsb (icatlist) + call fntclsb (ocatlist) + call sfree (sp) + return + } + + # Map the database. + cq = cq_map (Memc[catdb], READ_ONLY) + if (cq == NULL) { + if (verbose) { + call printf ("Cannot open catalog configuration file %s\n") + call pargstr (Memc[catdb]) + call flush (STDOUT) + } + call fntclsb (icatlist) + call fntclsb (ocatlist) + call sfree (sp) + return + } else { + if (verbose) { + call printf ("\nOpening catalog configuration file %s ...\n") + call pargstr (Memc[catdb]) + call flush (STDOUT) + } + } + + # Locate the dummy record, usually called "stext". + catno = cq_setcat (cq, Memc[catname]) + if (catno <= 0) { + if (verbose) { + call printf ("Cannot locate dummy catalog %s\n") + call pargstr (Memc[catname]) + call flush (STDOUT) + } + call cq_unmap (cq) + call fntclsb (icatlist) + call fntclsb (ocatlist) + call sfree (sp) + return + } else { + if (verbose) { + call printf ("Selecting dummy catalog %s\n") + call pargstr (Memc[catname]) + call flush (STDOUT) + } + } + + # Initilize the astrometry data structure. + call at_afinit (at) + + # Initialize the algorithm parameters. + call at_fapars (at) + + # Store the input and output templates. + call at_sets (at, CATALOGS, Memc[catname]) + call at_sets (at, INPUT, Memc[input]) + call at_sets (at, OUTPUT, Memc[output]) + call at_sets (at, CATDB, Memc[catdb]) + call at_sets (at, CATNAME, Memc[catname]) + + # Loop over the input and output files. + while (fntgfnb (icatlist, Memc[infname], SZ_FNAME) != EOF && + fntgfnb (ocatlist, Memc[outfname], SZ_FNAME) != EOF) { + + # Store the input and output catalog names. + call at_sets (at, INFNAME, Memc[infname]) + call at_sets (at, OUTFNAME, Memc[outfname]) + + # Create a temporary name and open the output file. + if (streq (Memc[infname], Memc[outfname])) + call mktemp ("tmp", Memc[tmpfname], SZ_FNAME) + else + call strcpy (Memc[outfname], Memc[tmpfname], SZ_FNAME) + iferr { + outfd = open (Memc[tmpfname], NEW_FILE, TEXT_FILE) + } then { + if (verbose) { + call printf (" Cannot open output file %s\n") + call pargstr (Memc[outfname]) + call flush (STDOUT) + next + } + } + + # Read the input catalog header. + infd = open (Memc[infname], READ_ONLY, TEXT_FILE) + nlines = at_gcathdr (infd, Memc[hdrtext], SZ_HDRTEXT) + call close (infd) + if (nlines <= 0) + nlines = at_pcathdr ("acatpars", Memc[hdrtext], SZ_HDRTEXT) + + # Read in the catalog and make it look like the results + # of a query. + if (nlines > 0) { + res = cq_fquery (cq, Memc[infname], Memc[hdrtext]) + if (res != NULL) { + if (filter) { + if (verbose) { + call printf ( + " Filtering catalog %s to catalog %s\n") + call pargstr (Memc[infname]) + call pargstr (Memc[outfname]) + } + call at_wfilrecs (outfd, at, res, standard) + } else { + if (verbose) { + call printf ( + " Copying catalog %s to catalog %s\n") + call pargstr (Memc[infname]) + call pargstr (Memc[outfname]) + } + call at_wnofilrecs (outfd, res, standard) + } + } else { + if (verbose) { + call printf (" Cannot read catalog %s\n") + call pargstr (Memc[infname]) + call flush (STDOUT) + } + } + } else { + if (verbose) { + call printf (" Cannot decode catalog %s\n") + call pargstr (Memc[infname]) + call flush (STDOUT) + } + res = NULL + } + + # Close the results structure. + if (res != NULL) + call cq_rclose (res) + + # Close the output file. + call close (outfd) + + # Replace the existing file with the temporary one. + if (streq (Memc[infname], Memc[outfname])) { + call delete (Memc[infname]) + call rename (Memc[tmpfname], Memc[infname]) + } + } + + # Free the database. + call cq_unmap (cq) + + # Update the algorithm parameters. + if (update) + call at_fppars (at) + + # Free the astrometry data structure. + call at_affree (at) + + # Free the input catalog list. + call fntclsb (icatlist) + call fntclsb (ocatlist) + + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_agetcat.x b/noao/astcat/src/agetcat/t_agetcat.x new file mode 100644 index 00000000..5f888eef --- /dev/null +++ b/noao/astcat/src/agetcat/t_agetcat.x @@ -0,0 +1,251 @@ +include "../../lib/astrom.h" + +define SZ_HDRTEXT (5 * SZ_LINE) + +procedure t_agetcat() + +pointer sp, output, hdrtext, str1, str2, at, cq, res, cres +int i, j, nfields, catlist, outlist, infd, outfd, nlines +bool standard, filter, update, verbose +pointer cq_map(), cq_query, cq_fquery(), at_tquery() +int at_rclist(), at_ocatlist(), at_catlist(), fntlenb(), cq_setcat() +int fntrfnb(), open(), at_rcquery(), access(), at_gcathdr() +bool clgetb() +errchk open() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (str1, SZ_FNAME, TY_CHAR) + call salloc (str2, SZ_FNAME, TY_CHAR) + call salloc (hdrtext, SZ_HDRTEXT, TY_CHAR) + + # Initalize the data structures + call at_aginit (at) + + # Get the iportant query parameters. + call clgstr ("regions", Memc[str1], SZ_FNAME) + call clgstr ("output", Memc[output], SZ_FNAME) + + # Get the mode parameters. + standard = clgetb ("standard") + filter = clgetb ("filter") + update = clgetb ("update") + verbose = clgetb ("verbose") + + # Allocate the astrometry structure and read in the algorithm + # parameters. This must be done before the field centers are + # decoded. + call at_gapars (at) + + # Print the field center parameters. + #call at_rcshow (at) + # Print the filtering parameters. + #call at_fsshow (at) + # Print the wcs parameters. + #call at_wcshow (at) + # Print the image parameters. + #call at_imshow (at) + + # Get the field center list. + nfields = at_rclist (at, Memc[str1]) + if (nfields <= 0) { + if (verbose) + call printf ("The field center list is empty\n") + call at_agfree (at) + call sfree (sp) + return + } + + # Print the field center symbol table. + #call at_stshow (at) + + # Get the catalog list. + call clgstr ("catalogs", Memc[str1], SZ_FNAME) + call clgstr ("catdb", Memc[str2], SZ_FNAME) + catlist = at_catlist (Memc[str1], Memc[str2]) + if (fntlenb (catlist) <= 0) { + if (verbose) + call printf ("The catalog list is empty\n") + call fntclsb (catlist) + call at_agfree (at) + call sfree (sp) + return + } + call at_sets (at, CATALOGS, Memc[str1]) + call at_sets (at, CATDB, Memc[str2]) + + + # Print the i/o parameters. + #call at_ioshow (at) + + # Create the output catalog file list. + outlist = at_ocatlist (at, catlist, Memc[output], "default", "cat", NO) + if (fntlenb (outlist) <= 0) { + if (verbose) + call printf ("The output file list is empty\n") + call fntclsb (outlist) + call fntclsb (catlist) + call at_agfree (at) + call sfree (sp) + return + } + call at_sets (at, OUTPUT, Memc[output]) + + # Open the catalog database. + cq = cq_map (Memc[str2], READ_ONLY) + if (verbose) { + call printf ("\nOpening catalog database %s\n") + call pargstr (Memc[str2]) + } + + # Loop over the catalog list. + do i = 1, fntlenb (catlist) { + + # Get the catalog name and save it. + if (fntrfnb (catlist, i, Memc[str2], SZ_FNAME) == EOF) + break + if (access (Memc[str2], READ_ONLY, TEXT_FILE) == YES) { + if (cq_setcat (cq, "filename@noao") <= 0) { + if (verbose) { + call printf ("Skipping catalog %s\n") + call pargstr (Memc[str2]) + call flush (STDOUT) + } + next + } else { + call at_sets (at, CATNAME, Memc[str2]) + if (verbose) { + call printf ("Selecting catalog %s\n") + call pargstr (Memc[str2]) + call flush (STDOUT) + } + } + } else if (cq_setcat (cq, Memc[str2]) <= 0) { + if (verbose) { + call printf ("Skipping catalog %s\n") + call pargstr (Memc[str2]) + call flush (STDOUT) + } + next + } else { + call at_sets (at, CATNAME, Memc[str2]) + if (verbose) { + call printf ("Selecting catalog %s\n") + call pargstr (Memc[str2]) + call flush (STDOUT) + } + } + + # Loop over the field centers. + do j = 1, nfields { + + # Get the output file name. + if (fntrfnb (outlist, (i - 1) * nfields + j, Memc[str1], + SZ_FNAME) == EOF) + break + call at_sets (at, OUTFNAME, Memc[str1]) + + # Open the output file. + iferr { + outfd = open (Memc[str1], NEW_FILE, TEXT_FILE) + } then { + if (verbose) { + call printf (" Unable to open output file %s\n") + call pargstr (Memc[str1]) + } + break + } + + if (access (Memc[str2], READ_ONLY, TEXT_FILE) == YES) { + + # Read the catalog header. + infd = open (Memc[str2], READ_ONLY, TEXT_FILE) + nlines = at_gcathdr (infd, Memc[hdrtext], SZ_HDRTEXT) + call close (infd) + if (nlines <= 0) { + if (verbose) + call printf (" Unable to read catalog header\n") + break + } + + # Copy the catalog file into the query structure. + cres = cq_fquery (cq, Memc[str2], Memc[hdrtext]) + if (cres == NULL) { + call printf (" Catalog query failed\n") + break + } + + # Extract the requested data. + res = at_tquery (at, cq, cres, Memc[hdrtext], nlines, j) + if (res == NULL) { + if (verbose) + call printf (" Catalog query failed\n") + break + } + + } else { + + # Format the query. + if (at_rcquery (at, cq, j) == ERR) { + if (verbose) + call printf (" Unable to format network query\n") + break + } + + # Query the catalog. + res = cq_query (cq) + if (res == NULL) { + if (verbose) + call printf (" Network query failed\n") + break + } + } + + # Write the output file. + if (filter) { + if (verbose) { + call printf (" Filtering region %d to file %s\n") + call pargi (j) + call pargstr (Memc[str1]) + call flush (STDOUT) + } + call at_wfilrecs (outfd, at, res, standard) + } else { + if (verbose) { + call printf (" Copying region %d to file %s\n") + call pargi (j) + call pargstr (Memc[str1]) + call flush (STDOUT) + } + call at_wnofilrecs (outfd, res, standard) + } + + # Close the output file. + call close (outfd) + + # Close the query structure. + call cq_rclose (res) + + } + + } + + # Close the catalog database. + call cq_unmap (cq) + + # Update the algorithm parameters. + if (update) + call at_gppars (at) + + # Close the catalog and output file lists. + call fntclsb (outlist) + call fntclsb (catlist) + + # Free the astrometry structure. + call at_agfree (at) + + # Free the working memory. + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_agetim.x b/noao/astcat/src/agetcat/t_agetim.x new file mode 100644 index 00000000..687c1229 --- /dev/null +++ b/noao/astcat/src/agetcat/t_agetim.x @@ -0,0 +1,247 @@ +include "../../lib/astrom.h" +include <pkg/cq.h> + +define SZ_IMEXTN 10 + +procedure t_agetim() + +pointer sp, output, extn, str1, str2 +pointer cq, at, im, res +int i, j, index, nfields, svlist, imlist, addext +char period +bool wcsedit, hdredit, update, verbose +pointer cq_map(), immap(), cq_imquery() +int at_rclist(), at_svlist(), at_osvlist(), fntlenb(), cq_setcat() +int at_rcquery(), fntrfnb(), strldx(), imaccess(), imtlen(), imtrgetim() +int open() +bool clgetb(), streq() +data period /'.'/ +errchk open(), immap(), imaccess(), cq_fgstr() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (str1, SZ_FNAME, TY_CHAR) + call salloc (str2, SZ_FNAME, TY_CHAR) + call salloc (extn, SZ_IMEXTN, TY_CHAR) + + # Initalize the data structures + call at_aiginit (at) + + # Get the iportant query parameters. + call clgstr ("regions", Memc[str1], SZ_FNAME) + call clgstr ("images", Memc[output], SZ_FNAME) + + # Get the editing parameters. + wcsedit = clgetb ("wcsedit") + hdredit = clgetb ("hdredit") + update = clgetb ("update") + verbose = clgetb ("verbose") + + # Allocate the astrometry structure and read in the algorithm + # parameters. This must be done before the field centers are + # decoded. + call at_giapars (at) + + # Print the field center parameters. + #call at_rcshow (at) + # Print the default wcs parameters. + #call at_wcshow (at) + # Print the default image data parameters. + #call at_imshow (at) + + # Get the field center list. + nfields = at_rclist (at, Memc[str1]) + if (nfields <= 0) { + if (verbose) + call printf ("The field center list is empty\n") + call at_aigfree (at) + call sfree (sp) + return + } + + # Print the field center symbol table. + #call at_stshow (at) + + # Get the surverys list. + call clgstr ("imsurveys", Memc[str1], SZ_FNAME) + call clgstr ("imdb", Memc[str2], SZ_FNAME) + svlist = at_svlist (Memc[str1], Memc[str2]) + if (fntlenb (svlist) <= 0) { + if (verbose) + call printf ("The image surveys list is empty\n") + call at_aigfree (at) + call fntclsb (svlist) + call sfree (sp) + return + } + call at_sets (at, SURVEYS, Memc[str1]) + call at_sets (at, IMDB, Memc[str2]) + + # Print the i/o parameters. + #call at_ioshow (at) + + # Create the output image list. + imlist = at_osvlist (at, svlist, Memc[output], "default", "", NO) + if (imtlen (imlist) <= 0) { + if (verbose) + call printf ("The output images list is empty\n") + call at_aigfree (at) + call imtclose (imlist) + call fntclsb (svlist) + call sfree (sp) + return + } + call at_sets (at, IMAGES, Memc[output]) + + # Open the catalog database. + cq = cq_map (Memc[str2], READ_ONLY) + if (verbose) { + call printf ("\nOpening surveys database %s\n") + call pargstr (Memc[str2]) + } + + # Loop over the catalog list. + do i = 1, fntlenb (svlist) { + + # Get the catalog name and save it. + if (fntrfnb (svlist, i, Memc[str1], SZ_FNAME) == EOF) + break + if (cq_setcat (cq, Memc[str1]) <= 0) { + if (verbose) { + call printf ("Skipping survey %s\n") + call pargstr (Memc[str1]) + call flush (STDOUT) + } + next + } else { + call at_sets (at, SVNAME, Memc[str1]) + if (verbose) { + call printf ("Selecting survey %s\n") + call pargstr (Memc[str1]) + call flush (STDOUT) + } + } + + # Loop over the field centers. + do j = 1, nfields { + + # Get the output file name. + if (imtrgetim (imlist, (i - 1) * nfields + j, Memc[str1], + SZ_FNAME) == EOF) + break + + # If the file is a fits file tack on the user extension. This + # is not the correct way to do this but for the moment it will + # work. Not sure there is a totally clean way to do this since + # we are not going through imio. + + ifnoerr { + call cq_fgstr (cq, "type", Memc[extn], SZ_IMEXTN) + } then { + addext = YES + index = strldx (period, Memc[str1]) + if (index > 0) { + if (streq (Memc[extn], Memc[str1+index])) + addext = NO + else + addext = YES + } + if (addext == YES) { + call strcpy (Memc[str1], Memc[str2], SZ_FNAME) + call strcat (".", Memc[str2], SZ_FNAME) + call strcat (Memc[extn], Memc[str2], SZ_FNAME) + call strcpy (Memc[str2], Memc[str1], SZ_FNAME) + } + } else { + if (verbose) + call printf ( + " Warning the image format is undefined\n") + } + call at_sets (at, IMNAME, Memc[str1]) + + # Can the output file be opened ? + iferr { + im = open (Memc[str1], NEW_FILE, BINARY_FILE) + } then { + if (verbose) { + call printf (" Unable to write output image %s\n") + call pargstr (Memc[str1]) + } + break + } else { + call close (im) + call delete (Memc[str1]) + } + + # Format the query. + if (at_rcquery (at, cq, j) == ERR) { + if (verbose) + call printf (" Unable to format network query\n") + break + } + + # Query the image surveys. + if (verbose) { + call printf ("Getting image %s ...\n") + call pargstr (Memc[str1]) + call flush (STDOUT) + } + res = cq_imquery (cq, Memc[str1]) + if (res == NULL) { + if (verbose) + call printf (" Network query failed\n") + next + } + + # Open the output file. + iferr { + if (imaccess (Memc[str1], READ_WRITE) == YES) { + if (wcsedit || hdredit) { + im = immap (Memc[str1], READ_WRITE, 0) + if (wcsedit) + call at_wedit (im, res, NULL, CQ_WNONE, true, + verbose) + if (hdredit) + call at_hedit (im, res, NULL, true, verbose) + call imunmap (im) + } + } else + im = NULL + } then { + if (verbose) { + call printf ( + " Warning %s is not a valid image\n") + call pargstr (Memc[str1]) + } + im = NULL + } + + if (verbose) + call flush (STDOUT) + + # Close the query structure. + call cq_imclose (res) + + } + + } + + # Close the catalog database. + call cq_unmap (cq) + + # Update the algorithm parameters. + if (update) + call at_gippars (at) + + # Close the file and image lists. + call imtclose (imlist) + call fntclsb (svlist) + + # Free the astrometry structure. + call at_aigfree (at) + + # Free the working memory. + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_ahedit.x b/noao/astcat/src/agetcat/t_ahedit.x new file mode 100644 index 00000000..c75ce3a7 --- /dev/null +++ b/noao/astcat/src/agetcat/t_ahedit.x @@ -0,0 +1,175 @@ +include "../../lib/astrom.h" +include <pkg/cq.h> + +procedure t_ahedit() + +pointer sp, images, str1, str2 +pointer at, cq, res, im +int j, imlist, catno, wcstype +bool hupdate, wcsedit, hdredit, update, verbose +bool clgetb() +pointer cq_map(), cq_fimquery(), immap() +int imtopen(), imtlen(), cq_setcat(), imtrgetim(), imaccess(), strdic() +errchk immap(), imaccess() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (images, SZ_FNAME, TY_CHAR) + call salloc (str1, SZ_FNAME, TY_CHAR) + call salloc (str2, SZ_FNAME, TY_CHAR) + + # Get the iportant query parameters. + call clgstr ("images", Memc[images], SZ_FNAME) + + # Get the editing parameters. + hupdate = clgetb ("hupdate") + wcsedit = clgetb ("wcsedit") + call clgstr ("wcs", Memc[str1], SZ_FNAME) + wcstype = strdic (Memc[str1], Memc[str1], SZ_FNAME, CQ_WTYPESTR) + if (wcstype <= 0) + wcstype = CQ_WNONE + hdredit = clgetb ("hdredit") + update = clgetb ("update") + if (hupdate) + verbose = clgetb ("verbose") + else + verbose = true + + # Open the image list. + imlist = imtopen (Memc[images]) + if (imtlen (imlist) <= 0) { + if (verbose) + call printf ("The input image list is empty\n") + call imtclose (imlist) + call sfree (sp) + return + } + + # Initalize the data structures + call at_ahinit (at) + + # Allocate the astrometry structure and read in the algorithm + # parameters. + call at_hapars (at) + + # Print the default wcs parameters. + #call at_wcshow (at) + # Print the default image data parameters. + #call at_imshow (at) + + # Set the i/o parameters. + call clgstr ("imsurveys", Memc[str1], SZ_FNAME) + call clgstr ("imdb", Memc[str2], SZ_FNAME) + call at_sets (at, IMAGES, Memc[images]) + call at_sets (at, SURVEYS, Memc[str1]) + call at_sets (at, IMDB, Memc[str2]) + + # Print the i/o parameters. + #call at_ioshow (at) + + # Open the catalog database. + cq = cq_map (Memc[str2], READ_ONLY) + if (cq == NULL) { + if (verbose) { + call printf ("\nCannot opening surveys database %s\n") + call pargstr (Memc[str2]) + } + call at_ahfree (at) + call imtclose (imlist) + call sfree (sp) + return + } else if (verbose) { + call printf ("\nOpening surveys database %s\n") + call pargstr (Memc[str2]) + call flush (STDOUT) + } + + # Get the catalog name and save it. + catno = cq_setcat (cq, Memc[str1]) + if (Memc[str1] == EOS) { + catno = ERR + } else if (catno == 0) { + if (verbose) { + call printf ("Cannot locate survey %s\n") + call pargstr (Memc[str1]) + } + } else { + if (verbose) { + call printf ("Selecting survey %s\n") + call pargstr (Memc[str1]) + } + call at_sets (at, CATNAME, Memc[str1]) + } + + # Loop over the field centers. + do j = 1, imtlen (imlist) { + + # Get the output image name. + if (imtrgetim (imlist, j, Memc[str1], SZ_FNAME) == EOF) + break + call at_sets (at, IMNAME, Memc[str1]) + + # Query the image survey to get the header info even though the + # image already exists. + if (verbose) { + call printf ("Getting image %s ...\n") + call pargstr (Memc[str1]) + call flush (STDOUT) + } + if (catno <= 0) + res = NULL + else + res = cq_fimquery (cq, Memc[str1]) + + # Open the output file. + iferr { + if (imaccess (Memc[str1], READ_WRITE) == YES) { + im = immap (Memc[str1], READ_WRITE, 0) + if (wcsedit) { + if (res != NULL) + call at_wedit (im, res, NULL, wcstype, hupdate, + verbose) + else + call at_wedit (im, NULL, at, wcstype, hupdate, + verbose) + } + if (hdredit) { + if (res != NULL) + call at_hedit (im, res, NULL, hupdate, verbose) + else + call at_hedit (im, NULL, at, hupdate, verbose) + } + call imunmap (im) + } else + im = NULL + } then { + if (verbose) { + call printf (" Warning %s is not a valid image name\n") + call pargstr (Memc[str1]) + } + im = NULL + } + + # Close the query structure. + if (res != NULL) + call cq_imclose (res) + + } + + # Close the catalog database. + call cq_unmap (cq) + + # Update the algorithm parameters. + if (update) + call at_hppars (at) + + # Close the image lists. + call imtclose (imlist) + + # Free the astrometry structure. + call at_ahfree (at) + + # Free the working memory. + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_aimfind.x b/noao/astcat/src/agetcat/t_aimfind.x new file mode 100644 index 00000000..58601bef --- /dev/null +++ b/noao/astcat/src/agetcat/t_aimfind.x @@ -0,0 +1,318 @@ +include <pkg/cq.h> +include "../../lib/astrom.h" + +define SZ_HDRTEXT (5 * SZ_LINE) + +procedure t_aimfind() + +pointer sp, images, output, imfile, catalog, catdb, hdrtext, str1 +pointer at, cq, cres, res, sym, im +int i, j, nfields, nout, nlines, catlist, outlist, imfd, infd, outfd +bool standard, filter, append, update, verbose +pointer cq_map(), cq_query(), cq_fquery(), at_tquery(), at_rcsym(), immap() +int at_rclist(), at_catlist(), at_ocatlist(), open(), access() +int fntlenb(), fntrfnb(), cq_setcat(), at_rcquery() +int at_gcathdr(), cq_rstati() +bool clgetb() +errchk open() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (images, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (imfile, SZ_FNAME, TY_CHAR) + call salloc (catalog, SZ_FNAME, TY_CHAR) + call salloc (catdb, SZ_FNAME, TY_CHAR) + call salloc (hdrtext, SZ_HDRTEXT, TY_CHAR) + call salloc (str1, SZ_FNAME, TY_CHAR) + + # Initalize the data structures + call at_aginit (at) + + # Get the iportant query parameters. + call clgstr ("images", Memc[images], SZ_FNAME) + call clgstr ("output", Memc[output], SZ_FNAME) + call clgstr ("imfile", Memc[imfile], SZ_FNAME) + + # Get the mode parameters. + standard = clgetb ("standard") + filter = clgetb ("filter") + append = clgetb ("append") + update = clgetb ("update") + verbose = clgetb ("verbose") + + # Allocate the astrometry structure and read in the algorithm + # parameters. If filtering is turned off then the filtering + # parameters are set to their default values. Probably need to + # make high level wrapper routines for the parameter defaults + # routines at some point. + if (! filter) + call at_dfspset(at) + else + call at_iapars (at) + + # Set the new field values and new field descriptions. At some + # point these may be input from parameters. At present they are + # hardwired. + call at_nflist (at, 2, "xp,yp", "d,d", "pixels,pixels", + "%10.3f,%10.3f", append) + + # Print the field center parameters. + #call at_rcshow (at) + # Print the filtering parameters. + #call at_fsshow (at) + # Print the wcs parameters. + #call at_wcshow (at) + # Print the image parameters. + #call at_imshow (at) + + # Create the region list from the image list. If an image does not + # have a valid fits wcs it will not be included in the valid + # region list. + nfields = at_rclist (at, Memc[images]) + if (nfields <= 0) { + if (verbose) + call printf ("The image list is empty\n") + call at_agfree (at) + call sfree (sp) + return + } + +# # Print the field center symbol table. +# #call at_stshow (at) + + # Get the catalog. The catalog may be a catalog server or an + # astrometry file. + call clgstr ("catalogs", Memc[catalog], SZ_FNAME) + call clgstr ("catdb", Memc[catdb], SZ_FNAME) + catlist = at_catlist (Memc[catalog], Memc[catdb]) + if (fntlenb (catlist) != 1) { + if (verbose) { + if (fntlenb (catlist) <= 0) + call printf ("The catalog is undefined\n") + else + call printf ("More than one catalog is specified\n") + + } + call fntclsb (catlist) + call at_agfree (at) + call sfree (sp) + return + } + call at_sets (at, CATALOGS, Memc[catalog]) + call at_sets (at, CATDB, Memc[catdb]) + + # Open the output image list file. If the output image file name + # is imdefined then no image list file is written. + if (Memc[imfile] == EOS) { + imfd = NULL + } else { + iferr (imfd = open (Memc[imfile], NEW_FILE, TEXT_FILE)) + imfd = NULL + } + +# # Print the i/o parameters. +# #call at_ioshow (at) + + # Create the output astrometry file list. If the output astrometry + # file list is empty no astrometry file is written. + outlist = at_ocatlist (at, catlist, Memc[output], "default", "coo", NO) + call at_sets (at, OUTPUT, Memc[output]) + + # Open the catalog database. + cq = cq_map (Memc[catdb], READ_ONLY) + if (verbose) { + call printf ("\nOpening catalog database %s\n") + call pargstr (Memc[catdb]) + } + + # Loop over the catalog list. + nout = 0 + do i = 1, fntlenb (catlist) { + + # Get the catalog name and save it. + if (fntrfnb (catlist, i, Memc[catalog], SZ_FNAME) == EOF) + break + + # Set the catalog. + if (access (Memc[catalog], READ_ONLY, TEXT_FILE) == YES) { + if (cq_setcat (cq, "filename@noao") <= 0) { + if (verbose) { + call printf ("Skipping catalog %s\n") + call pargstr (Memc[catalog]) + call flush (STDOUT) + } + next + } else { + call at_sets (at, CATNAME, Memc[catalog]) + if (verbose) { + call printf ("Selecting catalog %s\n") + call pargstr (Memc[catalog]) + call flush (STDOUT) + } + } + } else if (cq_setcat (cq, Memc[catalog]) <= 0) { + if (verbose) { + call printf ("Skipping catalog %s\n") + call pargstr (Memc[catalog]) + call flush (STDOUT) + } + next + } else { + call at_sets (at, CATNAME, Memc[catalog]) + if (verbose) { + call printf ("Selecting catalog %s\n") + call pargstr (Memc[catalog]) + call flush (STDOUT) + } + } + + # Loop over the field centers. + do j = 1, nfields { + + # Get the output file name. + if (fntrfnb (outlist, (i - 1) * nfields + j, Memc[str1], + SZ_FNAME) == EOF) + call at_sets (at, OUTFNAME, "") + #break + else + call at_sets (at, OUTFNAME, Memc[str1]) + + # Query the catalog. + if (access (Memc[catalog], READ_ONLY, TEXT_FILE) == YES) { + + # Read the catalog header. + infd = open (Memc[catalog], READ_ONLY, TEXT_FILE) + nlines = at_gcathdr (infd, Memc[hdrtext], SZ_HDRTEXT) + call close (infd) + if (nlines <= 0) { + if (verbose) + call printf (" Unable to read catalog header\n") + break + } + + # Copy the standard star catalog into the query structure. + cres = cq_fquery (cq, Memc[catalog], Memc[hdrtext]) + if (cres == NULL) { + if (verbose) + call printf (" Catalog query failed\n") + break + } + + # Extract the requested data. + res = at_tquery (at, cq, cres, Memc[hdrtext], nlines, j) + if (res == NULL) { + if (verbose) + call printf (" Catalog query failed\n") + break + } + call cq_rclose (cres) + + + } else { + + # Format the network query. + if (at_rcquery (at, cq, j) == ERR) { + if (verbose) + call printf (" Unable to format network query\n") + break + } + + # Query the catalog. + res = cq_query (cq) + if (res == NULL) { + if (verbose) + call printf (" Network query failed\n") + break + } + } + + # Get the region symbol. + sym = at_rcsym (at, j) + + # If at least one object was detected in the image then + # write out the catalog for that image, and add the image + # name to the image file. + if (cq_rstati (res, CQRNRECS) > 0) { + + # Print the number of objects found. + if (verbose) { + call printf ( + " Image %s contains %d catalog objects\n") + call pargstr (AT_RCSTNAME(sym)) + call pargi (cq_rstati(res, CQRNRECS)) + call flush (STDOUT) + } + + # Write the query results to the astrometry file. + outfd = NULL + if (fntlenb (outlist) > 0) { + + # Open the output file. + outfd = open (Memc[str1], NEW_FILE, TEXT_FILE) + if (verbose) { + call printf (" Writing catalog file %s\n") + call pargstr (Memc[str1]) + } + + im = immap (AT_RCSTNAME(sym), READ_ONLY, 0) + call at_wifilrecs (outfd, im, at, res, standard) + call imunmap (im) + + # Close the output file. + call close (outfd) + } + + # Write the image name to the image file. + if (imfd != NULL) { + if (outfd != NULL) { + call fprintf (imfd, "%s %s\n") + call pargstr (AT_RCSTNAME(sym)) + call pargstr (Memc[str1]) + } else { + call fprintf (imfd, "%s\n") + call pargstr (AT_RCSTNAME(sym)) + } + } + + # Count the number of non-empty files + nout = nout + 1 + + } else if (verbose) { + call printf (" Image %s contains no catalog objects\n") + call pargstr (AT_RCSTNAME(sym)) + call flush (STDOUT) + } + + # Close the query structure. + call cq_rclose (res) + + } + + } + + # Close the catalog database. + call cq_unmap (cq) + + # Update the algorithm parameters. + if (update) + call at_ippars (at) + + # Close the catalog and output file lists. + call fntclsb (outlist) + call fntclsb (catlist) + + # Close the image list file. Delete it if it is empty. + if (imfd != NULL) { + call close (imfd) + if (nout <= 0) + call delete (Memc[imfile]) + } + + # Free the astrometry structure. + call at_agfree (at) + + # Free the working memory. + call sfree (sp) +end diff --git a/noao/astcat/src/agetcat/t_aslist.x b/noao/astcat/src/agetcat/t_aslist.x new file mode 100644 index 00000000..19f7a55d --- /dev/null +++ b/noao/astcat/src/agetcat/t_aslist.x @@ -0,0 +1,102 @@ +# T_ASLIST -- List the support image surveys. + +procedure t_aslist() + +pointer sp, str1, str2, line, cq +int i, j, svlist, nwcs, nkeys +bool verbose +pointer cq_map() +int at_svlist(), fntlenb(), fntrfnb(), cq_setcat(), cq_fgeti(), cq_scan() +bool clgetb() +errchk cq_fgstr(), cq_fgeti() + +begin + # Allocate some working memory. + call smark (sp) + call salloc (str1, SZ_FNAME, TY_CHAR) + call salloc (str2, SZ_FNAME, TY_CHAR) + call salloc (line, SZ_LINE, TY_CHAR) + + # Get the parameters. + call clgstr ("imsurveys", Memc[str1], SZ_FNAME) + call clgstr ("imdb", Memc[str2], SZ_FNAME) + verbose = clgetb ("verbose") + + # Get the catalog list. + svlist = at_svlist (Memc[str1], Memc[str2]) + if (fntlenb (svlist) <= 0) { + if (verbose) + call printf ("The image surveys list is empty\n") + call fntclsb (svlist) + call sfree (sp) + return + } + + # Open the catalog database. + cq = cq_map (Memc[str2], READ_ONLY) + if (verbose) { + call printf ("\nScanning image surveys database %s\n") + call pargstr (Memc[str2]) + } + + # Loop over the catalogs. + if (verbose) + call printf ("Listing the supported image surveys\n") + do i = 1, fntlenb (svlist) { + + # Get the catalog name and set the current catalog. + if (fntrfnb (svlist, i, Memc[str1], SZ_FNAME) == EOF) + break + if (cq_setcat (cq, Memc[str1]) <= 0) { + next + } else { + call printf ("%s\n") + call pargstr (Memc[str1]) + } + + # Do a detailed listing. + if (verbose) { + iferr (call cq_fgstr (cq, "wcs", Memc[line], SZ_LINE)) + call strcpy ("none", Memc[line], SZ_LINE) + call printf ("wcs %s\n") + call pargstr (Memc[line]) + iferr (nwcs = cq_fgeti (cq, "nwcs")) + nwcs = 0 + call printf ("nwcs %d\n") + call pargi (nwcs) + if (nwcs > 0) { + do j = 1, nwcs { + if (cq_scan (cq) == EOF) + break + call gargstr (Memc[line], SZ_LINE) + call printf ("%s\n") + call pargstr (Memc[line]) + } + } + iferr (nkeys = cq_fgeti (cq, "nkeys")) + nkeys = 0 + call printf ("nkeys %d\n") + call pargi (nkeys) + if (nkeys > 0) { + do j = 1, nkeys { + if (cq_scan (cq) == EOF) + break + call gargstr (Memc[line], SZ_LINE) + call printf ("%s\n") + call pargstr (Memc[line]) + } + } + if (nwcs > 0 || nkeys > 0) + call printf ("\n") + } + } + + # Close the image surveys database. + call cq_unmap (cq) + + # Close the image surveys list. + call fntclsb (svlist) + + # Free working memory. + call sfree (sp) +end |