/*** File libwcs/ctgread.c *** September 27, 2001 *** By Doug Mink, dmink@cfa.harvard.edu *** Harvard-Smithsonian Center for Astrophysics */ /* int ctgread() Read catalog stars in specified region of the sky * int ctgrnum() Read catalog stars with specified numbers * int ctgopen() Open catalog, returning number of entries * int ctgstar() Get ASCII catalog entry for one star * int ctgsize() Return length of file in bytes * int isacat() Return 1 if file is ASCII catalog, else 0 */ #include #include #include #include #include #include #include "wcs.h" #include "wcscat.h" /* default pathname for catalog, used if catalog file not found in current working directory, but overridden by WCS_CATDIR environment variable */ char catdir[64]="/data/catalogs"; #define MAX_LTOK 80 static double ctg2ra(); static double ctg2dec(); static int ctgsize(); double dt2ep(); /* Julian Date to epoch (fractional year) */ static char newline = 10; /* CTGREAD -- Read ASCII stars in specified region */ int ctgread (catfile, refcat, distsort, cra, cdec, dra, ddec, drad, sysout, eqout, epout, mag1, mag2, sortmag, nsmax, starcat, tnum, tra, tdec, tpra, tpdec, tmag, tc, tobj, nlog) char *catfile; /* Name of reference star catalog file */ int refcat; /* Catalog code from wcctg.h */ 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) */ 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 nsmax; /* Maximum number of stars to be returned */ struct StarCat **starcat; /* Catalog data structure */ double *tnum; /* Array of ID 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 *tc; /* Array of fluxes (returned) */ char **tobj; /* Array of object names (returned) */ int nlog; { double ra1,ra2; /* Limiting right ascensions of region in degrees */ double dec1,dec2; /* Limiting declinations of region in degrees */ 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 */ char cstr[32]; struct Star *star; struct StarCat *sc; /* Catalog data structure */ char *objname; int nameobj; /* Save object name if 1, else do not */ int lname; int wrap; int imag; int jstar; int nstar; int magsort; double ra,dec,rapm,decpm; double mag; double num; int peak, i; int istar; int verbose; int isp; nstar = 0; /* Call the appropriate search program if not TDC ASCII catalog */ if (refcat != TXTCAT) { if (refcat == GSC || refcat == GSCACT) nstar = gscread (refcat,cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1,mag2,nsmax, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == GSC2) nstar = gsc2read (cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1,mag2,sortmag,nsmax, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == USAC || refcat == USA1 || refcat == USA2 || refcat == UAC || refcat == UA1 || refcat == UA2) nstar = uacread (catfile,distsort, cra,cdec,dra,ddec,drad,sysout,eqout,epout,mag1, mag2,sortmag,nsmax,tnum,tra,tdec,tmag,tc,nlog); else if (refcat == UJC || refcat == USNO) nstar = ujcread (catfile,cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1,mag2,nsmax, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == TMPSC) nstar = tmcread (cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1,mag2,sortmag,nsmax, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == ACT) nstar = actread (cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1, mag2,sortmag,nsmax, tnum,tra,tdec,tpra,tpdec,tmag,tc,nlog); else if (refcat == TYCHO2) nstar = ty2read (cra,cdec,dra,ddec,drad,distsort, sysout,eqout,epout,mag1,mag2,nsmax,sortmag, tnum,tra,tdec,tpra,tpdec,tmag,tc,nlog); else if (refcat == SAO) nstar = binread ("SAOra", distsort,cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == PPM) nstar = binread ("PPMra",distsort,cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == IRAS) nstar = binread ("IRAS", distsort, cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == TYCHO) nstar = binread ("tychora", distsort,cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == HIP) nstar = binread ("hipparcosra", distsort, cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == BSC) nstar = binread ("BSC5ra", distsort, cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == BINCAT) nstar = binread (catfile, distsort, cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,tobj,nlog); else if (refcat == TABCAT || refcat == WEBCAT) nstar = tabread (catfile, distsort,cra,cdec,dra,ddec,drad, sysout,eqout,epout,mag1,mag2,sortmag,nsmax,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,tobj,nlog); return (nstar); } star = NULL; sc = *starcat; if (nlog > 0) verbose = 1; else verbose = 0; wcscstr (cstr, sysout, eqout, epout); SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose); /* If RA range includes zero, split it in two */ wrap = 0; if (ra1 > ra2) wrap = 1; else wrap = 0; /* Search zones which include the poles cover 360 degrees in RA */ if (cdec - ddec < -90.0) { if (dec1 > dec2) dec2 = dec1; dec1 = -90.0; ra1 = 0.0; ra2 = 359.99999; wrap = 0; } if (cdec + ddec > 90.0) { if (dec2 < dec1) dec1 = dec2; dec2 = 90.0; ra1 = 0.0; ra2 = 359.99999; wrap = 0; } /* mag1 is always the smallest magnitude */ if (mag2 < mag1) { mag = mag2; mag2 = mag1; mag1 = mag; } /* Allocate catalog entry buffer */ star = (struct Star *) calloc (1, sizeof (struct Star)); star->num = 0.0; /* Logging interval */ tdist = (double *) malloc (nsmax * sizeof (double)); /* Open catalog file */ if (sc == NULL) { if ((sc = ctgopen (catfile, refcat)) == NULL) { fprintf (stderr,"CTGRNUM: Cannot read catalog %s\n", catfile); *starcat = sc; return (0); } } *starcat = sc; if (sc->nstars <= 0) { free (sc); if (star != NULL) free (star); sc = NULL; return (0); } jstar = 0; if (tobj == NULL || sc->ignore) nameobj = 0; else nameobj = 1; /* Loop through catalog */ for (istar = 1; istar <= sc->nstars; istar++) { if (ctgstar (istar, sc, star)) { fprintf (stderr,"\nCTGREAD: Cannot read %s star %d\n", sc->isfil, istar); break; } /* Set coordinate system for this star */ sysref = star->coorsys; eqref = star->equinox; epref = star->epoch; /* Extract selected fields */ num = star->num; ra = star->ra; dec = star->dec; rapm = star->rapm; decpm = star->decpm; if (sc->inform != 'X') { if (sc->mprop == 1) wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); else wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); } mag = star->xmag[magsort]; peak = 0; if (drad > 0 || distsort) { if (sc->inform == 'X') dist = sqrt ((ra-cra)*(ra-cra) + (dec-cdec)*(dec-cdec)); else dist = wcsdist (cra,cdec,ra,dec); } else dist = 0.0; /* If catalog is RA-sorted, stop reading if past highest RA */ if (sc->rasorted && !wrap && ra > ra2) break; /* Check magnitude and position limits */ if ((mag1 == mag2 || (mag >= mag1 && mag <= mag2)) && ((wrap && (ra >= ra1 || ra <= ra2)) || (!wrap && (ra >= ra1 && ra <= ra2))) && ((drad > 0.0 && dist < drad) || (drad == 0.0 && dec >= dec1 && dec <= dec2))) { /* Spectral type */ if (sc->sptype) isp = (1000 * (int) star->isp[0]) + (int)star->isp[1]; /* Save star position and magnitude in table */ if (nstar < nsmax) { tnum[nstar] = num; tra[nstar] = ra; tdec[nstar] = dec; for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][nstar] = star->xmag[imag]; } if (sc->sptype) tc[nstar] = isp; tdist[nstar] = dist; if (nameobj) { lname = strlen (star->objname) + 1; objname = (char *)calloc (lname, 1); strcpy (objname, star->objname); tobj[nstar] = objname; } if (dist > maxdist) { maxdist = dist; farstar = nstar; } if (mag > faintmag) { faintmag = mag; faintstar = nstar; } } /* If too many stars and distance sorting, replace furthest star */ else if (distsort) { if (dist < maxdist) { tnum[farstar] = num; tra[farstar] = ra; tdec[farstar] = dec; for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][farstar] = star->xmag[imag]; } if (sc->sptype) tc[farstar] = isp; tdist[farstar] = dist; if (nameobj) { free (tobj[farstar]); lname = strlen (star->objname) + 1; objname = (char *)calloc (lname, 1); strcpy (objname, star->objname); tobj[farstar] = objname; } /* Find new farthest star */ maxdist = 0.0; for (i = 0; i < nsmax; i++) { if (tdist[i] > maxdist) { maxdist = tdist[i]; farstar = i; } } } } /* If too many stars, replace faintest star */ else if (mag < faintmag) { tnum[faintstar] = num; tra[faintstar] = ra; tdec[faintstar] = dec; for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][faintstar] = star->xmag[imag]; } if (sc->sptype) tc[faintstar] = isp; tdist[faintstar] = dist; if (nameobj) { free (tobj[faintstar]); lname = strlen (star->objname) + 1; objname = (char *)calloc (lname, 1); strcpy (objname, star->objname); tobj[faintstar] = objname; } faintmag = 0.0; /* Find new faintest star */ for (i = 0; i < nsmax; i++) { if (tmag[magsort][i] > faintmag) { faintmag = tmag[magsort][i]; faintstar = i; } } } nstar++; jstar++; if (nlog == 1) fprintf (stderr,"CTGREAD: %11.6f: %9.5f %9.5f %s %5.2f %d \n", num,ra,dec,cstr,mag,peak); /* End of accepted star processing */ } /* Log operation */ if (nlog > 0 && istar%nlog == 0) fprintf (stderr,"CTGREAD: %5d / %5d / %5d sources catalog %s\r", jstar,istar,sc->nstars,catfile); /* End of star loop */ } /* Summarize search */ if (nlog > 0) { fprintf (stderr,"CTGREAD: Catalog %s : %d / %d / %d found\n", catfile,jstar,istar,sc->nstars); if (nstar > nsmax) fprintf (stderr,"CTGREAD: %d stars found; only %d returned\n", nstar,nsmax); } free ((char *)tdist); free (star); return (nstar); } /* CTGRNUM -- Read ASCII stars with specified numbers */ int ctgrnum (catfile,refcat, nnum,sysout,eqout,epout,match,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,tobj,nlog) char *catfile; /* Name of reference star catalog file */ int refcat; /* Catalog code from wcctg.h */ 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) */ int match; /* 1 to match star number exactly, else sequence num.*/ struct StarCat **starcat; /* Catalog data structure */ 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 *tc; /* Array of fluxes (returned) */ char **tobj; /* Array of object names (returned) */ int nlog; { int jnum; int nstar; double ra,dec; double rapm, decpm; int peak; int istar; int sysref; /* Catalog coordinate system */ double eqref; /* Catalog equinox */ double epref; /* Catalog epoch */ char cstr[32]; struct StarCat *sc; struct Star *star; char *objname; int imag; int lname; int starfound; int nameobj; nstar = 0; /* Call the appropriate search program if not TDC ASCII catalog */ if (refcat != TXTCAT) { if (refcat == GSC || refcat == GSCACT) nstar = gscrnum (refcat,nnum,sysout,eqout,epout, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == GSC2) nstar = gsc2rnum (nnum,sysout,eqout,epout, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == USAC || refcat == USA1 || refcat == USA2 || refcat == UAC || refcat == UA1 || refcat == UA2) nstar = uacrnum (catfile,nnum,sysout,eqout,epout, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == UJC || refcat == USNO) nstar = ujcrnum (catfile,nnum,sysout,eqout,epout, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == TMPSC) nstar = tmcrnum (nnum,sysout,eqout,epout, tnum,tra,tdec,tmag,tc,nlog); else if (refcat == SAO) nstar = binrnum ("SAO",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == PPM) nstar = binrnum ("PPM",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == IRAS) nstar = binrnum ("IRAS",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == TYCHO) nstar = binrnum ("tycho",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == HIP) nstar = binrnum ("hipparcos",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == BSC) nstar = binrnum ("BSC5",nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,NULL,nlog); else if (refcat == ACT) nstar = actrnum (nnum,sysout,eqout,epout, tnum,tra,tdec,tpra,tpdec,tmag,tc,nlog); else if (refcat == TYCHO2) nstar = ty2rnum (nnum,sysout,eqout,epout, tnum,tra,tdec,tpra,tpdec,tmag,tc,nlog); else if (refcat == TABCAT || refcat == WEBCAT) nstar = tabrnum (catfile,nnum,sysout,eqout,epout,starcat, tnum,tra,tdec,tpra,tpdec,tmag,tc,tobj,nlog); else if (refcat == BINCAT) nstar = binrnum (catfile,nnum,sysout,eqout,epout,match, tnum,tra,tdec,tpra,tpdec,tmag,tc,tobj,nlog); return (nstar); } sc = *starcat; if (sc == NULL) { if ((sc = ctgopen (catfile, refcat)) == NULL) { fprintf (stderr,"CTGRNUM: Cannot read catalog %s\n", catfile); *starcat = sc; return (0); } } *starcat = sc; sysref = sc->coorsys; eqref = sc->equinox; epref = sc->epoch; if (!sysout) sysout = sysref; if (!eqout) eqout = eqref; if (!epout) epout = epref; /* Allocate catalog entry buffer */ star = (struct Star *) calloc (1, sizeof (struct Star)); star->num = 0.0; if (tobj == NULL || sc->ignore) nameobj = 0; else nameobj = 1; /* Loop through star list */ for (jnum = 0; jnum < nnum; jnum++) { /* Loop through catalog to star */ starfound = 0; if (match && sc->stnum > 0) { for (istar = 1; istar <= sc->nstars; istar++) { if (ctgstar (istar, sc, star)) { fprintf (stderr,"CTGRNUM: Cannot read star %d\n", istar); break; } if (star->num == tnum[jnum]) { starfound = 1; break; } } } else { istar = (int) (tnum[jnum] + 0.5); if (ctgstar (istar, sc, star)) { fprintf (stderr,"CTGRNUM: Cannot read star %d\n", istar); continue; } starfound = 1; } /* If star has been found in catalog */ if (starfound) { /* Extract selected fields */ ra = star->ra; dec = star->dec; rapm = star->rapm; decpm = star->decpm; /* Set coordinate system for this star */ sysref = star->coorsys; eqref = star->equinox; epref = star->epoch; if (sc->inform != 'X') { if (sc->mprop == 1) wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); else wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); } peak = 0; /* Save star position and magnitude in table */ tnum[jnum] = star->num; tra[jnum] = ra; tdec[jnum] = dec; for (imag = 0; imag < sc->nmag; imag++) { if (tmag[imag] != NULL) tmag[imag][nstar] = star->xmag[imag]; } /* Spectral type */ if (sc->sptype) tc[jnum] = (1000 * (int) star->isp[0]) + (int)star->isp[1]; if (nameobj) { lname = strlen (star->objname) + 1; objname = (char *)calloc (lname, 1); strcpy (objname, star->objname); tobj[nstar] = objname; } nstar++; if (nlog == 1) fprintf (stderr,"CTGRNUM: %11.6f: %9.5f %9.5f %s %5.2f %d \n", star->num,ra,dec,cstr,star->xmag[0],peak); /* End of accepted star processing */ } /* Log operation */ if (nlog > 0 && jnum%nlog == 0) fprintf (stderr,"CTGRNUM: %5d / %5d / %5d sources catalog %s\r", nstar,jnum,sc->nstars,catfile); /* End of star loop */ } /* Summarize search */ if (nlog > 0) fprintf (stderr,"CTGRNUM: Catalog %s : %d / %d found\n", catfile,nstar,sc->nstars); free (star); return (nstar); } /* CTGOPEN -- Open ASCII catalog, returning number of entries */ struct StarCat * ctgopen (catfile, refcat) char *catfile; /* ASCII catalog file name */ int refcat; /* Catalog code from wcctg.h (TXTCAT,BINCAT,TABCAT) */ { struct StarCat *sc; struct Tokens tokens; FILE *fcat; char header[80]; char catpath[128]; /* Full pathname for catalog file */ char *catname; char *str; int nr, lfile, lhead, ldesc; char *catnew, *catdesc; char ctemp, *line, *linend, *cdot; char token[MAX_LTOK]; int ntok; int ltok; if (refcat != TXTCAT) { if (refcat == BINCAT) sc = binopen (catfile); else if (refcat == TABCAT) sc = tabcatopen (catfile, NULL, 0); else sc = NULL; return (sc); } /* Find length of ASCII catalog */ lfile = ctgsize (catfile); /* If catalog is not in current directory, look elsewhere */ if (lfile < 2) { /* Prepend directory name file not in working directory */ if ((str = getenv("WCS_CATDIR")) != NULL ) strcpy (catpath, str); else strcpy (catpath, catdir); strcat (catpath, "/"); strcat (catpath, catfile); lfile = ctgsize (catpath); if (lfile < 2) { fprintf (stderr,"CTGOPEN: ASCII catalog %s has no entries\n",catfile); return (NULL); } } else strcpy (catpath, catfile); /* Allocate catalog data structure */ sc = (struct StarCat *) calloc (1, sizeof (struct StarCat)); /* Open ASCII catalog */ if (!(fcat = fopen (catpath, "r"))) { fprintf (stderr,"CTGOPEN: ASCII catalog %s cannot be read\n",catfile); free (sc); return (NULL); } /* Allocate buffer to hold entire catalog and read it */ if ((sc->catbuff = malloc (lfile+2)) == NULL) { fprintf (stderr,"CTGOPEN: Cannot allocate memory for ASCII catalog %s\n", catfile); fclose (fcat); free (sc); return (NULL); } sc->catbuff[lfile] = (char) 0; sc->catbuff[lfile+1] = (char) 0; /* Separate filename from pathname and save in structure */ catname = strrchr (catfile,'/'); if (catname) catname = catname + 1; else catname = catfile; if (strlen (catname) < 24) strcpy (sc->isfil, catname); else strncpy (sc->isfil, catname, 23); /* Read entire catalog into memory at once */ nr = fread (sc->catbuff, 1, lfile, fcat); if (nr < lfile) { fprintf (stderr,"CTGOPEN: read only %d / %d bytes of file %s\n", nr, lfile, catfile); (void) fclose (fcat); free (sc); return (NULL); } /* Extract catalog information from first line */ sc->inform = 'H'; sc->coorsys = WCS_B1950; sc->epoch = 1950.0; sc->equinox = 1950.0; sc->nmag = 1; sc->mprop = 0; sc->rasorted = 0; sc->sptype = 0; sc->stnum = 1; sc->entepoch = 0; sc->entrv = 0; catdesc = strchr (sc->catbuff, newline) + 1; lhead = catdesc - sc->catbuff; if (lhead > 79) lhead = 79; if (sc->catbuff[0] == '#' && sc->catbuff[1] == ' ') { strncpy (header, sc->catbuff+2, lhead-2); header[lhead-2] = (char) 0; } else if (sc->catbuff[0] == '#') { strncpy (header, sc->catbuff+1, lhead-1); header[lhead-1] = (char) 0; } else { strncpy (header, sc->catbuff, lhead); header[lhead] = (char) 0; } /* Catalog positions are in B1950 (FK4) coordinates */ if (strsrch (header, "/b") || strsrch (header, "/B")) { sc->coorsys = WCS_B1950; sc->epoch = 1950.0; sc->equinox = 1950.0; } /* Catalog positions are in degrees */ if (strsrch (header, "/d") || strsrch (header, "/D")) sc->inform = 'D'; /* Catalog positions are in ecliptic coordinates */ if (strsrch (header, "/e") || strsrch (header, "/E")) { sc->coorsys = WCS_ECLIPTIC; sc->inform = 'D'; sc->epoch = 2000.0; sc->equinox = 2000.0; } /* Catalog positions are in hhmmss.s ddmmss.s */ if (strsrch (header, "/f") || strsrch (header, "/F")) sc->inform = 'F'; /* Catalog positions are galactic coordinates */ if (strsrch (header, "/g") || strsrch (header, "/G")) { sc->coorsys = WCS_GALACTIC; sc->inform = 'D'; sc->epoch = 2000.0; sc->equinox = 2000.0; } /* Catalog positions are in hh.mmsssss dd.mmssss format */ if (strsrch (header, "/h") || strsrch (header, "/H")) sc->inform = 'H'; /* Ignore information after position and magnitude */ if (strsrch (header, "/i") || strsrch (header, "/I")) sc->ignore = 1; else sc->ignore = 0; /* Catalog positions are J2000 (FK5) coordinates */ if (strsrch (header, "/j") || strsrch (header, "/J")) { sc->coorsys = WCS_J2000; sc->epoch = 2000.0; sc->equinox = 2000.0; } /* No magnitude */ if (strsrch (header, "/m") || strsrch (header, "/M")) sc->nmag = 0; else if (strsrch (header, "/2")) sc->nmag = 2; else if (strsrch (header, "/3")) sc->nmag = 3; else if (strsrch (header, "/4")) sc->nmag = 4; /* No number in first column, RA or object name first */ if (strsrch (header, "/n") || strsrch (header, "/N")) sc->stnum = 0; /* Object name instead of number in first column */ if (strsrch (header, "/o") || strsrch (header, "/O")) sc->stnum = -16; /* Proper motion */ if (strsrch (header, "/p") || strsrch (header, "/P")) sc->mprop = 1; if (strsrch (header, "/q") || strsrch (header, "/Q")) { sc->coorsys = 0; sc->epoch = 0.0; sc->equinox = 0.0; } /* RA-sorted catalog */ if (strsrch (header, "/r") || strsrch (header, "/R")) sc->rasorted = 1; /* Spectral type present */ if (strsrch (header, "/s") || strsrch (header, "/S")) sc->sptype = 1; /* Table format (hh mm ss dd mm ss) */ if (strsrch (header, "/t") || strsrch (header, "/T")) sc->inform = 'T'; /* Radial velocity is included after magnitude(s) */ if (strsrch (header, "/v") || strsrch (header, "/V")) { sc->mprop = 2; sc->entrv = 1; } /* X Y format (x.xxxxxx y.yyyyy) */ if (strsrch (header, "/x") || strsrch (header, "/X")) { sc->inform = 'X'; sc->coorsys = WCS_XY; } /* Epoch per catalog object as yyyy.mmdd */ if (strsrch (header, "/y") || strsrch (header, "/Y")) sc->nepoch = 1; else sc->nepoch = 0; /* Second line is description */ sc->catdata = strchr (catdesc, newline) + 1; ldesc = sc->catdata - catdesc; if (ldesc > 63) ldesc = 63; if (catdesc[0] == '#' && catdesc[1] == ' ') { strncpy (sc->isname, catdesc+2, ldesc-2); sc->isname[ldesc-2] = (char) 0; } else if (catdesc[0] == '#') { strncpy (sc->isname, catdesc+1, ldesc-1); sc->isname[ldesc-1] = (char) 0; } else { strncpy (sc->isname, catdesc, ldesc); sc->isname[ldesc] = (char) 0; } if (sc->entrv > 0 && sc->nmag < 2) { sc->nmag = sc->nmag + 1; if (sc->nmag == 2) strcpy (sc->keymag2, "velocity"); else if (sc->nmag == 3) strcpy (sc->keymag3, "velocity"); else if (sc->nmag == 4) strcpy (sc->keymag4, "velocity"); else strcpy (sc->keymag1, "velocity"); } /* Enumerate entries in ASCII catalog by counting newlines */ catnew = sc->catdata; sc->nstars = 0; while ((catnew = strchr (catnew, newline)) != NULL) { catnew = catnew + 1; if (*catnew != '#') sc->nstars = sc->nstars + 1; } sc->catline = sc->catdata; sc->catlast = sc->catdata + lfile; sc->istar = 1; /* Check number of decimal places in star number, if present */ if (sc->stnum == 1) { /* Temporarily terminate line with 0 */ line = sc->catline; linend = strchr (sc->catline, newline); if (linend == NULL) linend = sc->catlast; ctemp = *linend; *linend = (char) 0; /* Extract information from line of catalog */ ntok = setoken (&tokens, line, NULL); ltok = nextoken (&tokens, token, MAX_LTOK); sc->nnfld = ltok; if (ltok > 0) { sc->nndec = 0; if ((cdot = strchr (token,'.')) != NULL) { while (*(++cdot) != (char)0) sc->nndec++; } } else sc->nndec = 0; *linend = ctemp; } sc->refcat = TXTCAT; (void) fclose (fcat); return (sc); } /* CTGCLOSE -- Close ASCII catalog and free associated data structures */ void ctgclose (sc) struct StarCat *sc; { if (sc == NULL) return; else if (sc->refcat == BINCAT) binclose (sc); else if (sc->refcat == TABCAT) tabcatclose (sc); else if (sc->refcat == TXTCAT) { free (sc->catbuff); free (sc); } else free (sc); sc = NULL; return; } /* CTGSTAR -- Get ASCII catalog entry for one star; return 0 if successful */ int ctgstar (istar, sc, st) int istar; /* Star sequence number in ASCII catalog */ struct StarCat *sc; /* Star catalog data structure */ struct Star *st; /* Star data structure, updated on return */ { struct Tokens tokens; double ydate, dtemp; char *line; char *nextline; char token[80]; char ctemp, *linend; int ntok, itok; int ltok; int imag, nmag; /* Return error if requested number beyond catalog */ if (istar > sc->nstars) { fprintf (stderr, "CTGSTAR: %d is not in catalog\n",istar); return (-1); } /* If star is 0, read next star in catalog */ else if (istar < 1) { line = strchr (sc->catline, newline) + 1; while (*line == '#') line = strchr (line, newline) + 1; if (line == NULL) return (-1); else sc->catline = line; } /* If star is before current star, read from start of catalog */ else if (istar < sc->istar) { sc->istar = 1; sc->catline = sc->catdata; nextline = sc->catline; while (sc->istar < istar) { nextline = strchr (nextline, newline) + 1; while (*nextline == '#') nextline = strchr (nextline, newline) + 1; if (nextline == NULL) return (-1); sc->catline = nextline; sc->istar++; } } /* If star is after current star, read forward to it */ else if (istar > sc->istar) { nextline = sc->catline; while (sc->istar < istar) { nextline = strchr (nextline, newline) + 1; while (*nextline == '#') nextline = strchr (nextline, newline) + 1; if (nextline == NULL) return (-1); sc->catline = nextline; sc->istar++; } } /* Temporarily terminate line with 0 */ line = sc->catline; linend = strchr (sc->catline, newline); if (linend == NULL) linend = sc->catlast; ctemp = *linend; *linend = (char) 0; st->objname[0] = (char) 0; /* Extract information from line of catalog */ ntok = setoken (&tokens, line, NULL); /* Source number */ if (sc->stnum > 0) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 0) { st->num = atof (token); if (st->num - ((double) ((int) st->num)) < 0.000001) sc->stnum = 1; else sc->stnum = 2; } } else st->num = (double) sc->istar; /* Object name, if at start of line */ if (sc->stnum < 0) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 31) { strncpy (st->objname, token, 31); st->objname[31] = 0; } else if (ltok > 0) strcpy (st->objname, token); } /* Right ascension or longitude */ ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); /* Translate 3-token right ascension (hh mm ss.ss) */ if (sc->inform == 'T') { int hr, mn; double sec; hr = (int) (atof (token) + 0.5); ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); mn = (int) (atof (token) + 0.5); ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); sec = atof (token); st->ra = 15.0 * ((double) hr + ((double) mn / 60.0) + (sec / 3600.0)); } /* Translate single-token right ascension from hhmmss.ss */ else if (sc->inform == 'F') { dtemp = 0.0001 * atof (token); sprintf (token, "%.6f", dtemp); st->ra = ctg2ra (token); } /* Translate single-token right ascension as degrees */ else if (sc->inform == 'D') st->ra = atof (token); /* Translate single-token right ascension as X image coordinate */ else if (sc->inform == 'X') st->ra = atof (token); /* Translate single-token right ascension (hh:mm:ss.ss or hh.mmssss) */ else st->ra = ctg2ra (token); /* Declination or latitude */ ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); /* Translate 3-token declination (sdd mm ss.ss) */ if (sc->inform == 'T') { int deg, min; int decsgn = 0; double sec; deg = (int) (atof (token) + 0.5); if (strchr (token, '-') != NULL) decsgn = 1; ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); min = (int) (atof (token) + 0.5); ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); sec = atof (token); st->dec = (double) deg + ((double) min / 60.0) + (sec / 3600.0); if (deg > 0 && decsgn) st->dec = -st->dec; } /* Translate single-token declination as Y image coordinate */ else if (sc->inform == 'X') st->dec = atof (token); /* Translate single-token declination as degrees */ else if (sc->inform == 'D') st->dec = atof (token); /* Translate single-token declination from ddmmss.ss */ else if (sc->inform == 'F') { dtemp = 0.0001 * atof (token); sprintf (token, "%.6f", dtemp); st->dec = ctg2dec (token); } /* Translate single-token declination (dd:mm:ss.ss or dd.mmssss) */ else st->dec = ctg2dec (token); /* Equinox, if not set by header flag */ if (sc->coorsys == 0) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok < 1) return (-1); st->coorsys = wcscsys (token); st->equinox = wcsceq (token); st->epoch = sc->epoch; if (st->epoch == 0.0) st->epoch = st->equinox; } else { st->coorsys = sc->coorsys; st->equinox = sc->equinox; st->epoch = sc->epoch; } /* Magnitude, if present */ if (sc->nmag > 0 && sc->mprop == 2) nmag = sc->nmag - 1; else nmag = sc->nmag; if (nmag > 0) { for (imag = 0; imag < nmag; imag++) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 0) st->xmag[imag] = atof (token); else st->xmag[imag] = 0.0; } } /* Spectral type, if present */ if (sc->sptype > 0) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 0) { if (token[0] == ' ' && token[1] == ' ') { st->isp[0] = '_'; st->isp[1] = '_'; } else { st->isp[0] = token[0]; st->isp[1] = token[1]; } } } /* Radial velocity, if present */ if (sc->entrv > 0) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 0) st->radvel = atof (token); else st->radvel = 0.0; st->xmag[nmag] = st->radvel; } /* Proper motion, if present */ if (sc->mprop == 1) { ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 1) st->rapm = atof (token) / 3600.0; ltok = nextoken (&tokens, token, MAX_LTOK); if (ltok > 1) st->decpm = atof (token) / 3600.0; } /* Epoch, if present */ if (sc->nepoch) { ltok = nextoken (&tokens, token, MAX_LTOK); ydate = atof (token); st->epoch = dt2ep (ydate, 12.0); } /* Object name */ itok = tokens.itok; if (sc->stnum > 0 && itok < ntok && !sc->ignore) { itok = -(itok+1); ltok = getoken (&tokens, itok, token, MAX_LTOK); if (ltok > 31) { strncpy (st->objname, token, 31); st->objname[31] = 0; } else if (ltok > 0) strcpy (st->objname, token); } *linend = ctemp; return (0); } /* CTGSIZE -- return size of ASCII catalog file in bytes */ static int ctgsize (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); } /* ISACAT -- Return 1 if file is likely to be an ASCII catalog */ int isacat (catpath) char *catpath; { char buf[100]; char *errc; FILE *fcat; /* Open ASCII catalog */ if (!(fcat = fopen (catpath, "r"))) { return (0); } errc = fgets (buf, 100, fcat); fclose (fcat); if (isnum (buf)) return (0); else return (1); } /* CTG2RA -- Read the right ascension, ra, in sexagesimal hours from in[] */ static double ctg2ra (in) char *in; /* Character string */ { double ra; /* Right ascension in degrees (returned) */ ra = ctg2dec (in); ra = ra * 15.0; return (ra); } /* CTG2DEC -- Read the declination, dec, in sexagesimal degrees from in[] */ static double ctg2dec (in) char *in; /* Character string */ { double dec; /* Declination in degrees (returned) */ double deg, min, sec, sign; char *value, *c1; dec = 0.0; /* Translate value from ASCII colon-delimited string to binary */ if (!in[0]) return (dec); else value = in; /* Set sign */ if (!strchr (value,'-')) sign = 1.0; else { sign = -1.0; value = strchr (value,'-') + 1; } /* Translate value from ASCII colon-delimited string to binary */ if ((c1 = strchr (value,':')) != NULL) { *c1 = 0; deg = (double) atoi (value); *c1 = ':'; value = c1 + 1; if ((c1 = strchr (value,':')) != NULL) { *c1 = 0; min = (double) atoi (value); *c1 = ':'; value = c1 + 1; sec = atof (value); } else { sec = 0.0; if ((c1 = strchr (value,'.')) != NULL) min = atof (value); if (strlen (value) > 0) min = (double) atoi (value); } dec = sign * (deg + (min / 60.0) + (sec / 3600.0)); } /* Translate value from dd.mmssss */ else if ((c1 = strchr (value,'.')) != NULL) { double xnum, deg, min, sec; xnum = atof (value); deg = (double)((int) (xnum + 0.000000001)); xnum = (xnum - deg) * 100.0; min = (double)((int) (xnum + 0.000000001)); sec = (xnum - min) * 100.0; dec = sign * (deg + (min / 60.0) + (sec / 3600.0)); } /* Translate integer */ else dec = sign * (double) atoi (value); return (dec); } /* Oct 16 1998 New subroutines * Oct 20 1998 Clean up error messages * Oct 26 1998 Return object names in catread() and catrnum() * Oct 29 1998 Correctly assign numbers when too many stars are found * Oct 30 1998 Fix epoch and equinox for J2000 * Nov 9 1998 Drop out of star loop if rasorted catalog and past max RA * Dec 8 1998 Do not declare catsize() static * Dec 21 1998 Fix parsing so first character of line is not dropped * Jan 20 1999 Use strchr() instead of strsrch() for single char searches * Jan 29 1999 Default to star id number present * Feb 1 1999 Rewrite tokenizing subroutines for clarity * Feb 1 1999 Add match argument to catrnum() * Feb 2 1999 Add code to count decimal places in numbers * Feb 2 1999 Set sysout, eqout, and epout in catrnum() if not set * Feb 10 1999 Implement per star coordinate system * Feb 11 1999 Change starcat.insys to starcat.coorsys for consistency * Feb 17 1999 Fix per star coordinate system bugs * May 20 1999 Add option to read epoch of coordinates * 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 catread() * Sep 13 1999 Use these subroutines for general catalog access * 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 5 1999 Move token subroutines to catutil.c * Oct 15 1999 Fix calls to catclose(); eliminate dno in catstar() * Oct 22 1999 Fix declarations after lint * Oct 22 1999 Rename subroutines from cat* to ctg* to avoid system conflict * Nov 16 1999 Transfer dec degree sign if table format * * Jan 10 2000 Add second magnitude to tabread() and tabrnum() * Mar 10 2000 Add proper motions to ctgread(), ctgrnum(), tabread(), tabrnum() * Mar 15 2000 Add proper motions to binread(), binrnum(), actread(), actrnum() * Apr 3 2000 Implement fractional degrees (/d) for positions * Apr 3 2000 Add /i option to ignore stuff at end of line * Apr 3 2000 Ignore leading # on first two lines of ASCII catalog file * May 20 2000 Add Tycho 2 catalog support * Jun 23 2000 Ignore any line with # as first character, not just heading * Jun 26 2000 Implement /X catalog flag to avoid conversions * Jun 26 2000 Add coordinate system to SearchLim() arguments * Jun 26 2000 Ignore lines starting with # when counting stars in catalog * Jul 12 2000 Add star catalog data structure to ctgread() argument list * Jul 25 2000 Pass star catalog address of data structure address * Sep 20 2000 Implement multiple catalog magnitudes; return only first 2 * Sep 20 2000 Add isacat() subroutine to detect ASCII catalogs * Sep 25 2000 Add spectral type with flag /s * Oct 17 2000 Add missing arguments in generic binrnum() call * Oct 24 2000 Use ujcread() and ujcrnum() for USNO plate catalogs * Nov 21 2000 Add WEBCAT as tab catalog results returned from the Web * Nov 28 2000 Add starcat structure to *rnum() calls * Dec 18 2000 Include math.h for sqrt() * * Feb 14 2001 Search all around RA if either pole is included * Feb 16 2001 Update comments * May 22 2001 Add GSC-ACT catalog; pass refcat to gscread() and gscrnum() * May 23 2001 Add 2MASS Point Source Catalog * May 29 2001 Save length of star i.d. number in ctgopen() * Jun 18 2001 Add maximum length of returned string to getoken(), nextoken() * Jun 20 2001 Add GSC II to ctgread() and ctgrnum() * Aug 8 2001 Add /v option to return radial velocity * Aug 8 2001 Add /f option for hhmmss.s ddmmss.s coordinates * Aug 24 2001 Use STNUM < 0 instead of STNUM ==5 for object name not number * Sep 11 2001 Allow an arbitrary number of magnitudes * Sep 11 2001 Add sort magnitude as argument to *read() subroutines * Sep 19 2001 Drop fitshead.h; it is in wcs.h * Sep 27 2001 Fix bug which reset number of magnitudes to 1 if /m not last */