/*** File libwcs/tabread.c *** August 29, 2014 *** By Jessica Mink, jmink@cfa.harvard.edu *** Harvard-Smithsonian Center for Astrophysics *** Copyright (C) 1996-2014 *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Correspondence concerning WCSTools should be addressed as follows: Internet email: jmink@cfa.harvard.edu Postal address: Jessica Mink Smithsonian Astrophysical Observatory 60 Garden St. Cambridge, MA 02138 USA */ /* int tabread() Read tab table stars in specified region * int tabrnum() Read tab table stars with specified numbers * int tabbin () Bin tab table stars in specified region into an image * int tabxyread() Read x, y, and magnitude from tab table star list * int tabrkey() Read single keyword from specified tab table stars * struct StarCat tabcatopen() Open tab table catalog, return number of entries * struct TabTable *tabopen() Open tab table, returning number of entries * char *tabline() Get tab table entry for one line * double tabgetra() Return double right ascension in degrees * double tabgetdec() Return double declination in degrees * double tabgetr8() Return 8-byte floating point number from tab table line * int tabgeti4() Return 4-byte integer from tab table line * int tabgetk() Return character entry from tab table line for column * int tabgetc() Return n'th character entry from tab table line * int tabhgetr8() Return 8-byte floating point keyword value from header * int tabhgeti4() Return 4-byte integer keyword value from header * int tabhgetc() Return character keyword value from header * int tabparse() Make a table of column headings * int tabcol() Find entry in a table of column headings (case-dependent) * int tabccol() Find entry in a table of column headings (case-independent) * int tabsize() Return length of file in bytes * int istab() Return 1 if first line of file contains a tab, else 0 */ #include #include #include #include #include #include #include "wcs.h" #include "fitsfile.h" #include "wcscat.h" #define ABS(a) ((a) < 0 ? (-(a)) : (a)) static int tabhgetr8(); static int tabhgeti4(); static int tabhgetc(); static int tabcont(); static int tabccont(); static int tabsize(); static int nndec = 0; static int verbose = 0; static char *taberr; static struct Tokens startok; char *gettaberr () { return (taberr); } int gettabndec() { return (nndec); } static char *kwo = NULL; /* Keyword returned by tabread(), tabrnum() */ void settabkey (keyword0) char *keyword0; { kwo = keyword0; return; } /* TABREAD -- Read tab table stars in specified region */ int tabread (tabcatname,distsort,cra,cdec,dra,ddec,drad,dradi, sysout,eqout,epout,mag1,mag2,sortmag,nstarmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tpeak,tkey,nlog) char *tabcatname; /* Name of reference star catalog file */ int distsort; /* 1 to sort stars by distance from center */ double cra; /* Search center J2000 right ascension in degrees */ double cdec; /* Search center J2000 declination in degrees */ double dra; /* Search half width in right ascension in degrees */ double ddec; /* Search half-width in declination in degrees */ double drad; /* Limiting separation in degrees (ignore if 0) */ double dradi; /* Inner edge of annulus in degrees (ignore if 0) */ int sysout; /* Search coordinate system */ double eqout; /* Search coordinate equinox */ double epout; /* Proper motion epoch (0.0 for no proper motion) */ double mag1,mag2; /* Limiting magnitudes (none if equal) */ int sortmag; /* Magnitude by which to sort (1 to nmag) */ int nstarmax; /* Maximum number of stars to be returned */ struct StarCat **starcat; /* Star catalog data structure */ double *tnum; /* Array of UJ numbers (returned) */ double *tra; /* Array of right ascensions (returned) */ double *tdec; /* Array of declinations (returned) */ double *tpra; /* Array of right ascension proper motions (returned) */ double *tpdec; /* Array of declination proper motions (returned) */ double **tmag; /* 2-D Array of magnitudes (returned) */ int *tpeak; /* Array of peak counts (returned) */ char **tkey; /* Array of values of additional keyword */ int nlog; { double ra1,ra2; /* Limiting right ascensions of region in degrees */ double dec1,dec2; /* Limiting declinations of region in degrees */ double rra1,rra2; /* Catalog coordinate limiting right ascension (deg) */ double rdec1,rdec2; /* Catalog coordinate limiting declination (deg) */ double dist = 0.0; /* Distance from search center in degrees */ double faintmag=0.0; /* Faintest magnitude */ double maxdist=0.0; /* Largest distance */ int faintstar=0; /* Faintest star */ int farstar=0; /* Most distant star */ double *tdist; /* Array of distances to stars */ int sysref; /* Catalog coordinate system */ double eqref; /* Catalog equinox */ double epref; /* Catalog epoch */ double secmarg = 0.0; /* Arcsec/century margin for proper motion */ double magt; double rdist, ddist; int pass; char cstr[32]; struct Star *star; struct StarCat *sc; /* Star catalog data structure */ int wrap; int jstar; int magsort; int nstar; char *objname; int lname; int imag; double ra,dec, rapm, decpm; double mag, parallax, rv; double num; int peak, i; int istar, nstars, lstar; sc = *starcat; if (nlog > 0) verbose = 1; else verbose = 0; SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose); /* mag1 is always the smallest magnitude limit */ if (mag2 < mag1) { mag = mag2; mag2 = mag1; mag1 = mag; } /* Logging interval */ nstar = 0; tdist = (double *) calloc (nstarmax, sizeof (double)); lstar = sizeof (struct Star); star = (struct Star *) calloc (1, lstar); if (sc == NULL) sc = tabcatopen (tabcatname, NULL, 0); *starcat = sc; if (sc == NULL || sc->nstars <= 0) { if (taberr != NULL) fprintf (stderr,"%s\n", taberr); fprintf (stderr,"TABREAD: Cannot read catalog %s\n", tabcatname); free (star); sc = NULL; return (0); } nstars = sc->nstars; jstar = 0; if (sortmag > 0 && sortmag <= sc->nmag) magsort = sortmag - 1; else magsort = 0; /* Set catalog coordinate system */ if (sc->equinox != 0.0) eqref = sc->equinox; else eqref = eqout; if (sc->epoch != 0.0) epref = sc->epoch; else epref = epout; if (sc->coorsys) sysref = sc->coorsys; else sysref = sysout; wcscstr (cstr, sysout, eqout, epout); rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; if (sc->mprop > 0) secmarg = 60.0; RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,secmarg, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); /* Loop through catalog */ for (istar = 1; istar <= nstars; istar++) { /* Read position of next star */ if (tabstar (istar, sc, star, verbose)) { if (verbose) fprintf (stderr,"TABREAD: Cannot read star %d\n", istar); break; } /* Set magnitude to test */ if (sc->nmag > 0) { magt = star->xmag[magsort]; if (sortmag < 1) { imag = 0; while (magt == 99.90 && imag < sc->nmag) magt = star->xmag[imag++]; if (magt > 100.0) magt = magt - 100.0; } } else magt = mag1; /* Check magnitude limits */ pass = 1; if (mag1 != mag2 && (magt < mag1 || magt > mag2)) pass = 0; /* Check rough position limits */ ra = star->ra; dec = star->dec; if ((!wrap && (ra < rra1 || ra > rra2)) || (wrap && (ra < rra1 && ra > rra2)) || dec < rdec1 || dec > rdec2) pass = 0; /* Convert coordinate system for this star and test it*/ if (pass) { sysref = star->coorsys; eqref = star->equinox; epref = star->epoch; /* Extract selected fields */ num = star->num; rapm = star->rapm; decpm = star->decpm; parallax = star->parallax; rv = star->radvel; /* Convert from catalog to search coordinate system */ if (sc->entpx || sc->entrv) wcsconv (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm, ¶llax, &rv); else if (sc->mprop == 1) wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); else wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); if (sc->sptype) peak = (1000 * (int) star->isp[0]) + (int)star->isp[1]; else peak = star->peak; /* Compute distance from search center */ if (drad > 0 || distsort) dist = wcsdist (cra,cdec,ra,dec); else dist = 0.0; /* Check radial distance to search center */ if (drad > 0) { if (dist > drad) pass = 0; if (dradi > 0.0 && dist < dradi) pass = 0; } /* Check distance along RA and Dec axes */ else { ddist = wcsdist (cra,cdec,cra,dec); if (ddist > ddec) pass = 0; rdist = wcsdist (cra,dec,ra,dec); if (rdist > dra) pass = 0; } } /* Add this star's information to the list if still OK */ if (pass) { /* Save star position and magnitude in table */ if (nstar < nstarmax) { tnum[nstar] = num; tra[nstar] = ra; tdec[nstar] = dec; if (sc->mprop == 1) { tpra[nstar] = rapm; tpdec[nstar] = decpm; } for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][nstar] = star->xmag[imag]; } if (tpeak) tpeak[nstar] = peak; tdist[nstar] = dist; lname = strlen (star->objname); if (lname > 0) { objname = (char *)calloc (lname+1, 1); strcpy (objname, star->objname); if (tkey[nstar]) free(tkey[nstar]); tkey[nstar] = objname; } if (dist > maxdist) { maxdist = dist; farstar = nstar; } if (sc->nmag > 0 && magt > faintmag) { faintmag = magt; faintstar = nstar; } } /* If radial search & too many stars, replace furthest star */ else if (distsort) { if (dist < maxdist) { tnum[farstar] = num; tra[farstar] = ra; tdec[farstar] = dec; if (sc->mprop == 1) { tpra[farstar] = rapm; tpdec[farstar] = decpm; } for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][farstar] = star->xmag[imag]; } tpeak[farstar] = peak; tdist[farstar] = dist; lname = strlen (star->objname); if (lname > 0) { objname = (char *)calloc (lname+1, 1); strcpy (objname, star->objname); if (tkey[farstar]) free(tkey[farstar]); tkey[farstar] = objname; } /* Find new farthest star */ maxdist = 0.0; for (i = 0; i < nstarmax; i++) { if (tdist[i] > maxdist) { maxdist = tdist[i]; farstar = i; } } } } /* Otherwise if too many stars, replace faintest star */ else if (sc->nmag > 0 && magt < faintmag) { tnum[faintstar] = num; tra[faintstar] = ra; tdec[faintstar] = dec; if (sc->mprop == 1) { tpra[faintstar] = rapm; tpdec[faintstar] = decpm; } for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][faintstar] = star->xmag[imag]; } tpeak[faintstar] = peak; tdist[faintstar] = dist; lname = strlen (star->objname); if (lname > 0) { objname = (char *)calloc (lname+1, 1); strcpy (objname, star->objname); if (tkey[faintstar]) free(tkey[faintstar]); tkey[faintstar] = objname; } faintmag = 0.0; /* Find new faintest star */ for (i = 0; i < nstarmax; i++) { magt = tmag[magsort][i]; imag = 0; while (magt == 99.90 && imag < sc->nmag) magt = tmag[imag++][i]; if (magt > 100.0) magt = magt - 100.0; if (magt > faintmag) { faintmag = magt; faintstar = i; } } } nstar++; jstar++; if (nlog == 1) fprintf (stderr,"TABREAD: %11.6f: %9.5f %9.5f %s %5.2f %d \n", num,ra,dec,cstr,magt,peak); /* End of accepted star processing */ } /* Log operation */ if (nlog > 0 && istar%nlog == 0) fprintf (stderr,"TABREAD: %5d / %5d / %5d sources catalog %s\r", jstar,istar,nstars,tabcatname); /* End of star loop */ } /* Summarize search */ if (nlog > 0) { fprintf (stderr,"TABREAD: Catalog %s : %d / %d / %d found\n",tabcatname, jstar,istar,nstars); if (nstar > nstarmax) fprintf (stderr,"TABREAD: %d stars found; only %d returned\n", nstar,nstarmax); } free ((char *) tdist); free ((char *) star); return (nstar); } /* TABRNUM -- Read tab table stars with specified numbers */ int tabrnum (tabcatname, nnum, sysout, eqout, epout, starcat, match, tnum,tra,tdec,tpra,tpdec,tmag,tpeak,tkey,nlog) char *tabcatname; /* Name of reference star catalog file */ int nnum; /* Number of stars to look for */ int sysout; /* Search coordinate system */ double eqout; /* Search coordinate equinox */ double epout; /* Proper motion epoch (0.0 for no proper motion) */ struct StarCat **starcat; /* Star catalog data structure */ int match; /* 1 to match star number exactly, else sequence num.*/ double *tnum; /* Array of star numbers to look for */ double *tra; /* Array of right ascensions (returned) */ double *tdec; /* Array of declinations (returned) */ double *tpra; /* Array of right ascension proper motions (returned) */ double *tpdec; /* Array of declination proper motions (returned) */ double **tmag; /* 2-D array of magnitudes (returned) */ int *tpeak; /* Array of peak counts (returned) */ char **tkey; /* Array of additional keyword values */ int nlog; { int jnum; int nstar; double ra,dec, rapm, decpm; double mag, parallax, rv; double num; int peak; int istar, istar0, nstars; int imag; char *line; char numstr[32]; /* Catalog number */ int sysref; /* Catalog coordinate system */ double eqref; /* Catalog equinox */ double epref; /* Catalog epoch */ char cstr[32]; char str[32]; char *objname; char *datestring; int lname, lstar; int ireg, inum, nnfld, i; char rastr[32], decstr[32]; struct TabTable *startab; struct StarCat *sc; struct Star *star; line = NULL; nnfld = 0; nstar = 0; nndec = 0; /* Allocate catalog entry buffer */ lstar = sizeof (struct Star); star = (struct Star *) calloc (1, lstar); /* Open star catalog */ sc = *starcat; if (sc == NULL) sc = tabcatopen (tabcatname, NULL, 0); *starcat = sc; if (sc == NULL || sc->nstars <= 0) { if (taberr != NULL) fprintf (stderr,"%s\n", taberr); fprintf (stderr,"TABRNUM: Cannot read catalog %s\n", tabcatname); free (star); return (0); } startab = sc->startab; nstars = sc->nstars; /* Set catalog coordinate system */ if (sc->equinox != 0.0) eqref = sc->equinox; else eqref = eqout; if (sc->epoch != 0.0) epref = sc->epoch; else epref = epout; if (sc->coorsys) sysref = sc->coorsys; else sysref = sysout; wcscstr (cstr, sysout, eqout, epout); /* Write header if printing star entries as found */ if (nlog < 0) { char *revmessage; revmessage = getrevmsg(); printf ("catalog %s\n", tabcatname); printf ("radecsys %s\n", cstr); printf ("equinox %.3f\n", eqout); printf ("epoch %.3f\n", epout); printf ("program scat %s\n", revmessage); if (sc->nnfld > 0) nnfld = sc->nnfld; else nnfld = CatNumLen (TABCAT, tnum[nnum-1], sc->nndec); printf ("id ra dec "); for (i = 1; i < sc->nmag+1; i++) { if (i == sc->nmag && sc->entepoch) printf (" epoch "); else printf (" %s", sc->keymag[i-1]); } if (kwo != NULL) printf (" %s\n", kwo); else printf ("\n"); printf ("----- ------------ ------------"); for (i = 1; i < sc->nmag+1; i++) { if (i == sc->nmag && sc->entepoch) printf (" --------"); else printf (" -----"); } if (kwo != NULL) printf (" ---------\n"); else printf ("\n"); } star->num = 0.0; istar0 = 0; num = 0.0; /* Loop through star list */ line = startab->tabdata; for (jnum = 0; jnum < nnum; jnum++) { /* Read forward from the last star if possible */ inum = (int) (tnum[jnum] + 0.5); if ((double)inum != tnum[jnum] || inum < istar0 || istar0 == 0) { istar0 = 1; } /* Loop through catalog to star */ for (istar = istar0; istar <= nstars; istar++) { if (!match && istar == inum) break; if (num < tnum[jnum]) { if ((line = gettabline (startab, istar)) == NULL) { if (nlog) fprintf (stderr,"TABRNUM: Cannot read star %d\n", istar); break; } } /* Check ID number first */ (void) setoken (&startok, line, "tab"); if (!strcmp (sc->isfil,"gsc-server")) { if (tabgetc (&startok, sc->entid, str, 24)) num = 0.0; else { num = atof (str+3) * 0.00001; ireg = (int) num; inum = (int) (((num - (double)ireg) * 100000.0) + 0.5); num = (double) ireg + 0.0001 * (double) inum; } } else num = tabgetr8 (&startok,sc->entid); if (num == 0.0) num = (double) istar; if (num > tnum[jnum]) { break; } if (num == tnum[jnum]) break; } /* If star has been found in table, read rest of entry */ if ((match && num == tnum[jnum]) || (!match && inum == istar)) { istar0 = istar; sc->istar = startab->iline; if (tabstar (istar, sc, star, nlog)) fprintf (stderr,"TABRNUM: Cannot read star %d\n", istar); /* If star entry has been read successfully */ else { /* Set coordinate system for this star */ sysref = star->coorsys; eqref = star->equinox; /* Extract selected fields */ num = star->num; ra = star->ra; dec = star->dec; rapm = star->rapm; decpm = star->decpm; parallax = star->parallax; rv = star->radvel; if (sc->entrv > 0) star->xmag[sc->nmag-1] = rv; if (sc->entpx || sc->entrv) wcsconv (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm, ¶llax, &rv); else if (sc->mprop == 1) wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); else wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); if (sc->nmag > 0) mag = star->xmag[0]; else mag = 99.99; if (sc->sptype) peak = (1000 * (int) star->isp[0]) + (int)star->isp[1]; else peak = star->peak; if (nlog < 0) { CatNum (TABCAT, -nnfld, sc->nndec, num, numstr); ra2str (rastr, 31, ra, 3); dec2str (decstr, 31, dec, 2); printf ("%s %s %s", numstr,rastr,decstr); for (imag = 0; imag < sc->nmag; imag++) if (imag == sc->nmag-1 && sc->entepoch) { datestring = DateString (star->xmag[imag], 1); printf ("%s", datestring); free (datestring); } else printf (" %.2f", star->xmag[imag]); if (kwo != NULL) printf (" %s\n", star->objname); else printf ("\n"); continue; } /* Save star position and magnitude in table */ tnum[jnum] = num; tra[jnum] = ra; tdec[jnum] = dec; if (sc->mprop == 1) { tpra[jnum] = rapm; tpdec[jnum] = decpm; } for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][jnum] = star->xmag[imag]; } tpeak[jnum] = peak; lname = strlen (star->objname); if (lname > 0) { objname = (char *)calloc (lname+1, 1); strcpy (objname, star->objname); if (tkey[jnum]) free(tkey[jnum]); tkey[jnum] = objname; } nstar++; if (nlog == 1) fprintf (stderr,"TABRNUM: %11.6f: %9.5f %9.5f %s %5.2f %d \n", num,ra,dec,cstr,mag,peak); /* End of accepted star processing */ } } else { nstar++; istar0 = istar; if (nlog < 0) { CatNum (TABCAT, -nnfld, sc->nndec, tnum[jnum], numstr); ra = 0.0; ra2str (rastr, 31, ra, 3); dec = 0.0; dec2str (decstr, 31, dec, 2); printf ("%s %s %s", numstr,rastr,decstr); for (imag = 0; imag < sc->nmag; imag++) printf (" 99.0"); printf ("\n"); } } /* Log operation */ if (nlog > 0 && jnum%nlog == 0) fprintf (stderr,"TABRNUM: %5d / %5d / %5d sources catalog %s\r", nstar,jnum,nstars,tabcatname); /* End of star loop */ } /* Summarize search */ if (nlog > 0) fprintf (stderr,"TABRNUM: Catalog %s : %d / %d found\n", tabcatname,nstar,nstars); free ((char *) star); return (nstar); } /* TABBIN -- Bin tab table stars in specified region into an image */ int tabbin (tabcatname, wcs, header, image, mag1, mag2, sortmag, magscale, nlog) char *tabcatname; /* Name of reference star catalog file */ struct WorldCoor *wcs; /* World coordinate system for image */ char *header; /* FITS header for output image */ char *image; /* Output FITS image */ double mag1,mag2; /* Limiting magnitudes (none if equal) */ int sortmag; /* Magnitude by which to sort (1 to nmag) */ double magscale; /* Scaling factor for magnitude to pixel flux * (number of catalog objects per bin if 0) */ int nlog; { double cra; /* Search center J2000 right ascension in degrees */ double cdec; /* Search center J2000 declination in degrees */ double dra; /* Search half width in right ascension in degrees */ double ddec; /* Search half-width in declination in degrees */ int sysout; /* Search coordinate system */ double eqout; /* Search coordinate equinox */ double epout; /* Proper motion epoch (0.0 for no proper motion) */ double ra1,ra2; /* Limiting right ascensions of region in degrees */ double dec1,dec2; /* Limiting declinations of region in degrees */ double rra1,rra2; /* Catalog coordinate limiting right ascension (deg) */ double rdec1,rdec2; /* Catalog coordinate limiting declination (deg) */ int sysref; /* Catalog coordinate system */ double eqref; /* Catalog equinox */ double epref; /* Catalog epoch */ double secmarg = 0.0; /* Arcsec/century margin for proper motion */ double magt; double rdist, ddist; int ix, iy; int pass; char cstr[32]; struct Star *star; struct StarCat *sc; /* Star catalog data structure */ int wrap; int jstar; int magsort; int nstar; int imag; double ra,dec, rapm, decpm; double mag, parallax, rv; double num; int peak; int istar, nstars, lstar; double xpix, ypix, flux; int offscl; int bitpix, w, h; /* Image bits/pixel and pixel width and height */ double logt = log(10.0); if (nlog > 0) verbose = 1; else verbose = 0; /* Set image parameters */ bitpix = 0; (void)hgeti4 (header, "BITPIX", &bitpix); w = 0; (void)hgeti4 (header, "NAXIS1", &w); h = 0; (void)hgeti4 (header, "NAXIS2", &h); /* Set catalog search limits from image WCS information */ sysout = wcs->syswcs; eqout = wcs->equinox; epout = wcs->epoch; wcscstr (cstr, sysout, eqout, epout); wcssize (wcs, &cra, &cdec, &dra, &ddec); SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose); /* mag1 is always the smallest magnitude limit */ if (mag2 < mag1) { mag = mag2; mag2 = mag1; mag1 = mag; } /* Logging interval */ nstar = 0; lstar = sizeof (struct Star); star = (struct Star *) calloc (1, lstar); sc = tabcatopen (tabcatname, NULL, 0); if (sc == NULL || sc->nstars <= 0) { if (taberr != NULL) fprintf (stderr,"%s\n", taberr); fprintf (stderr,"TABBIN: Cannot read catalog %s\n", tabcatname); free (star); sc = NULL; return (0); } nstars = sc->nstars; jstar = 0; if (sortmag > 0 && sortmag <= sc->nmag) magsort = sortmag - 1; else magsort = 0; /* Set catalog coordinate system */ if (sc->equinox != 0.0) eqref = sc->equinox; else eqref = eqout; if (sc->epoch != 0.0) epref = sc->epoch; else epref = epout; if (sc->coorsys) sysref = sc->coorsys; else sysref = sysout; wcscstr (cstr, sysout, eqout, epout); rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; if (sc->mprop > 0) secmarg = 60.0; RefLim(cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,secmarg, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); /* If RA range includes zero, split it in two */ if (rra1 > rra2) wrap = 1; else wrap = 0; /* Loop through catalog */ for (istar = 1; istar <= nstars; istar++) { /* Read position of next star */ if (tabstar (istar, sc, star, verbose)) { if (verbose) fprintf (stderr,"TABBIN: Cannot read star %d\n", istar); break; } /* Set magnitude to test */ if (sc->nmag > 0) { magt = star->xmag[magsort]; imag = 0; while (magt == 99.90 && imag < sc->nmag) magt = star->xmag[imag++]; if (magt > 100.0) magt = magt - 100.0; } else magt = mag1; /* Check magnitude limits */ pass = 1; if (mag1 != mag2 && (magt < mag1 || magt > mag2)) pass = 0; /* Check rough position limits */ ra = star->ra; dec = star->dec; if ((!wrap && (ra < rra1 || ra > rra2)) || (wrap && (ra < rra1 && ra > rra2)) || dec < rdec1 || dec > rdec2) pass = 0; /* Convert coordinate system for this star and test it*/ if (pass) { sysref = star->coorsys; eqref = star->equinox; epref = star->epoch; /* Extract selected fields */ num = star->num; rapm = star->rapm; decpm = star->decpm; parallax = star->parallax; rv = star->radvel; /* Convert from catalog to search coordinate system */ if (sc->entpx || sc->entrv) wcsconv (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm, ¶llax, &rv); else if (sc->mprop == 1) wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); else wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); if (sc->sptype) peak = (1000 * (int) star->isp[0]) + (int)star->isp[1]; else peak = star->peak; /* Check distance along RA and Dec axes */ ddist = wcsdist (cra,cdec,cra,dec); if (ddist > ddec) pass = 0; rdist = wcsdist (cra,dec,ra,dec); if (rdist > dra) pass = 0; } /* Save star in FITS image */ if (pass) { wcs2pix (wcs, ra, dec, &xpix, &ypix, &offscl); if (!offscl) { if (magscale > 0.0) flux = magscale * exp (logt * (-magt / 2.5)); else flux = 1.0; ix = (int) (xpix + 0.5); iy = (int) (ypix + 0.5); addpix1 (image, bitpix, w,h, 0.0,1.0, xpix,ypix, flux); nstar++; jstar++; } else { ix = 0; iy = 0; } if (nlog == 1) { fprintf (stderr,"TABBIN: %11.6f: %9.5f %9.5f %s", num,ra,dec,cstr); if (magscale > 0.0) fprintf (stderr, " %5.2f", mag); if (!offscl) flux = getpix1 (image, bitpix, w, h, 0.0, 1.0, ix, iy); else flux = 0.0; fprintf (stderr," (%d,%d): %f\n", ix, iy, flux); } /* End of accepted star processing */ } /* Log operation */ if (nlog > 0 && istar%nlog == 0) fprintf (stderr,"TABBIN: %5d / %5d / %5d sources catalog %s\r", jstar,istar,nstars,tabcatname); /* End of star loop */ } /* Summarize search */ if (nlog > 0) { fprintf (stderr,"TABBIN: Catalog %s : %d / %d / %d found\n",tabcatname, jstar,istar,nstars); } return (nstar); } /* TABXYREAD -- Read X, Y, and magnitude of tab table stars */ int tabxyread (tabcatname, xa, ya, ba, pa, nlog) char *tabcatname; /* Name of reference star catalog file */ double **xa; /* Array of x coordinates (returned) */ double **ya; /* Array of y coordinates (returned) */ double **ba; /* Array of magnitudes (returned) */ int **pa; /* Array of fluxes (returned) */ int nlog; { double xi, yi, magi, flux; char *line; int istar, nstars; struct TabTable *startab; int entx, enty, entmag; /* Open tab table file */ nndec = 0; startab = tabopen (tabcatname, 0); if (startab == NULL || startab->nlines <= 0) { fprintf (stderr,"TABXYREAD: Cannot read catalog %s\n", tabcatname); return (0); } /* Find columns for X and Y */ entx = tabccol (startab, "x"); enty = tabccol (startab, "y"); /* Find column for magnitude */ if (!(entmag = tabccol (startab, "mag"))) { if (!(entmag = tabccol (startab, "magv"))) { if (!(entmag = tabccol (startab, "magj"))) entmag = tabccol (startab, "magr"); } } /* Allocate vectors for x, y, magnitude, and flux */ nstars = startab->nlines; *xa = (double *) realloc(*xa, nstars*sizeof(double)); if (*xa == NULL) { fprintf (stderr,"TABXYREAD: Cannot allocate memory for x\n"); return (0); } *ya = (double *) realloc(*ya, nstars*sizeof(double)); if (*ya == NULL) { fprintf (stderr,"TABXYREAD: Cannot allocate memory for y\n"); return (0); } *ba = (double *) realloc(*ba, nstars*sizeof(double)); if (*ba == NULL) { fprintf (stderr,"TABXYREAD: Cannot allocate memory for mag\n"); return (0); } *pa = (int *) realloc(*pa, nstars*sizeof(int)); if (*pa == NULL) { fprintf (stderr,"TABXYREAD: Cannot allocate memory for flux\n"); return (0); } /* Loop through catalog */ for (istar = 0; istar < nstars; istar++) { /* Read line for next star */ if ((line = gettabline (startab, istar+1)) == NULL) { fprintf (stderr,"TABXYREAD: Cannot read star %d\n", istar); break; } /* Extract x, y, and magnitude */ (void) setoken (&startok, line, "tab"); xi = tabgetr8 (&startok, entx); yi = tabgetr8 (&startok, enty); magi = tabgetr8 (&startok, entmag); (*xa)[istar] = xi; (*ya)[istar] = yi; (*ba)[istar] = magi; flux = 1000000000.0 * pow (10.0, (-magi / 2.5)); (*pa)[istar] = (int) flux; if (nlog == 1) fprintf (stderr,"DAOREAD: %6d/%6d: %9.5f %9.5f %6.2f %15.4f\n", istar,nstars,xi,yi,magi,flux); /* Log operation */ if (nlog > 1 && istar%nlog == 0) fprintf (stderr,"TABXYREAD: %5d / %5d sources catalog %s\r", istar,nstars,tabcatname); /* End of star loop */ } /* Summarize search */ if (nlog > 0) fprintf (stderr,"TABXYREAD: Catalog %s : %d / %d found\n",tabcatname, istar,nstars); /* Free table */ tabclose (startab); if (istar < nstars-1) return (istar + 1); else return (nstars); } #define TABMAX 64 /* TABRKEY -- Read single keyword from tab table stars with specified numbers */ int tabrkey (tabcatname, starcat, nnum, tnum, keyword, tval) char *tabcatname; /* Name of reference star catalog file */ struct StarCat **starcat; /* Star catalog data structure */ int nnum; /* Number of stars to look for */ double *tnum; /* Array of star numbers to look for */ char *keyword; /* Keyword for which to return values */ char **tval; /* Returned values for specified keyword */ { int jnum, lval; int nstar; int istar, nstars; double num; char *line; char *tvalue; char value[TABMAX]; struct TabTable *startab; struct StarCat *sc; /* Star catalog data structure */ nstar = 0; /* Open star catalog */ sc = *starcat; if (sc == NULL) sc = tabcatopen (tabcatname, NULL, 0); *starcat = sc; if (sc == NULL || sc->nstars <= 0) { if (taberr != NULL) fprintf (stderr,"%s\n", taberr); fprintf (stderr,"TABRKEY: Cannot read catalog %s\n", tabcatname); return (0); } startab = sc->startab; if (startab == NULL || startab->nlines <= 0) { fprintf (stderr,"TABRKEY: Cannot read catalog %s\n", tabcatname); return (0); } /* Loop through star list */ nstars = startab->nlines; for (jnum = 0; jnum < nnum; jnum++) { /* Loop through catalog to star */ for (istar = 1; istar <= nstars; istar++) { if ((line = gettabline (startab, istar)) == NULL) { fprintf (stderr,"TABRKEY: Cannot read star %d\n", istar); num = 0.0; break; } /* Check ID number */ (void) setoken (&startok, line, "tab"); if ((num = tabgetr8 (&startok,sc->entid)) == 0.0) num = (double) istar; if (num == tnum[jnum]) break; } /* If star has been found in table */ if (num == tnum[jnum]) { nstar++; /* Extract selected field */ (void) tabgetk (startab, &startok, keyword, value, TABMAX); lval = strlen (value); if (lval > 0) { tvalue = (char *) calloc (1, lval+1); strcpy (tvalue, value); } else tvalue = NULL; if (tval[jnum]) free(tval[jnum]); tval[jnum] = tvalue; } } return (nstars); } static char newline = 10; static char tab = 9; /* TABCATOPEN -- Open tab table catalog, returning number of entries */ struct StarCat * tabcatopen (tabpath, tabtable, nbbuff) char *tabpath; /* Tab table catalog file pathname */ struct TabTable *tabtable; int nbbuff; /* Number of bytes in buffer; 0=read whole file */ { char cstr[32]; char keyword[16]; char *tabname; struct TabTable *startab; struct StarCat *sc; int i, lnum, ndec, istar, nbsc, icol, j, nnfld; int entpmq, entnid; char *line; double dnum; /* Open the tab table file */ if (tabtable != NULL) startab = tabtable; else if ((startab = tabopen (tabpath, nbbuff)) == NULL) return (NULL); /* Allocate catalog data structure */ nbsc = sizeof (struct StarCat); sc = (struct StarCat *) calloc (1, nbsc); sc->startab = startab; /* Save name of catalog */ tabname = strrchr (tabpath, '/'); if (tabname) tabname = tabname + 1; else tabname = tabpath; if (strlen (tabname) < 24) strcpy (sc->isfil, tabname); else { strncpy (sc->isfil, tabname, 23); sc->isfil[23] = (char) 0; } /* Find column and name of object identifier */ sc->entid = -1; sc->keyid[0] = (char) 0; if ((sc->entid = tabccol (startab, "id"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } else if ((sc->entid = tabccont (startab, "_id"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } else if ((sc->entid = tabccont (startab, "id"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } else if ((sc->entid = tabccont (startab, "num"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } else if ((sc->entid = tabccont (startab, "name"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } else if ((sc->entid = tabccont (startab, "obj"))) { i = sc->entid - 1; strncpy (sc->keyid, startab->colname[i], startab->lcol[i]); } sc->nndec = nndec; /* Figure out what the coordinate system is */ if (tabhgetc (startab, "radecsys", cstr)) { if (!strcmp (cstr, "galactic")) sc->coorsys = WCS_GALACTIC; else if (!strcmp (cstr, "ecliptic")) sc->coorsys = WCS_ECLIPTIC; else if (!strcmp (cstr, "fk4")) sc->coorsys = WCS_B1950; else sc->coorsys = WCS_J2000; } else sc->coorsys = WCS_J2000; /* Find column and name of object right ascension */ sc->entra = -1; sc->keyra[0] = (char) 0; if (sc->coorsys == WCS_GALACTIC) { if ((sc->entra = tabccol (startab, "long_gal"))) strcpy (sc->keyra, "long_gal"); else if ((sc->entra = tabccol (startab, "long"))) strcpy (sc->keyra, "long_gal"); else if ((sc->entra = tabccont (startab, "long"))) { i = sc->entra - 1; strncpy (sc->keyra, startab->colname[i], startab->lcol[i]); } } else if (sc->coorsys == WCS_ECLIPTIC) { if ((sc->entra = tabccol (startab, "long_ecl"))) strcpy (sc->keyra, "long_ecl"); else if ((sc->entra = tabccol (startab, "long"))) strcpy (sc->keyra, "long_ecl"); else if ((sc->entra = tabccont (startab, "long"))) { i = sc->entra - 1; strncpy (sc->keyra, startab->colname[i], startab->lcol[i]); } } else { if ((sc->entra = tabccol (startab, "ra"))) { i = sc->entra - 1; strncpy (sc->keyra, startab->colname[i], startab->lcol[i]); } else if ((sc->entra = tabccont (startab, "ra"))) { i = sc->entra - 1; strncpy (sc->keyra, startab->colname[i], startab->lcol[i]); } } /* Find column and name of object declination */ sc->entdec = -1; sc->keydec[0] = (char) 0; if (sc->coorsys == WCS_GALACTIC) { if ((sc->entdec = tabccol (startab, "lat_gal"))) strcpy (sc->keydec, "lat_gal"); else if ((sc->entdec = tabccol (startab, "lat"))) strcpy (sc->keydec, "lat_gal"); else if ((sc->entdec = tabccont (startab, "lat"))) { i = sc->entdec - 1; strncpy (sc->keydec, startab->colname[i], startab->lcol[i]); } } else if (sc->coorsys == WCS_ECLIPTIC) { if ((sc->entdec = tabccol (startab, "lat_ecl"))) strcpy (sc->keydec, "lat_ecl"); else if ((sc->entdec = tabccol (startab, "lat"))) strcpy (sc->keydec, "lat_ecl"); else if ((sc->entdec = tabccont (startab, "lat"))) { i = sc->entdec - 1; strncpy (sc->keydec, startab->colname[i], startab->lcol[i]); } } else { if ((sc->entdec = tabccol (startab, "de(deg)"))) { i = sc->entdec - 1; strcpy (sc->keydec, "dec"); } else if ((sc->entdec = tabccol (startab, "dec"))) { i = sc->entdec - 1; strncpy (sc->keydec, startab->colname[i], startab->lcol[i]); } else if ((sc->entdec = tabccont (startab, "dec"))) { i = sc->entdec - 1; strncpy (sc->keydec, startab->colname[i], startab->lcol[i]); } } /* Check columns for magnitudes and save columns and names */ sc->nmag = 0; icol = 0; for (i = 0; i < startab->ncols; i++) { icol = icol + 1; if (icol == sc->entid) continue; if (icol == sc->entra) continue; if (icol == sc->entdec) continue; for (j = 0; j < 16; j++) keyword[j] = (char)0; if (startab->lcol[i] < 16) strncpy (keyword, startab->colname[i], startab->lcol[i]); else strncpy (keyword, startab->colname[i], 15); if (strcsrch (keyword, "mag") && !strcsrch (keyword, "err") && !strcsrch (keyword, "code")) { strcpy (sc->keymag[sc->nmag], keyword); sc->entmag[sc->nmag] = icol; sc->nmag++; } else if (strcsrch (keyword, "exp")) { strcpy (sc->keymag[sc->nmag], keyword); sc->entmag[sc->nmag] = icol; sc->nmag++; } else if (strcsrch (keyword, "dist")) { strcpy (sc->keymag[sc->nmag], keyword); sc->entmag[sc->nmag] = icol; sc->nmag++; } else if ((keyword[0] == 'M' || keyword[0] == 'm') && !strcsrch (keyword, "err") && !strcsrch (keyword, "code") && !strcsrch (keyword, "flag")) { strcpy (sc->keymag[sc->nmag], keyword); sc->entmag[sc->nmag] = icol; sc->nmag++; } } /* Find column and name of object right ascension proper motion */ sc->entrpm = -1; sc->keyrpm[0] = (char) 0; if ((sc->entrpm = tabccol (startab, "ura"))) strcpy (sc->keyrpm, "ura"); else if ((sc->entrpm = tabccol (startab, "rapm"))) strcpy (sc->keyrpm, "rapm"); else if ((sc->entrpm = tabccol (startab, "pmra"))) strcpy (sc->keyrpm, "pmra"); else if ((sc->entrpm = tabccol (startab, "dra"))) strcpy (sc->keyrpm, "dra"); else if ((sc->entrpm = tabccol (startab, "ux"))) strcpy (sc->keyrpm, "ux"); /* Find column and name of object declination proper motion */ sc->entdpm = -1; sc->keydpm[0] = (char) 0; if ((sc->entdpm = tabccol (startab, "udec"))) strcpy (sc->keydpm, "udec"); else if ((sc->entdpm = tabccol (startab, "decpm"))) strcpy (sc->keydpm, "decpm"); else if ((sc->entdpm = tabccol (startab, "pmdec"))) strcpy (sc->keydpm, "pmdec"); else if ((sc->entrpm = tabccol (startab, "ddec"))) strcpy (sc->keyrpm, "ddec"); else if ((sc->entdpm = tabccol (startab, "uy"))) strcpy (sc->keydpm, "uy"); /* Find units for RA proper motion */ sc->mprop = 0; cstr[0] = 0; if (!tabhgetc (startab,"RPMUNIT", cstr)) { if (!tabhgetc (startab,"rpmunit", cstr)) { if (!tabhgetc (startab,"pmunit", cstr)) { if (sc->entdpm > 0 && sc->entrpm > 0) strcpy (cstr,"mas/yr"); } } } if (strlen (cstr) > 0) { sc->mprop = 1; if (!strcmp (cstr, "mas/yr") || !strcmp (cstr, "mas/year")) sc->rpmunit = PM_MASYR; else if (!strcmp (cstr, "ms/yr") || !strcmp (cstr, "millisec/year")) sc->rpmunit = PM_MTSYR; else if (!strcmp (cstr, "arcsec/yr") || !strcmp (cstr, "arcsec/year")) sc->rpmunit = PM_ARCSECYR; else if (!strcmp (cstr, "arcsec/hr") || !strcmp (cstr, "arcsec/hour")) sc->rpmunit = PM_ARCSECHR; else if (!strcmp (cstr, "arcsec/cen") || !strcmp (cstr, "arcsec/century")) sc->rpmunit = PM_ARCSECCEN; else if (!strcmp (cstr, "rad/yr") || !strcmp (cstr, "rad/year")) sc->rpmunit = PM_RADYR; else if (!strcmp (cstr, "sec/yr") || !strcmp (cstr, "sec/year")) sc->rpmunit = PM_TSECYR; else if (!strcmp (cstr, "tsec/yr") || !strcmp (cstr, "tsec/year")) sc->rpmunit = PM_TSECYR; else if (!strcmp (cstr, "tsec/cen") || !strcmp (cstr, "tsec/century")) sc->rpmunit = PM_TSECCEN; else sc->rpmunit = PM_DEGYR; } else if (sc->entrpm > 0) sc->rpmunit = PM_TSECCEN; /* Find units for Dec proper motion */ cstr[0] = 0; if (!tabhgetc (startab,"DPMUNIT", cstr)) { if (!tabhgetc (startab,"dpmunit", cstr)) { if (!tabhgetc (startab,"pmunit", cstr)) tabhgetc (startab,"pmunit", cstr); } } if (strlen (cstr) > 0) { sc->mprop = 1; if (!strcmp (cstr, "mas/yr") || !strcmp (cstr, "mas/year")) sc->dpmunit = PM_MASYR; else if (!strcmp (cstr, "sec/yr") || !strcmp (cstr, "sec/year")) sc->dpmunit = PM_ARCSECYR; else if (!strcmp (cstr, "tsec/yr") || !strcmp (cstr, "tsec/year")) sc->dpmunit = PM_ARCSECYR; else if (!strcmp (cstr, "arcsec/hr") || !strcmp (cstr, "arcsec/hour")) sc->dpmunit = PM_ARCSECHR; else if (!strcmp (cstr, "arcsec/cen") || !strcmp (cstr, "arcsec/century")) sc->dpmunit = PM_ARCSECCEN; else if (!strcmp (cstr, "arcsec/yr") || !strcmp (cstr, "arcsec/year")) sc->dpmunit = PM_ARCSECYR; else if (!strcmp (cstr, "rad/yr") || !strcmp (cstr, "rad/year")) sc->dpmunit = PM_RADYR; else sc->dpmunit = PM_DEGYR; } else if (sc->entdpm > 0) sc->dpmunit = PM_ARCSECCEN; /* Find column for parallax */ sc->entpx = 0; sc->entpx = tabccol (startab, "px"); /* Find column for parallax error */ sc->entpxe = 0; sc->entpxe = tabccol (startab, "pxerr"); /* Find column for radial velocity */ sc->entrv = 0; sc->keyrv[0] = (char) 0; if ((sc->entrv = tabccol (startab, "rv"))) strcpy (sc->keyrv, "rv"); else if ((sc->entrv = tabccol (startab, "cz"))) strcpy (sc->keyrv, "cz"); if (sc->entrv > 0 && sc->nmag < 10) { strcpy (sc->keymag[sc->nmag], sc->keyrv); sc->entmag[sc->nmag] = sc->entrv; sc->nmag++; } /* Find column for epoch */ sc->entepoch = 0; sc->keyepoch[0] = (char) 0; if ((sc->entepoch = tabccol (startab, "epoch"))) strcpy (sc->keyepoch, "epoch"); else if ((sc->entepoch = tabccol (startab, "ep"))) strcpy (sc->keyepoch, "ep"); if (sc->entepoch > 0 && sc->nmag < 10) { strcpy (sc->keymag[sc->nmag], sc->keyepoch); sc->entmag[sc->nmag] = sc->entepoch; sc->nmag++; sc->nepoch = 1; } /* Find column for date (in FITS format as epoch alternate) */ sc->entdate = 0; sc->entdate = tabccol (startab, "date"); /* Find column and name of object peak or plate number */ sc->entpeak = -1; sc->keypeak[0] = (char) 0; if ((sc->entpeak = tabccol (startab, "PEAK"))) strcpy (sc->keypeak, "PEAK"); else if ((sc->entpeak = tabccol (startab, "peak"))) strcpy (sc->keypeak, "peak"); else if ((sc->entpeak = tabccol (startab, "plate"))) { strcpy (sc->keypeak, "plate"); sc->plate = 1; } else if ((sc->entpeak = tabccol (startab, "field"))) { strcpy (sc->keypeak, "field"); sc->plate = 1; } else if ((sc->entpeak = tabccol (startab, "Class"))) { strcpy (sc->keypeak, "class"); sc->plate = 1; } else if ((sc->entpeak = tabcol (startab, "c"))) { strcpy (sc->keypeak, "class"); sc->plate = 1; } else if ((entpmq = tabccol (startab, "pm")) > 0 && (entnid = tabccol (startab, "ni")) > 0) { sc->entpeak = (entpmq * 100) + entnid; } /* Find column and name of object spectral type */ sc->enttype = 0; sc->keytype[0] = (char) 0; if ((sc->enttype = tabccol (startab, "SpT"))) strcpy (sc->keytype, "spt"); else if ((sc->enttype = tabccol (startab, "TYPE"))) strcpy (sc->keytype, "type"); else if ((sc->enttype = tabccont (startab, "typ"))) { i = sc->enttype - 1; strncpy (sc->keytype, startab->colname[i], startab->lcol[i]); } if (sc->enttype > 0) sc->sptype = 1; sc->entadd = -1; sc->keyadd[0] = (char) 0; if (kwo != NULL) { sc->entadd = tabccol (startab, kwo); strcpy (sc->keyadd, kwo); } /* Set catalog coordinate system */ sc->coorsys = 0; sc->equinox = 0.0; sc->epoch = 0.0; if (tabhgetc (startab,"RADECSYS", cstr)) { sc->coorsys = wcscsys (cstr); if (!tabhgetr8 (startab,"EQUINOX", &sc->equinox)) sc->equinox = wcsceq (cstr); if (!tabhgetr8 (startab,"EPOCH",&sc->epoch)) sc->epoch = sc->equinox; } if (tabhgetc (startab,"radecsys", cstr)) { sc->coorsys = wcscsys (cstr); if (!tabhgetr8 (startab,"equinox", &sc->equinox)) sc->equinox = wcsceq (cstr); if (!tabhgetr8 (startab,"epoch",&sc->epoch)) sc->epoch = sc->equinox; } else if (tabhgetr8 (startab,"EQUINOX", &sc->equinox)) { if (!tabhgetr8 (startab,"EPOCH",&sc->epoch)) sc->epoch = sc->equinox; if (sc->equinox == 1950.0) sc->coorsys = WCS_B1950; else sc->coorsys = WCS_J2000; } else if (tabhgetr8 (startab,"equinox", &sc->equinox)) { if (!tabhgetr8 (startab,"epoch",&sc->epoch)) sc->epoch = sc->equinox; if (sc->equinox == 1950.0) sc->coorsys = WCS_B1950; else sc->coorsys = WCS_J2000; } else if (tabhgetr8 (startab,"EPOCH", &sc->epoch)) { sc->equinox = sc->epoch; if (sc->equinox == 1950.0) sc->coorsys = WCS_B1950; else sc->coorsys = WCS_J2000; } else if (tabhgetr8 (startab,"epoch", &sc->epoch)) { sc->equinox = sc->epoch; if (sc->equinox == 1950.0) sc->coorsys = WCS_B1950; else sc->coorsys = WCS_J2000; } /* Check whether catalog is sorted by right ascension */ if (tabhgetc (startab, "rasort", cstr)) sc->rasorted = 1; else sc->rasorted = 0; /* Set other stuff */ sc->nstars = startab->nlines; if (sc->entid) sc->stnum = 1; else sc->stnum = 0; /* Find out if ID is number, and if so, how many decimal places it has */ if (sc->entid) { istar = 1; if ((line = gettabline (startab, istar)) == NULL) { fprintf (stderr,"TABCATOPEN: Cannot read first star\n"); tabcatclose (sc); return (NULL); } if (setoken (&startok, line, "tab") < 2) { fprintf (stderr,"TABCATOPEN: First star has too few columns\n"); tabcatclose (sc); return (NULL); } tabgetc (&startok, sc->entid, cstr, 32); /* Find length of identifier */ if (tabhgeti4 (startab,"nfield", &nnfld)) sc->nnfld = nnfld; else sc->nnfld = strlen (cstr); /* Find number of decimal places in identifier */ if (tabhgeti4 (startab, "ndec", &nndec)) { sc->nndec = nndec; if (!isnum (cstr)) sc->stnum = -nndec; } else if (isnum (cstr)) { dnum = tabgetr8 (&startok,sc->entid); sprintf (cstr,"%.0f", (dnum * 100000000.0) + 0.1); lnum = strlen (cstr); for (i = 0; i < 8; i++) { if (cstr[lnum-i-1] != '0') { ndec = 8 - i; if (ndec > nndec) { nndec = ndec; sc->nndec = nndec; } break; } } sc->nndec = nndec; } else { sc->stnum = -strlen (cstr); sc->nndec = nndec; } } sc->refcat = TABCAT; return (sc); } /* TABCATCLOSE -- Close tab table catalog and free associated data structures */ void tabcatclose (sc) struct StarCat *sc; { tabclose (sc->startab); free (sc); return; } /* TABSTAR -- Get tab table catalog entry for one star; return 0 if successful, else -1 */ int tabstar (istar, sc, st, verbose) int istar; /* Star sequence number in tab table catalog */ struct StarCat *sc; /* Star catalog data structure */ struct Star *st; /* Star data structure, updated on return */ int verbose; /* 1 to print error messages */ { struct TabTable *startab = sc->startab; char *line; char *uscore; char cnum[64]; char temp[64]; double ydate; char *cn; int ndec, i, imag, ltok; int lnum, ireg, inum; int pmq, nid; int lcn; char str[24]; if ((line = gettabline (startab, istar)) == NULL) { if (verbose) fprintf (stderr,"TABSTAR: Cannot read star %d\n", istar); return (-1); } /* Parse line for rapid field extraction */ if (setoken (&startok, line, "tab") < 2) { fprintf (stderr,"TABSTAR: Star %d entry too short\n", istar); return (-1); } /* Extract ID */ st->objname[0] = (char) 0; if (sc->entid) { tabgetc (&startok, sc->entid, cnum, 64); if (!strcmp (sc->isfil,"usnoa-server")) { if ((uscore = strchr (cnum, '_')) != NULL) *uscore = '.'; cn = cnum + 1; } else if (!strcmp (sc->isfil,"gsc-server")) cn = cnum + 3; else cn = cnum; lcn = strlen (cn); /* if (lcn < 16 && (isnum (cn) || isnum (cn+1))) { */ if (lcn < 16 && isnum (cn)) { if (isnum(cnum)) { lcn = strlen (cn); if (lcn > 15) st->num = atof (cn + lcn - 15); else st->num = atof (cn); } else { if (cnum[0] == 'S') st->num = -atof (cn+1); else st->num = atof (cn+1); } if (!strcmp (sc->isfil,"gsc-server")) { ireg = atoi (cn) / 100000; inum = atoi (cn) % 100000; st->num = (double) ireg + 0.0001 * (double) inum; } /* Find field length of identifier */ sc->nnfld = strlen (cn); /* Find number of decimal places in identifier */ if (strchr (cn,'.') == NULL) { nndec = 0; sc->nndec = nndec; } else { sprintf (cn,"%.0f", (st->num * 100000000.0) + 0.1); lnum = strlen (cnum); for (i = 0; i < 8; i++) { if (cn[lnum-i-1] != '0') { ndec = 8 - i; if (ndec > nndec) { nndec = ndec; sc->nndec = nndec; } break; } } } } else { strcpy (st->objname, cnum); st->num = st->num + 1.0; } } else { st->num = (double) istar; nndec = 0; sc->nndec = nndec; } /* Right ascension */ st->ra = tabgetra (&startok, sc->entra); /* Declination */ st->dec = tabgetdec (&startok, sc->entdec); /* Magnitudes */ for (imag = 0; imag < sc->nmag; imag++) { if (sc->entmag[imag]) { if (tabgetc (&startok, sc->entmag[imag], str, 24)) return (0); ltok = strlen (str); if (str[ltok-1] == 'L') { str[ltok-1] = (char) 0; if (isnum (str)) st->xmag[imag] = 100.0 + atof (str); else st->xmag[imag] = 0.0; } else st->xmag[imag] = tabgetr8 (&startok, sc->entmag[imag]); } else st->xmag[imag] = 0.0; } /* Convert right ascension proper motion to degrees/year */ st->rapm = tabgetr8 (&startok, sc->entrpm); if (sc->rpmunit == PM_MASYR) st->rapm = (st->rapm / 3600000.0) / cosdeg(st->dec); else if (sc->rpmunit == PM_MTSYR) st->rapm = st->rapm / 240000.0; else if (sc->rpmunit == PM_ARCSECYR) st->rapm = (st->rapm / 3600.0) / cosdeg (st->dec); else if (sc->rpmunit == PM_ARCSECCEN) st->rapm = (st->rapm / 360000.0) / cosdeg (st->dec); else if (sc->rpmunit == PM_TSECYR) st->rapm = st->rapm / 240.0; else if (sc->rpmunit == PM_TSECCEN) st->rapm = st->rapm / 24000.0; else if (sc->rpmunit == PM_RADYR) st->rapm = raddeg (st->rapm); else st->rapm = 0.0; /* Convert declination proper motion to degrees/year */ st->decpm = tabgetr8 (&startok, sc->entdpm); if (sc->dpmunit == PM_MASYR) st->decpm = st->decpm / 3600000.0; else if (sc->dpmunit == PM_ARCSECYR) st->decpm = st->decpm / 3600.0; else if (sc->dpmunit == PM_ARCSECCEN) st->decpm = st->decpm / 360000.0; else if (sc->dpmunit == PM_RADYR) st->decpm = raddeg (st->decpm); else st->decpm = 0.0; /* Parallax */ if (sc->entpx) st->parallax = tabgetr8 (&startok, sc->entpx); else st->parallax = 0.0; /* Radial velocity */ if (sc->entrv) st->radvel = tabgetr8 (&startok, sc->entrv); else st->radvel = 0.0; /* Epoch */ if (sc->entepoch) { tabgetc (&startok, sc->entepoch, temp, 10); if (temp[0] == '_') { if (sc->entdate > 0) { tabgetc (&startok, sc->entdate, temp, 10); st->epoch = fd2ep (temp); } else st->epoch = sc->epoch; } else if (strchr (temp, '-') != NULL) st->epoch = fd2ep (temp); else { ydate = tabgetr8 (&startok, sc->entepoch); if (ydate < 3000.0 && ydate > 0.0) st->epoch = dt2ep (ydate, 12.0); else if (ydate < 100000.0) st->epoch = mjd2ep (ydate); else st->epoch = jd2ep (ydate); if (st->epoch > 2005.000) printf ("TABSTAR: %s = %.5f -> %.5f\n", temp, ydate, st->epoch); } st->xmag[sc->nmag-1] = st->epoch; } else st->epoch = sc->epoch; /* Peak counts */ if (sc->entpeak > 100) { pmq = tabgeti4 (&startok, sc->entpeak/100); nid = tabgeti4 (&startok, sc->entpeak%100); st->peak = (pmq * 100) + nid; } else if (sc->entpeak > 0) st->peak = tabgeti4 (&startok, sc->entpeak); else st->peak = 0; /* Spectral type */ if (sc->enttype > 0) { strcpy (st->isp, "__"); tabgetc (&startok, sc->enttype, st->isp, 24); } /* Extract selected field */ if (kwo != NULL) (void) tabgetk (startab, &startok, kwo, st->objname, 79); st->coorsys = sc->coorsys; st->equinox = sc->equinox; return (0); } /* TABOPEN -- Open tab table file, returning number of entries */ struct TabTable * tabopen (tabfile, nbbuff) char *tabfile; /* Tab table catalog file name */ int nbbuff; /* Number of bytes in buffer; 0=read whole file */ { FILE *fcat; int nr, lfile, lname, lline; char *tabnew, *tabline, *lastline; char *tabcomma, *nextline; char *thisname, *tabname; char tabcstr[4]; int thistab, itab, nchar, nbtab; int formfeed = (char) 12; struct TabTable *tabtable; tabcomma = NULL; if (taberr != NULL) { free (taberr); taberr = NULL; } tabname = NULL; if (!strcmp (tabfile, "stdin")) { lfile = 100000; fcat = stdin; } else { /* Separate table name from file name, if necessary */ if ((tabcomma = strchr (tabfile, ',')) != NULL) { tabname = (char *) calloc (1,64); strcpy (tabname, tabcomma+1); *tabcomma = (char) 0; } /* Find length of tab table catalog */ lfile = tabsize (tabfile); if (nbbuff > 1 && nbbuff < lfile) lfile = nbbuff; else if (nbbuff == 1 && lfile > 10000) lfile = 10000; if (lfile < 1) { taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: Tab table file %s has no entries", tabfile); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } /* Open tab table catalog */ if (!(fcat = fopen (tabfile, "r"))) { taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: Tab table file %s cannot be read", tabfile); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else if (verbose) { fprintf (stderr,"TABOPEN: tab table %s opened", tabfile); } } /* Allocate tab table structure */ nbtab = sizeof(struct TabTable); if ((tabtable=(struct TabTable *) calloc(1,nbtab)) == NULL){ taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: cannot allocate %d bytes for tab table structure for %s", nbtab, tabfile); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else if (verbose) { fprintf (stderr,"TABOPEN: allocated %d bytes for tab table structure for %s", nbtab, tabfile); } tabtable->tabname = tabname; /* Allocate space in structure for filename and save it */ lname = strlen (tabfile) + 2; if ((tabtable->filename = (char *)calloc (1, lname)) == NULL) { taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: cannot allocate filename %s in structure", tabfile); (void) fclose (fcat); tabclose (tabtable); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } strcpy (tabtable->filename, tabfile); /* Allocate buffer to hold entire catalog (or buffer length) and read it */ if ((tabtable->tabbuff = (char *) calloc (1, lfile+2)) == NULL) { taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: cannot allocate buffer for tab table %s", tabfile); (void) fclose (fcat); tabclose (tabtable); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else { if (verbose) { fprintf (stderr,"TABOPEN: allocated %d bytes for tab table for %s", lfile+2, tabfile); } nr = fread (tabtable->tabbuff, 1, lfile, fcat); if (fcat != stdin && nr < lfile) { fprintf (stderr,"TABOPEN: read only %d / %d bytes of file %s\n", nr, lfile, tabfile); (void) fclose (fcat); tabclose (tabtable); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else if (verbose) { fprintf (stderr,"TABOPEN: read %d byte tab table from %s", lfile, tabfile); } tabtable->tabbuff[lfile] = (char) 0; /* Check for named table within a file */ if (tabname != NULL) { if (isnum (tabname)) { itab = atoi (tabname); thisname = tabtable->tabbuff; thistab = 1; if (itab > 1) { while (thistab < itab && thisname != NULL) { thisname = strchr (thisname, formfeed); if (thisname != NULL) thisname++; thistab++; } } if (thisname == NULL) { fprintf (stderr, "GETTAB: There are < %d tables in %s\n", itab, tabfile); return (NULL); } while (*thisname==' ' || *thisname==newline || *thisname==formfeed || *thisname==(char)13) thisname++; tabline = strchr (thisname, newline); if (tabline != NULL) { nchar = tabline - thisname; if (strchr (thisname, tab) > tabline) strncpy (tabtable->tabname, thisname, nchar); } } else { lname = strlen (tabname); thisname = tabtable->tabbuff; while (*thisname != (char) 0) { while (*thisname==' ' || *thisname==newline || *thisname==formfeed || *thisname==(char)13) thisname++; if (!strncmp (tabname, thisname, lname)) break; else thisname = strchr (thisname, formfeed); } } if (thisname == NULL) { fprintf (stderr, "TABOPEN: table %s in file %s not found\n", tabname, tabfile); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else tabtable->tabheader = strchr (thisname, newline) + 1; } else tabtable->tabheader = tabtable->tabbuff; tabline = tabtable->tabheader; lastline = NULL; while (*tabline!='-' && tabline!=NULL && tabline < tabtable->tabbuff+lfile) { lastline = tabline; if ((tabline = strchr (lastline,newline)) != NULL) tabline++; else break; } if (tabline == NULL || *tabline != '-') { taberr = (char *) calloc (64 + strlen (tabfile), 1); sprintf (taberr,"TABOPEN: No - line in tab table %s",tabfile); (void) fclose (fcat); tabclose (tabtable); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } tabcstr[0] = tabline[1]; tabcstr[1] = (char) 0; if ( tabline[1] != '-' && tabline[1] != (char)9 ) { sprintf (taberr,"TABOPEN: No - line in tab table %s",tabfile); return (NULL); } tabtable->tabhead = lastline; tabtable->tabdash = tabline; tabtable->tabdata = strchr (tabline, newline) + 1; /* Extract positions of keywords we will want to use */ if (!tabparse (tabtable)) { fprintf (stderr,"TABOPEN: No columns in tab table %s\n",tabfile); (void) fclose (fcat); tabclose (tabtable); if (tabcomma != NULL) *tabcomma = ','; return (NULL); } else if (verbose) { fprintf (stderr,"TABOPEN: tab table %s header parsed", tabfile); } /* Enumerate entries in tab table catalog by counting newlines */ if (nbbuff > 0) tabtable->nlines = 10000000; else { tabnew = tabtable->tabdata; tabtable->nlines = 0; while ((tabnew = strchr (tabnew, newline)) != NULL) { tabnew = tabnew + 1; tabtable->nlines = tabtable->nlines + 1; if (*tabnew == formfeed) break; } } if (verbose) { fprintf (stderr,"TABOPEN: %d lines in tab table %s", tabtable->nlines, tabfile); } } /* Set up line buffer, and put first line in it */ if (nbbuff > 0) { tabtable->lhead = tabtable->tabdata - tabtable->tabheader; tabtable->tcat = fcat; nextline = strchr (tabtable->tabdata, newline) + 1; tabtable->lline = (nextline - tabtable->tabdata) * 2; tabtable->tabline = (char *) calloc (tabtable->lline, 1); fseek (tabtable->tcat, (long) tabtable->lhead, SEEK_SET); (void) fgets (tabtable->tabline, tabtable->lline, tabtable->tcat); lline = strlen (tabtable->tabline); if (tabtable->tabline[lline-1] < 32) tabtable->tabline[lline-1] = (char) 0; tabtable->tabdata = tabtable->tabline; } /* Close catalog file if not reading one line at a time */ else { tabtable->lhead = 0; (void) fclose (fcat); tabtable->tcat = NULL; } tabtable->tabline = tabtable->tabdata; tabtable->iline = 1; if (tabcomma != NULL) *tabcomma = ','; return (tabtable); } void tabclose (tabtable) struct TabTable *tabtable; { if (tabtable != NULL) { if (tabtable->filename != NULL) free (tabtable->filename); if (tabtable->tabname != NULL) free (tabtable->tabname); if (tabtable->tabbuff != NULL) free (tabtable->tabbuff); if (tabtable->colname != NULL) free (tabtable->colname); if (tabtable->lcol != NULL) free (tabtable->lcol); if (tabtable->lcfld != NULL) free (tabtable->lcfld); if (tabtable->tcat != NULL) fclose (tabtable->tcat); free (tabtable); } return; } /* TABLINE -- Get tab table entry for one line; return NULL if unsuccessful */ char * gettabline (tabtable, iline) struct TabTable *tabtable; /* Tab table structure */ int iline; /* Line sequence number in tab table */ { char *nextline = tabtable->tabline; char *next; int lline, i; /* Return NULL if tab table has not been opened */ if (tabtable == NULL) return (NULL); /* Read one line at a time if file is still open */ if (tabtable->tcat != NULL) { /* Return current line from buffer */ if (iline == tabtable->iline) return (tabtable->tabline); /* Read next line from file */ if (iline < 1 || iline > tabtable->iline) { for (i = tabtable->iline; i < iline; i++) { next = fgets (tabtable->tabline, tabtable->lline, tabtable->tcat); if (next == NULL || *next == EOF) return (NULL); tabtable->iline++; } lline = strlen (tabtable->tabline); if (lline < 2) return (NULL); if (tabtable->tabline[lline-1] < 32) tabtable->tabline[lline-1] = (char) 0; } else if (iline < tabtable->iline) { fseek (tabtable->tcat, (long) tabtable->lhead, SEEK_SET); tabtable->iline = 0; for (i = tabtable->iline; i < iline; i++) { (void) fgets (tabtable->tabline, tabtable->lline, tabtable->tcat); tabtable->iline++; } lline = strlen (tabtable->tabline); if (tabtable->tabline[lline-1] < 32) tabtable->tabline[lline-1] = (char) 0; } return (tabtable->tabline); } /* Return NULL if trying to read past last line */ if (iline > tabtable->nlines) { fprintf (stderr, "TABLINE: line %d is not in table\n",iline); return (NULL); } /* If iline is 0 or less, just read next line from table */ else if (iline < 1 && nextline) { tabtable->iline++; if (tabtable->iline > tabtable->nlines) { fprintf (stderr, "TABLINE: line %d is not in table\n",iline); return (NULL); } nextline = strchr (nextline, newline) + 1; } /* If iline is before current line, read from start of file */ else if (iline < tabtable->iline) { tabtable->iline = 1; tabtable->tabline = tabtable->tabdata; while (tabtable->iline < iline) { tabtable->tabline = strchr (tabtable->tabline, newline) + 1; tabtable->iline ++; } } /* If iline is after current line, read forward */ else if (iline > tabtable->iline) { while (tabtable->iline < iline) { tabtable->tabline = strchr (tabtable->tabline, newline) + 1; tabtable->iline ++; } } return (tabtable->tabline); } /* TABGETRA -- returns double right ascension in degrees */ double tabgetra (tabtok, ientry) struct Tokens *tabtok; /* Line token structure */ int ientry; /* sequence of entry on line */ { char str[24]; strcpy (str, "0.0"); if (tabgetc (tabtok, ientry, str, 24)) return (0.0); else return (str2ra (str)); } /* TABGETDEC -- returns double declination in degrees */ double tabgetdec (tabtok, ientry) struct Tokens *tabtok; /* Line token structure */ int ientry; /* sequence of entry on line */ { char str[24]; strcpy (str, "0.0"); if (tabgetc (tabtok, ientry, str, 24)) return (0.0); else return (str2dec (str)); } /* TABGETR8 -- returns 8-byte floating point number from tab table line */ double tabgetr8 (tabtok, ientry) struct Tokens *tabtok; /* Line token structure */ int ientry; /* sequence of entry on line */ { char str[24]; strcpy (str, "0.0"); if (tabgetc (tabtok, ientry, str, 24)) return (0.0); else if (isnum (str)) return (atof (str)); else return (0.0); } /* TABGETI4 -- returns a 4-byte integer from tab table line */ int tabgeti4 (tabtok, ientry) struct Tokens *tabtok; /* Line token structure */ int ientry; /* sequence of entry on line */ { char str[24]; strcpy (str, "0"); if (tabgetc (tabtok, ientry, str, 24)) return (0); else if (isnum (str)) return ((int) atof (str)); else return (0); } /* TABGETK -- returns a character entry from tab table line for named column */ int tabgetk (tabtable, tabtok, keyword, string, maxchar) struct TabTable *tabtable; /* Tab table structure */ struct Tokens *tabtok; /* Line token structure */ char *keyword; /* column header of desired value */ char *string; /* character string (returned) */ int maxchar; /* Maximum number of characters in returned string */ { int ientry = tabccol (tabtable, keyword); return (tabgetc (tabtok, ientry, string, maxchar)); } /* TABGETC -- returns n'th entry from tab table line as character string */ int tabgetc (tabtok, ientry, string, maxchar) struct Tokens *tabtok; /* Line token structure */ int ientry; /* Sequence of entry on line (1-ncol) */ char *string; /* Character string (returned) */ int maxchar; /* Maximum number of characters in returned string */ { if (ientry > tabtok->ntok) return (0); else if (getoken (tabtok, ientry, string, maxchar)) return (0); else return (-1); } /* TABHGETR8 -- read an 8-byte floating point number from a tab table header */ static int tabhgetr8 (tabtable, keyword, result) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* sequence of entry on line */ double *result; { char value[24]; if (tabhgetc (tabtable, keyword, value)) { *result = atof (value); return (1); } else return (0); } /* TABHGETI4 -- read a 4-byte integer from a tab table header */ static int tabhgeti4 (tabtable, keyword, result) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* sequence of entry on line */ int *result; { char value[24]; if (tabhgetc (tabtable, keyword, value)) { *result = (int) atof (value); return (1); } else return (0); } /* TABHGETC -- read a string from a tab table header */ static int tabhgetc (tabtable, keyword, result) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* sequence of entry on line */ char *result; { char *str0, *str1, *line, *head, keylow[24], keyup[24]; int ncstr, lkey, i; head = tabtable->tabbuff; str0 = 0; /* Make all-upper-case and all-lower-case versions of keyword */ lkey = strlen (keyword); if (lkey > 24) lkey = 24; for (i = 0; i < lkey; i++) { if (keyword[i] > 96 && keyword[i] < 123) keyup[i] = keyword[i] - 32; else keyup[i] = keyword[i]; if (keyword[i] > 64 && keyword[i] < 91) keylow[i] = keyword[i] + 32; else keylow[i] = keyword[i]; } keyup[lkey] = (char) 0; keylow[lkey] = (char) 0; /* Find keyword or all-upper-case or all-lower-case version in header */ while (head < tabtable->tabhead) { line = strsrch (head, keyword); if (line == NULL) line = strsrch (head, keylow); if (line == NULL) line = strsrch (head, keyup); if (line == NULL) break; if (line == tabtable->tabbuff || line[-1] == newline) { str0 = strchr (line, tab) + 1; str1 = strchr (str0, newline); break; } else head = line + 1; } /* Return value as a character string and 1 if found */ if (str0) { ncstr = str1 - str0; strncpy (result, str0, ncstr); result[ncstr] = (char)0; return (1); } else return (0); } /* TABPARSE -- Make a table of column headings */ int tabparse (tabtable) struct TabTable *tabtable; /* Tab table structure */ { char *colhead; /* Column heading first character */ char *endcol; /* Column heading last character */ char *headlast; char *hyphens; char *hyphlast; char *nextab; int icol; int nbytes, nba; /* Return if no column names in header */ if (tabtable->tabhead == NULL) return (0); headlast = strchr (tabtable->tabhead, newline); if (headlast == tabtable->tabhead) return (0); /* Count columns in table header */ tabtable->ncols = 1; for (colhead = tabtable->tabhead; colhead < headlast; colhead++) { if (*colhead == tab) tabtable->ncols++; } /* Tabulate column names */ nbytes = tabtable->ncols * sizeof (char *); nba = nbytes / 64; if (nbytes > nba * 64) nba = (nba + 1) * 64; else nba = nba * 64; tabtable->colname = (char **)calloc (tabtable->ncols, sizeof (char *)); tabtable->lcol = (int *) calloc (tabtable->ncols, sizeof (int)); colhead = tabtable->tabhead; for (icol = 0; icol < tabtable->ncols; icol++) { nextab = strchr (colhead, tab); if (nextab < headlast) endcol = nextab - 1; else endcol = headlast - 1; while (*endcol == ' ') endcol = endcol - 1; tabtable->lcol[icol] = (int) (endcol - colhead) + 1; tabtable->colname[icol] = colhead; colhead = nextab + 1; if (colhead > headlast) break; } /* Tabulate field widths */ hyphens = headlast + 1; hyphlast = strchr (hyphens, newline); if (hyphlast == hyphens) return (0); tabtable->lcfld = (int *) calloc (tabtable->ncols, sizeof (int)); colhead = hyphens; for (icol = 0; icol < tabtable->ncols; icol++) { if ((nextab = strchr (colhead, tab)) == NULL) endcol = hyphlast - 1; else endcol = nextab - 1; tabtable->lcfld[icol] = (int) (endcol - colhead) + 1; if (nextab != NULL) colhead = nextab + 1; else break; } return (tabtable->ncols); } /* Search table of column headings for a particlar entry (case-dependent) */ int tabcol (tabtable, keyword) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* Column heading to find */ { int i, lkey, lcol; lkey = strlen (keyword); for (i = 0; i < tabtable->ncols; i++) { lcol = tabtable->lcol[i]; if (lcol == lkey && !strncmp (keyword, tabtable->colname[i], lcol)) { return (i + 1); } } return (0); } /* Search table of column headings for a particlar entry (case-independent) */ int tabccol (tabtable, keyword) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* Column heading to find */ { int i, lkey, lcol; lkey = strlen (keyword); for (i = 0; i < tabtable->ncols; i++) { lcol = tabtable->lcol[i]; if (lcol == lkey && !strncasecmp (keyword, tabtable->colname[i], lcol)) { return (i + 1); } } return (0); } /* Search table of column headings for first with string (case-dependent) */ static int tabcont (tabtable, keyword) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* Part of column heading to find */ { int i; for (i = 0; i < tabtable->ncols; i++) { if (strnsrch (tabtable->colname[i], keyword, tabtable->lcol[i])) { return (i + 1); } } return (0); } /* Search table of column headings for first with string (case-independent) */ static int tabccont (tabtable, keyword) struct TabTable *tabtable; /* Tab table structure */ char *keyword; /* Part of column heading to find */ { int i; for (i = 0; i < tabtable->ncols; i++) { if (strncsrch (tabtable->colname[i], keyword, tabtable->lcol[i])) { return (i + 1); } } return (0); } /* TABSIZE -- return size of file in bytes */ static int tabsize (filename) char *filename; /* Name of file for which to find size */ { FILE *diskfile; long filesize; /* Open file */ if ((diskfile = fopen (filename, "r")) == NULL) return (-1); /* Move to end of the file */ if (fseek (diskfile, 0, 2) == 0) /* Position is the size of the file */ filesize = ftell (diskfile); else filesize = -1; fclose (diskfile); return (filesize); } /* ISTAB -- Return 1 if tab table file, else 0 */ int istab (filename) char *filename; /* Name of file to check */ { FILE *diskfile; char *line; int ncmax=63; /* First check file extension */ if (strsrch (filename, ".tab")) return (1); /* If no .tab file extension, look for tab on first line */ else { /* Open file */ if ((diskfile = fopen (filename, "r")) == NULL) return (0); /* Allocate buffer for first line of file */ line = (char *)calloc (1, 63); /* Return 0 if no characters in first line */ if (!next_line (diskfile, ncmax, line)) { fclose (diskfile); free (line); return (0); } /* Return 0 if there is a tab in the first line */ if (strchr (line, (char)9)) { fclose (diskfile); free (line); return (1); } /* Otherwise, return 0 */ else { fclose (diskfile); free (line); return (0); } } } /* Jul 18 1996 New subroutines * Aug 6 1996 Remove unused variables after lint * Aug 8 1996 Fix bugs in entry reading and logging * Oct 15 1996 Add comparison when testing an assignment * Nov 5 1996 Drop unnecessary static declarations * Nov 13 1996 Return no more than maximum star number * Nov 13 1996 Write all error messages to stderr with subroutine names * Nov 15 1996 Implement search radius; change input arguments * Nov 19 1996 Allow lower case column headings * Dec 18 1996 Add UJCRNUM to read specified catalog entries * Dec 18 1996 Keep closest stars, not brightest, if searching within radius * * Mar 20 1997 Clean up code in TABRNUM * May 7 1997 Set entry number to zero if column not found * May 29 1997 Add TABPARSE and TABCOL to more easily extract specific columns * May 29 1997 Add TABCLOSE to free memory from outside this file * Jun 4 1997 Set ID to sequence number in table if no ID/id entry present * * Jun 2 1998 Fix bug parsing last column of header * Jun 24 1998 Add string lengths to ra2str() and dec2str() calls * Sep 16 1998 Use limiting radius correctly * Sep 22 1998 Convert to output coordinate system * Oct 15 1998 Add tabsize() and istab() * Oct 21 1998 Add tabrkey() to read values of keyword for list of stars * Oct 29 1998 Correctly assign numbers when too many stars are found * Oct 30 1998 Fix istab() to check only first line of file * Dec 8 1998 Do not declare tabsize() static * Jan 20 1999 Add tabcatopen() and keep table info in structure, not global * Jan 25 1999 Add lcfld to structure to keep track of field widths * Jan 29 1999 Add current line number and pointer to table structure * Feb 2 1999 Allow for equinox other than 2000.0 in tab table header * Feb 2 1999 Add tabhgetc() to read char string values from tab table header * Feb 17 1999 Increase maximum line in istab() from 80 to 1024 * Mar 2 1999 Fix bugs calling tabhgetx() * Mar 2 1999 Rewrite tabhgetx() to use tabhgetc() for all header reading * May 28 1999 Add tabcatopen() and tabstar() and use them * Jun 3 1999 Fix bug so header parameters are read correctly * Jun 16 1999 Use SearchLim() * Aug 16 1999 Fix bug to fix failure to search across 0:00 RA * Aug 25 1999 Return real number of stars from tabread() * Sep 10 1999 Fix bug setting equinox and coordinate system in tabcatopen() * Sep 10 1999 Set additional keyword selection with subroutine * Sep 13 1999 Fix comment for tabstar() * Sep 16 1999 Fix bug which didn't always return closest stars * Sep 16 1999 Add distsort argument so brightest stars in circle works, too * Oct 21 1999 Clean up code after lint * Oct 25 1999 Fix subroutine declaration inconsistency * Oct 25 1999 Replace malloc() calls with calloc() * Oct 29 1999 Add tabxyread() for image catalogs * Nov 23 1999 Improve error checking on Starbase tables; istab() opens file * Nov 30 1999 Fix bugs found when compiling under SunOS 4.1.3 * * Jan 4 2000 Always close file and print error message on tabopen() failures * Jan 6 2000 If "id" not found, try heading with "_id" to catch scat output * Jan 10 2000 Add second magnitude; save column headers in catalog structure * Feb 10 2000 Implement proper motion in source catalogs * Feb 10 2000 Accept first mag-containing column as first magnitude * Feb 10 2000 Clean up id reading: no. decimals, non-numeric id * Feb 14 2000 Save table opening errors in string * Feb 16 2000 Lengthen short calloc() lengths * Feb 16 2000 Pad tabbuff with 2 nulls so end can be found * Mar 10 2000 Return proper motions from tabread() and tabrnum() * Mar 13 2000 Do not free tabtable structure if it is null * Mar 27 2000 Clean up code after lint * May 26 2000 Add ability to read named tables in a multi-table file * Jun 26 2000 Add coordinate system to SearchLim() arguments * Jul 5 2000 Check for additional column heading variations * Jul 10 2000 Deal with number of decimal places and name/number in tabcatopen() * Jul 12 2000 Add star catalog data structure to tabread() argument list * Jul 13 2000 Use nndec header parameter to optionally set decimal places * Jul 25 2000 Pass star catalog address of data structure address * Aug 3 2000 Skip first character of ID if rest is number * Aug 3 2000 If no decimal point in numeric ID, set ndec to zero * Sep 27 2000 Use first column with name if no id column is found * Oct 26 2000 Add proper motion in seconds and arcseconds per century * Oct 31 2000 Add proper motion in milliseconds of time per year * Nov 22 2000 Add tabtable argument to tabcatopen() for URL search returns * Nov 28 2000 Add starcat as argument to tabrnum() as well as tabread() * Nov 29 2000 Do not set tmagb if it is null; set type if present * Nov 30 2000 Add spectral type as possible column * Dec 1 2000 Add field as synonym for plate * Dec 18 2000 Clean up after lint * Dec 29 2000 Clean up after lint * * May 29 2001 Save length of identifier in catalog structure * May 30 2001 Rewrite to deal with up to 3 magnitudes * Jun 11 2001 Add one-line-at-a-time catalog reading from tabcatopen() * Jun 14 2001 In tabopen(), use actual file size if less than buffer size * Jun 14 2001 In tabcol() make sure lengths are equal, too * Jun 18 2001 Fix bug finding of length of token string extracted * Jun 19 2001 Use setoken() to parse each catalog line * Jun 20 2001 Add fourth magnitude for GSC II * Jun 25 2001 Check fgets() return in gettabline() * Jun 25 2001 Print star read errors only in verbose mode * Jun 28 2001 If no proper motion unit is set, set proper motion to 0.0 * Jun 28 2001 If first magnitude is 99.90, sort by second magnitude * Jun 28 2001 Up default line limit from 1 million to 10 million * Jul 2 2001 Fix order of 3rd and 4th magnitudes * Jul 20 2001 Return magnitude as well as flux from tabxyread() * Aug 8 2001 Return radial velocity as additional magnitude * Aug 17 2001 Fix bug reading radial velocity * Aug 21 2001 Check numbers using isnum() in tabgetr8() and tabgeti4() * Aug 21 2001 Add starcat to tabrkey() argument list * Aug 21 2001 Read object name into tkey if it is present * Sep 11 2001 Allow an arbitrary number of magnitudes * Sep 11 2001 Add sort magnitude argument to tabread() * Sep 19 2001 Drop fitshead.h; it is in wcs.h * Sep 20 2001 Clean up pointers for Alpha and Linux; drop tabgetpm() * Sep 20 2001 Change _ to . in ESO server USNO-A2.0 numbers * Sep 21 2001 Rearrange ESO server GSC numbers * Oct 12 2001 Check for additional column names for magnitude * Oct 15 2001 Read first star only once, compute relative flux in tabxyread() (bug fix) * Oct 16 2001 Add pointer to line of dashes to table structure * Dec 3 2001 Initialize keyword search variables in tabhgetc() * * May 6 2002 Allow object names to be up to 79 characters long * Aug 5 2002 Deal correctly with magnitude-less catalogs * Aug 5 2002 Add magu, magb, magv for UBV magnitudes * Aug 6 2002 Pass through magnitude keywords * Aug 6 2002 Return initial string if token not found by tabgetc() * Oct 30 2002 Add code to pass epoch as a final magnitude (but not RV, yet) * * Jan 26 2003 Add code to pass USNO-B1.0 pm quality and no. of ids * Jan 28 2003 Improve spatial position test * Feb 4 2003 Compare character to 0, not NULL * Mar 11 2003 Fix limit setting * Apr 3 2003 Drop unused variables after lint * May 28 2003 Read long and lat if radecsys is "galactic" or "ecliptic" * Jun 2 003 Divide by cos(dec) for arcsec and mas RA proper motions * Aug 22 2003 Add radi argument for inner edge of search annulus * Sep 25 2003 Add tabbin() to fill an image with sources * Oct 6 2003 Update tabread() and tabbin() for improved RefLim() * Nov 18 2003 Initialize image size and bits/pixel from header in tabbin() * Nov 22 2003 Add GSC II object class (c) as possible content for tpeak * Dec 4 2003 Add default of arcsec/century as proper motion unit * Dec 10 2003 If magnitude ends in L and is a number, add 100.0 * Dec 12 2003 Fix bug in wcs2pix() call in tabbin() * * Jan 5 2004 If more than 15 digits in numberic ID, drop excess off front * Mar 16 2004 Be more clever about reading by number * Mar 19 2004 Make verbose flag global * Nov 17 2004 Accept SpT and spt before type for spectral type * * May 12 2005 Add "num" as a possible substring to identify an ID column * Jul 26 2005 Fix bug to set magnitudes and keywords correctly in tabrnum() * Aug 3 2005 If nlog < 0 in tabrnum(), print objects as found * Aug 10 2005 Add "exp" as possible "magnitude" column heading * Aug 11 2005 Do not re-read lines if missing object numbers in tabrnum() * Aug 17 2005 If nmax is -1, print magnitude keywords from file * Sep 29 2005 Read first 10000 characters for istab() to capture long heads * Sep 29 2005 Add rasort header flag if catalog is sorted by RA * * Jun 15 2006 Read ID field length from header, if present * Jun 20 2006 Drop unused variables; initialize uninitialized variables * Jun 30 2006 Add match argument to tabrnum() to add sequential read option * * Jan 11 2007 Switch order of wcscat.h and fitsfile.h includes * Mar 13 2007 Return object name if first character of ID is non-numeric * Mar 13 2007 Check for "Class" for GSC2 object class * Jul 9 2007 Reject mflag as a magnitude * Jul 9 2007 Sort by designated magnitude only unless sort mag is 0 * Jul 18 2007 Add tabccol() and tabccont() case-insensitive header searches * Jul 19 2007 Add proper motion in arsec/hour for solar system objects * Jul 23 2007 Add ...obj... as possible identifier * Jul 23 2007 Add ...dist... as possible "magnitude" * * Aug 17 2009 Fix columns for declination column name * Sep 25 2009 Fix memory leaks found by Jessicalas Burke * Sep 30 2009 Fix bugs freeing object names for first pass and farthest star * * Sep 22 2010 Fix bug when checking file which is not tab table * * Jun 20 2014 Not a tab table if no tab in first line * Aug 28 2014 Fix istab() and tabread() * Aug 29 2014 Exclude "code" from magnitude keywords */