/*** File libwcs/gscread.c *** September 22, 2009 *** By Jessica Mink, jmink@cfa.harvard.edu *** Harvard-Smithsonian Center for Astrophysics *** Copyright (C) 1996-2009 *** 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 */ #include #include #include #include #include #include "fitsfile.h" #include "wcs.h" #include "wcscat.h" /* Pathname of northern hemisphere GSC CDROM or search engine URL */ static char cdn[64]="/data/astrocat/gsc1"; /* Uncomment following line to use ESO GSC server for GSC static char cdn[64]="http://archive.eso.org/skycat/servers/gsc-server"; */ /* Pathname of southern hemisphere GSC CDROM */ static char cds[64]="/data/astrocat/gsc2"; /* Pathname of northern hemisphere GSC-ACT CDROM or search engine URL */ static char cdna[64]="/data/astrocat/gscact1"; /* Pathname of southern hemisphere GSC-ACT CDROM */ static char cdsa[64]="/data/astrocat/gscact2"; static void gscpath(); static int gscreg(); static char *table = NULL; /* FITS table buffer */ static int ltab = 0; /* Length of FITS table buffer */ static double *gdist = NULL; /* Array of distances to stars */ static int ndist = 0; void gscfree() { free ((void *)table); ltab = 0; free ((void *)gdist); ndist = 0; return; } static int classd = -1; /* Desired object class (-2=all objects, -1=all, 0=stars, 3=nonstars) */ void setgsclass (class) int class; { classd = class; return; } /* GSCREAD -- Read HST Guide Star Catalog stars from CDROM */ int gscread (refcat,cra,cdec,dra,ddec,drad,dradi,distsort,sysout,eqout,epout, mag1,mag2,nstarmax,gnum,gra,gdec,gmag,gtype,nlog) int refcat; /* Catalog code (GSC or GSCACT) */ 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 separation in degrees for annulus (ignore if 0) */ int distsort; /* 1 to sort stars by distance from center */ 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 nstarmax; /* Maximum number of stars to be returned */ double *gnum; /* Array of Guide Star numbers (returned) */ double *gra; /* Array of right ascensions (returned) */ double *gdec; /* Array of declinations (returned) */ double **gmag; /* Array of magnitudes (returned) */ int *gtype; /* Array of object types (returned) */ int nlog; /* 1 for diagnostics */ { 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 */ int magsort=0; int nreg; /* Number of input FITS tables files */ double xnum; /* Guide Star number */ int rlist[100]; /* List of input FITS tables files */ char inpath[64]; /* Pathname for input FITS table file */ char entry[100]; /* Buffer for FITS table row */ int class; /* Object class (0>star, 3>other) */ int sysref=WCS_J2000; /* Catalog coordinate system */ double eqref=2000.0; /* Catalog equinox */ double epref=2000.0; /* Catalog epoch */ struct Keyword kw[8]; /* Keyword structure */ struct Keyword *kwn; int verbose; int pass; int rnum, num0, num, itot,ireg; int ik,nk,itable,ntable,jstar; int nbline,npos,nbhead; int nbr,nrmax,nstar,i; int ift; int wrap; double ra,rasum,dec,decsum,perr,perr2,perrsum,msum; double mag,merr,merr2,merrsum; int band0 = 0; int band = 0; int class0 = 0; double merr0 = 0.0; double mag0 = 0.0; double perr0 = 0.0; double ra0 = 0.0; double dec0 = 0.0; double rra1, rra2, rdec1, rdec2; double rdist, ddist; char *str; char *url; char *title; char cstr[32], numstr[32], rastr[32], decstr[32], catid[16]; itot = 0; if (nlog == 1) verbose = 1; else verbose = 0; verbose = nlog; /* If root pathname is a URL, search and return */ url = cdn; if (refcat == GSC) { url = cdn; if ((str = getenv("GSC_NORTH")) == NULL) str = getenv ("GSC_PATH"); if (str != NULL) url = str; } if (refcat == GSCACT) { url = cdna; if ((str = getenv("GSCACT_NORTH")) == NULL) str = getenv ("GSCACT_PATH"); if (str != NULL) url = str; } if (!strncmp (url, "http:",5)) return (webread (url,"gsc",distsort,cra,cdec,dra,ddec,drad,dradi, sysout,eqout,epout,mag1,mag2,magsort,nstarmax, gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog)); /* Allocate FITS table buffer which is saved between calls */ if (ltab < 1) { ltab = 10000; table = (char *)calloc (ltab, sizeof (char)); if (table == NULL) { fprintf (stderr, "GSCREAD: cannot allocate FITS table buffer\n"); return (0); } } for (i = 0; i < 100; i++) entry[i] = 0; /* Set path to Guide Star Catalog */ if (refcat == GSCACT) { if ((str = getenv("GSCACT_NORTH")) != NULL ) strcpy (cdna,str); if ((str = getenv("GSCACT_SOUTH")) != NULL ) strcpy (cdsa,str); } else { if ((str = getenv("GSC_NORTH")) != NULL ) strcpy (cdn,str); if ((str = getenv("GSC_SOUTH")) != NULL ) strcpy (cds,str); } wcscstr (cstr, sysout, eqout, epout); SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose); /* Make mag1 always the smallest magnitude */ if (mag2 < mag1) { mag = mag2; mag2 = mag1; mag1 = mag; } /* Find Guide Star Catalog regions in which to search */ nrmax = 100; rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,0.0, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); nreg = gscreg (refcat,rra1,rra2,rdec1,rdec2,table,nrmax,rlist,verbose); if (nreg <= 0) { fprintf (stderr,"GSCREAD: no Guide Star regions found\n"); return (0); } /* Allocate array for distances from search center */ if (nstarmax > ndist) { if (ndist > 0) free ((void *)gdist); gdist = (double *) malloc (nstarmax * sizeof (double)); if (gdist == NULL) { fprintf (stderr,"GSCREAD: cannot allocate separation array\n"); return (0); } ndist = nstarmax; } /* Set keyword list */ nk = 8; strcpy (kw[0].kname,"GSC_ID"); strcpy (kw[1].kname,"RA_DEG"); strcpy (kw[2].kname,"DEC_DEG"); strcpy (kw[3].kname,"POS_ERR"); strcpy (kw[4].kname,"MAG"); strcpy (kw[5].kname,"MAG_ERR"); strcpy (kw[6].kname,"MAG_BAND"); strcpy (kw[7].kname,"CLASS"); for (ik = 0; ik < nk; ik++) { kw[ik].lname = (int) strlen (kw[ik].kname); kw[ik].kn = 0; kw[ik].kf = 0; kw[ik].kl = 0; } nstar = 0; /* Write header if printing star entries as found */ if (nstarmax < 1) { char *revmessage; revmessage = getrevmsg(); title = CatName (refcat, NULL); printf ("catalog %s\n", title); free (title); ra2str (rastr, 31, cra, 3); printf ("ra %s\n", rastr); dec2str (decstr, 31, cdec, 2); printf ("dec %s\n", decstr); if (drad != 0.0) { printf ("radmin %.1f\n", drad*60.0); if (dradi > 0) printf ("radimin %.1f\n", dradi*60.0); } else { printf ("dramin %.1f\n", dra*60.0* cosdeg (cdec)); printf ("ddecmin %.1f\n", ddec*60.0); } printf ("radecsys %s\n", cstr); printf ("equinox %.3f\n", eqout); printf ("epoch %.3f\n", epout); printf ("program scat %s\n", revmessage); CatID (catid, refcat); printf ("%s ra dec ", catid); printf ("magv class band n arcmin\n"); printf ("--------- ------------ ------------ "); printf ("----- ----- ---- - ------\n"); } /* Loop through region list */ for (ireg = 0; ireg < nreg; ireg++) { gscpath (refcat, rlist[ireg], inpath); /* Read size and keyword info from FITS table header */ kwn = kw; ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead); rnum = rlist[ireg]; num0 = 0; rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; num = 0; fitsrtlset(); jstar = 0; class = 0; /* Loop through FITS table for this region */ for (itable = 0; itable <= ntable; itable++) { if (itable < ntable) { nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry); if (nbr < nbline) { fprintf (stderr,"GSCREAD: %d / %d bytes read, line %d / %d, region %d\n", nbr,nbline,itable,ntable,rnum); break; } /* Extract selected fields */ /* Star number within region */ num0 = ftgeti4 (entry, &kw[0]); /* Right ascension in degrees */ ra0 = ftgetr8 (entry, &kw[1]); /* Declination in degrees */ dec0 = ftgetr8 (entry, &kw[2]); /* Position error */ perr0 = ftgetr8 (entry, &kw[3]); /* Magnitude */ mag0 = ftgetr8 (entry, &kw[4]); /* Magnitude error */ merr0 = ftgetr8 (entry, &kw[5]); /* Bandpass code */ band0 = ftgeti4 (entry, &kw[6]); /* Object class code */ class0 = ftgeti4 (entry, &kw[7]); } else num0 = 0; /* Compute mean position and magnitude for object */ if (itable > 0 && npos > 0 && ((classd < -1 && band != band0) || (classd < -1 && class != class0) || num != num0)) { pass = 1; if (perrsum == 0.0 || merrsum == 0) pass = 0; else { ra = rasum / perrsum; dec = decsum / perrsum; mag = msum / merrsum; } if (pass > 0 && classd > -1 && class != classd) pass = 0; /* Check magnitude and position limits */ if (pass > 0 && mag1 != mag2 && (mag < mag1 || mag > mag2)) pass = 0; if (pass) { wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); /* Compute distance from search center */ if (drad > 0 || distsort) dist = wcsdist (cra,cdec,ra,dec); else dist = 0.0; /* Check position limits */ if (drad > 0) { if (dist > drad) pass = 0; if (dradi > 0.0 && dist < dradi) pass = 0; } else { ddist = wcsdist (cra,cdec,cra,dec); if (ddist > ddec) pass = 0; rdist = wcsdist (cra,dec,ra,dec); if (rdist > dra) pass = 0; } } if (pass) { class = class + (band * 100) + (npos * 10000); xnum = (double)rnum + (0.0001 * (double) num); /* Write star position and magnitudes to stdout */ if (nstarmax < 1) { CatNum (refcat, -9, 0, xnum, numstr); ra2str (rastr, 31, ra, 3); dec2str (decstr, 31, dec, 2); dist = wcsdist (cra,cdec,ra,dec) * 60.0; printf ("%s %s %s", numstr,rastr,decstr); printf (" %.2f %d %d %d %.2f\n", mag, class-(band*100)-(npos*10000), band, npos, dist); } /* Save star position in table */ else if (nstar < nstarmax) { gnum[nstar] = xnum; gra[nstar] = ra; gdec[nstar] = dec; gmag[0][nstar] = mag; gtype[nstar] = class; gdist[nstar] = dist; 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) { gnum[farstar] = xnum; gra[farstar] = ra; gdec[farstar] = dec; gmag[0][farstar] = mag; gtype[farstar] = class; gdist[farstar] = dist; maxdist = 0.0; /* Find new farthest star */ for (i = 0; i < nstarmax; i++) { if (gdist[i] > maxdist) { maxdist = gdist[i]; farstar = i; } } } } /* If too many stars, replace faintest star */ else if (mag < faintmag) { gnum[faintstar] = xnum; gra[faintstar] = ra; gdec[faintstar] = dec; gmag[0][faintstar] = mag; gtype[faintstar] = class; gdist[faintstar] = dist; faintmag = 0.0; /* Find new faintest star */ for (i = 0; i < nstarmax; i++) { if (gmag[0][i] > faintmag) { faintmag = gmag[0][i]; faintstar = i; } } } nstar++; jstar++; if (nlog == 1) fprintf (stderr,"GSCREAD: %04d.%04d: %9.5f %9.5f %s %5.2f %d %d\n", rnum,num,ra,dec,cstr,mag,class,npos); } /* Reset star position for averaging */ rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; } /* Add information from current line to current object */ /* Check object class */ if ((classd > -1 && class0 == classd) || classd < -2 || (classd < 0 && class0 != 5)) { perr = perr0; perr2 = perr * perr; if (perr2 <= 0.0) perr2 = 0.01; rasum = rasum + (ra0 / perr2); decsum = decsum + (dec0 / perr2); perrsum = perrsum + (1.0 / perr2); if (merr0 <= 0.0) merr0 = 0.01; merr = merr0; merr2 = merr * merr; msum = msum + (mag0 / merr2); merrsum = merrsum + (1.0 / merr2); num = num0; class = class0; band = band0; npos++; } /* Log operation */ if (nlog > 0 && itable%nlog == 0) fprintf (stderr,"GSCREAD: %4d / %4d: %5d / %5d / %5d sources, region %4d.%04d\r", ireg,nreg,jstar,itable,ntable,rlist[ireg],num0); /* End of region */ } /* Close region input file */ (void) close (ift); itot = itot + itable; if (nlog > 0) fprintf (stderr,"GSCREAD: %4d / %4d: %5d / %5d / %5d sources from region %4d \n", ireg+1,nreg,jstar,itable,ntable,rlist[ireg]); } /* Close output file and summarize transfer */ if (nlog > 0) { if (nreg > 1) fprintf (stderr,"GSCREAD: %d regions: %d / %d found\n",nreg,nstar,itot); else fprintf (stderr,"GSCREAD: 1 region: %d / %d found\n",nstar,itable); if (nstar > nstarmax ) fprintf (stderr,"GSCREAD: %d stars found; only %d returned\n", nstar,nstarmax); } return (nstar); } /* GSCRNUM -- Read HST Guide Star Catalog stars from CDROM */ int gscrnum (refcat, nstars, sysout, eqout, epout, gnum,gra,gdec,gmag,gtype,nlog) int refcat; /* Catalog code (GSC or GSCACT) */ int nstars; /* Number of stars to find */ int sysout; /* Search coordinate system */ double eqout; /* Search coordinate equinox */ double epout; /* Proper motion epoch (0.0 for no proper motion) */ double *gnum; /* Array of Guide Star numbers (returned) */ double *gra; /* Array of right ascensions (returned) */ double *gdec; /* Array of declinations (returned) */ double **gmag; /* Array of magnitudes (returned) */ int *gtype; /* Array of object types (returned) */ int nlog; /* 1 for diagnostics */ { char *table; /* FITS table */ char inpath[64]; /* Pathname for input FITS table file */ char entry[100]; /* Buffer for FITS table row */ int class; /* Object class (0>star, 3>other) */ int sysref=WCS_J2000; /* Catalog coordinate system */ double eqref=2000.0; /* Catalog equinox */ struct Keyword kw[8]; /* Keyword structure */ struct Keyword *kwn; int rnum, num0, num, itot; int ik,nk,itable,ntable,jstar; int nbline,npos,nbhead; int nbr,nstar,i, snum; int ift; double dnum; double ra,rasum,dec,decsum,perr,perr2,perrsum,msum; double mag,merr,merr2,merrsum; int band0 = 0; int band = 0; int class0 = 0; double merr0 = 0.0; double mag0 = 0.0; double perr0 = 0.0; double ra0 = 0.0; double dec0 = 0.0; char *str; char *url; itot = 0; /* If root pathname is a URL, search and return */ url = cdn; if (refcat == GSC) { url = cdn; if ((str = getenv("GSC_NORTH")) == NULL) str = getenv ("GSC_PATH"); if (str != NULL) url = str; if (!strncmp (url, "http:",5)) return (webrnum (str,"gsc",nstars,sysout,eqout,epout,1, gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog)); } if (refcat == GSCACT) { url = cdna; if ((str = getenv("GSCACT_NORTH")) == NULL) str = getenv ("GSCACT_PATH"); if (str != NULL) url = str; } if (!strncmp (url, "http:",5)) return (webrnum (str,"gsc",nstars,sysout,eqout,epout,1, gnum,gra,gdec,NULL,NULL,gmag,gtype,nlog)); if (ltab < 1) { ltab = 10000; table = (char *)calloc (ltab, sizeof (char)); } for (i = 0; i < 100; i++) entry[i] = 0; /* Set path to Guide Star Catalog */ if ((str = getenv("GSC_NORTH")) != NULL ) strcpy (cdn,str); if ((str = getenv("GSC_SOUTH")) != NULL ) strcpy (cds,str); /* Set keyword list */ nk = 8; strcpy (kw[0].kname,"GSC_ID"); strcpy (kw[1].kname,"RA_DEG"); strcpy (kw[2].kname,"DEC_DEG"); strcpy (kw[3].kname,"POS_ERR"); strcpy (kw[4].kname,"MAG"); strcpy (kw[5].kname,"MAG_ERR"); strcpy (kw[6].kname,"MAG_BAND"); strcpy (kw[7].kname,"CLASS"); for (ik = 0; ik < nk; ik++) { kw[ik].lname = (int) strlen (kw[ik].kname); kw[ik].kn = 0; kw[ik].kf = 0; kw[ik].kl = 0; } nstar = 0; /* Loop through star list */ for (jstar = 0; jstar < nstars; jstar++) { rnum = (int) gnum[jstar]; gscpath (refcat, rnum, inpath); dnum = (gnum[jstar] - (double)rnum) * 10000.0; snum = (int) (dnum + 0.5); /* Read size and keyword info from FITS table header */ kwn = kw; ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead); if (ift < 0) { fprintf (stderr,"GSCRNUM: File %s not found\n",inpath); return (0); } num0 = 0; rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; num = 0; fitsrtlset(); class = 0; /* Loop through FITS table for this region */ for (itable = 0; itable <= ntable; itable++) { if (itable < ntable) { nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry); if (nbr < nbline) { fprintf (stderr,"GSCRNUM: %d / %d bytes read, line %d / %d, region %d\n", nbr,nbline,itable,ntable,rnum); break; } /* Extract selected fields */ /* Star number within region */ num0 = ftgeti4 (entry, &kw[0]); if (num0 < snum) continue; else if (num == 0) num = num0; /* Right ascension in degrees */ ra0 = ftgetr8 (entry, &kw[1]); /* Declination in degrees */ dec0 = ftgetr8 (entry, &kw[2]); /* Position error */ perr0 = ftgetr8 (entry, &kw[3]); /* Magnitude */ mag0 = ftgetr8 (entry, &kw[4]); /* Magnitude error */ merr0 = ftgetr8 (entry, &kw[5]); /* Bandpass code */ band0 = ftgeti4 (entry, &kw[6]); /* Object class code */ class0 = ftgeti4 (entry, &kw[7]); } else num0 = 0; /* Compute mean position and magnitude for object */ if (num != num0 && itable > 0 && npos > 0) { ra = rasum / perrsum; dec = decsum / perrsum; wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); mag = msum / merrsum; class = class + (band * 100) + (npos * 10000); /* Save star position in table */ gra[nstar] = ra; gdec[nstar] = dec; gmag[0][nstar] = mag; gtype[nstar] = class; nstar++; if (nlog == 1) fprintf (stderr,"GSCRNUM: %04d.%04d: %9.5f %9.5f %5.2f %d %d\n", rnum,num,ra,dec,mag,class,npos); /* Reset star position for averaging */ rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; break; } /* Add information from current line to current object */ perr = perr0; perr2 = perr * perr; if (perr2 <= 0.0) perr2 = 0.01; rasum = rasum + (ra0 / perr2); decsum = decsum + (dec0 / perr2); perrsum = perrsum + (1.0 / perr2); if (merr0 <= 0.0) merr0 = 0.01; merr = merr0; merr2 = merr * merr; msum = msum + (mag0 / merr2); merrsum = merrsum + (1.0 / merr2); num = num0; class = class0; band = band0; npos = npos + 1; /* Log operation */ if (nlog > 0 && itable%nlog == 0) fprintf (stderr,"GSCRNUM: %4d / %4d: %5d / %5d sources, region %4d.%04d\r", jstar+1,nstars,itable,ntable,rnum,snum); /* End of region */ } /* Close region input file */ (void) close (ift); itot = itot + itable; if (nlog > 0) fprintf (stderr,"GSCRNUM: %4d / %4d: %5d / %5d sources, region %4d.%04d\n", jstar+1,nstars,itable,ntable,rnum,snum); } return (nstars); } /* GSCBIN -- Fill FITS WCS image with HST Guide Star Catalog objects */ int gscbin (refcat, wcs, header, image, mag1, mag2, magscale, nlog) int refcat; /* Catalog code (GSC or GSCACT) */ 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) */ double magscale; /* Scaling factor for magnitude to pixel flux * (number of catalog objects per bin if 0) */ int nlog; /* 1 for diagnostics */ { 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 = wcs->syswcs; /* Image coordinate system */ double eqout = wcs->equinox; /* Image coordinate equinox */ double epout = wcs->epoch; /* Image epoch */ double ra1,ra2; /* Limiting right ascensions of image in degrees */ double dec1,dec2; /* Limiting declinations of image in degrees */ int nreg; /* Number of input FITS tables files */ double xnum; /* Guide Star number */ int rlist[100]; /* List of input FITS tables files */ char inpath[64]; /* Pathname for input FITS table file */ char entry[100]; /* Buffer for FITS table row */ int class, class0; /* Object class (0>star, 3>other) */ int sysref=WCS_J2000; /* Catalog coordinate system */ double eqref=2000.0; /* Catalog equinox */ struct Keyword kw[8]; /* Keyword structure */ struct Keyword *kwn; int verbose; int pass; int wrap; int rnum, num0, num, itot,ireg; int ik,nk,itable,ntable,jstar; int nbline,npos,nbhead; int nbr,nrmax,nstar,i; int ift, band0, band; double ra,ra0,rasum,dec,dec0,decsum,perr,perr0,perr2,perrsum,msum; double mag,mag0,merr,merr0,merr2,merrsum; double rra1, rra2, rdec1, rdec2; double rdist, ddist; char *str; char cstr[32]; int bitpix, w, h; /* Image bits/pixel and pixel width and height */ double logt = log(10.0); double xpix, ypix, flux; int ix, iy; int offscl; itot = 0; if (nlog == 1) verbose = 1; else verbose = 0; verbose = nlog; /* Set image parameters */ bitpix = 0; (void)hgeti4 (header, "BITPIX", &bitpix); w = 0; (void)hgeti4 (header, "NAXIS1", &w); h = 0; (void)hgeti4 (header, "NAXIS2", &h); if (refcat == GSCACT) { if ((str = getenv("GSCACT_NORTH")) == NULL) str = getenv ("GSCACT_PATH"); } /* Allocate FITS table buffer which is saved between calls */ if (ltab < 1) { ltab = 10000; table = (char *)calloc (ltab, sizeof (char)); if (table == NULL) { fprintf (stderr, "GSCBIN: cannot allocate FITS table buffer\n"); return (0); } } for (i = 0; i < 100; i++) entry[i] = 0; /* Set path to Guide Star Catalog */ if (refcat == GSCACT) { if ((str = getenv("GSCACT_NORTH")) != NULL ) strcpy (cdna,str); if ((str = getenv("GSCACT_SOUTH")) != NULL ) strcpy (cdsa,str); } else { if ((str = getenv("GSC_NORTH")) != NULL ) strcpy (cdn,str); if ((str = getenv("GSC_SOUTH")) != NULL ) strcpy (cds,str); } wcscstr (cstr, sysout, eqout, epout); wcssize (wcs, &cra, &cdec, &dra, &ddec); SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose); /* Make mag1 always the smallest magnitude */ if (mag2 < mag1) { mag = mag2; mag2 = mag1; mag1 = mag; } /* Find Guide Star Catalog regions in which to search */ nrmax = 100; rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epout,0.0, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); nreg = gscreg (refcat,rra1,rra2,rdec1,rdec2,table,nrmax,rlist,verbose); if (nreg <= 0) { fprintf (stderr,"GSCBIN: no Guide Star regions found\n"); return (0); } /* Set keyword list */ nk = 8; strcpy (kw[0].kname,"GSC_ID"); strcpy (kw[1].kname,"RA_DEG"); strcpy (kw[2].kname,"DEC_DEG"); strcpy (kw[3].kname,"POS_ERR"); strcpy (kw[4].kname,"MAG"); strcpy (kw[5].kname,"MAG_ERR"); strcpy (kw[6].kname,"MAG_BAND"); strcpy (kw[7].kname,"CLASS"); for (ik = 0; ik < nk; ik++) { kw[ik].lname = (int) strlen (kw[ik].kname); kw[ik].kn = 0; kw[ik].kf = 0; kw[ik].kl = 0; } nstar = 0; /* Loop through region list */ for (ireg = 0; ireg < nreg; ireg++) { gscpath (refcat, rlist[ireg], inpath); /* Read size and keyword info from FITS table header */ kwn = kw; ift = fitsrtopen (inpath,&nk,&kwn,&ntable,&nbline,&nbhead); rnum = rlist[ireg]; num0 = 0; rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; num = 0; fitsrtlset(); jstar = 0; class = 0; /* Loop through FITS table for this region */ for (itable = 0; itable <= ntable; itable++) { if (itable < ntable) { nbr = fitsrtline (ift,nbhead,ltab,table,itable,nbline,entry); if (nbr < nbline) { fprintf (stderr,"GSCBIN: %d / %d bytes read, line %d / %d, region %d\n", nbr,nbline,itable,ntable,rnum); break; } /* Extract selected fields */ /* Star number within region */ num0 = ftgeti4 (entry, &kw[0]); /* Right ascension in degrees */ ra0 = ftgetr8 (entry, &kw[1]); /* Declination in degrees */ dec0 = ftgetr8 (entry, &kw[2]); /* Position error */ perr0 = ftgetr8 (entry, &kw[3]); /* Magnitude */ mag0 = ftgetr8 (entry, &kw[4]); /* Magnitude error */ merr0 = ftgetr8 (entry, &kw[5]); /* Bandpass code */ band0 = ftgeti4 (entry, &kw[6]); /* Object class code */ class0 = ftgeti4 (entry, &kw[7]); } else num0 = 0; /* Compute mean position and magnitude for object */ if (itable > 0 && npos > 0 && ((classd < -1 && band != band0) || (classd < -1 && class != class0) || num != num0)) { pass = 1; if (perrsum == 0.0 || merrsum == 0) pass = 0; else { ra = rasum / perrsum; dec = decsum / perrsum; mag = msum / merrsum; } if (pass > 0 && classd > -1 && class != classd) pass = 0; /* Check magnitude and position limits */ if (pass > 0 && mag1 != mag2 && (mag < mag1 || mag > mag2)) pass = 0; if (pass) { wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout); /* Check position limits */ 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) { class = class + (band * 100) + (npos * 10000); xnum = (double)rnum + (0.0001 * (double) num); wcs2pix (wcs, ra, dec, &xpix, &ypix, &offscl); if (!offscl) { if (magscale > 0.0) flux = magscale * exp (logt * (-mag / 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,"GSCBIN: %04d.%04d: %9.5f %9.5f %s %d %d", rnum, num, ra, dec, cstr, class, npos); 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); } } /* Reset star position for averaging */ rasum = 0.0; decsum = 0.0; msum = 0.0; perrsum = 0.0; merrsum = 0.0; npos = 0; } /* Add information from current line to current object */ /* Check object class */ if ((classd > -1 && class0 == classd) || classd < -2 || (classd < 0 && class0 != 5)) { perr = perr0; perr2 = perr * perr; if (perr2 <= 0.0) perr2 = 0.01; rasum = rasum + (ra0 / perr2); decsum = decsum + (dec0 / perr2); perrsum = perrsum + (1.0 / perr2); if (merr0 <= 0.0) merr0 = 0.01; merr = merr0; merr2 = merr * merr; msum = msum + (mag0 / merr2); merrsum = merrsum + (1.0 / merr2); num = num0; class = class0; band = band0; npos++; } /* Log operation */ if (nlog > 0 && itable%nlog == 0) fprintf (stderr,"GSCBIN: %4d / %4d: %5d / %5d / %5d sources, region %4d.%04d\r", ireg,nreg,jstar,itable,ntable,rlist[ireg],num0); /* End of region */ } /* Close region input file */ (void) close (ift); itot = itot + itable; if (nlog > 0) fprintf (stderr,"GSCBIN: %4d / %4d: %5d / %5d / %5d sources from region %4d \n", ireg+1,nreg,jstar,itable,ntable,rlist[ireg]); } /* Close output file and summarize transfer */ if (nlog > 0) { if (nreg > 1) fprintf (stderr,"GSCBIN: %d regions: %d / %d found\n",nreg,nstar,itot); else fprintf (stderr,"GSCBIN: 1 region: %d / %d found\n",nstar,itable); } return (nstar); } /* First region in each declination zone */ int zreg1[24]={1,594,1178,1729,2259,2781,3246,3652,4014,4294, 4492,4615, 4663,5260,5838,6412,6989,7523,8022,8464,8840,9134,9346,9490}; /* Last region in each declination zone */ int zreg2[24]={593,1177,1728,2258,2780,3245,3651,4013,4293,4491,4614,4662, 5259,5837,6411,6988,7522,8021,8463,8839,9133,9345,9489,9537}; /* Directory for each declination zone */ char zdir[24][8]={"n0000","n0730","n1500","n2230","n3000","n3730","n4500", "n5230","n6000","n6730","n7500","n8230","s0000","s0730", "s1500","s2230","s3000","s3730","s4500","s5230","s6000", "s6730","s7500","s8230"}; static struct Keyword rkw[15]; static int nrkw = 13; /* GSCREG -- search the HST Guide Star Catalog index table for fields * in the specified range of coordinates and magnitudes. Build a * list containing the pathnames of the files on the cd-rom. */ static int gscreg (refcat, ra1, ra2, dec1, dec2, table, nrmax, rgns, verbose) int refcat; /* Catalog code (GSC or GSCACT) */ double ra1, ra2; /* Right ascension limits in degrees */ double dec1, dec2; /* Declination limits in degrees */ char *table; /* Table data buffer */ int nrmax; /* Maximum number of regions to find */ int *rgns; /* Region numbers (returned)*/ int verbose; /* 1 for diagnostics */ { int nrgn; /* Number of regions found (returned) */ char tabpath[64]; /* Pathname for regions table */ int nrows; /* Number of entries in table */ int nchar; /* Number of characters per line in table */ int nwrap; /* 1 if 0h included in RA span*/ int iwrap; struct Keyword *kwn; double gscra(), gscdec(); int gsczone(); char fitsline[120]; int irow,iz1,iz2,ir1,ir2,jr1,jr2,i; int nsrch,nsrch1,nbhead,nbr; int ift; double ralow, rahi; double declow, dechi, decmin, decmax; int regnum; /* Set up keyword list for table entries to extract */ strcpy (rkw[0].kname,"REG_NO"); strcpy (rkw[1].kname,"RA_H_LOW"); strcpy (rkw[2].kname,"RA_M_LOW"); strcpy (rkw[3].kname,"RA_S_LOW"); strcpy (rkw[4].kname,"RA_H_HI"); strcpy (rkw[5].kname,"RA_M_HI"); strcpy (rkw[6].kname,"RA_S_HI"); strcpy (rkw[7].kname,"DECSI_LO"); strcpy (rkw[8].kname,"DEC_D_LO"); strcpy (rkw[9].kname,"DEC_M_LO"); strcpy (rkw[10].kname,"DECSI_HI"); strcpy (rkw[11].kname,"DEC_D_HI"); strcpy (rkw[12].kname,"DEC_M_HI"); rkw[13].kname[0] = 0; /* Add lengths of keywords to keyword structures */ for (i = 0; i < 13; i++) rkw[i].lname = strlen (rkw[i].kname); /* Initialize region list to zeroes */ for (i = 0; i < nrmax; i++) rgns[i] = 0; /* Initialize FITS table line buffer to zeroes */ for (i = 0; i < 120; i++) fitsline[i] = (char) 0; nrgn = 0; /* Set pathnames to Guide Star Catalog CDROM */ if (refcat == GSCACT) strcpy (tabpath,cdna); else strcpy (tabpath,cdn); /* Set pathname for region table file */ strcat (tabpath,"/tables/regions.tbl"); /* Open the index table */ kwn = rkw; ift = fitsrtopen (tabpath,&nrkw,&kwn, &nrows, &nchar, &nbhead); /* If the northern hemisphere CDROM cannot be read, try the southern */ if (ift < 0) { if (refcat == GSCACT) strcpy (tabpath,cdsa); else strcpy (tabpath,cds); strcat (tabpath,"/tables/regions.tbl"); ift = fitsrtopen (tabpath,&nrkw,&kwn,&nchar,&nrows,&nbhead); if (ift < 0) { fprintf (stderr,"GSCREG: error reading region table %s\n",tabpath); return (0); } } /* Find region range to search based on declination */ iz1 = gsczone (dec1); iz2 = gsczone (dec2); jr1 = 0; jr2 = 0; nwrap = 1; /* Search region northern hemisphere or only one region */ if (iz2 >= iz1) { ir1 = zreg1[iz1]; ir2 = zreg2[iz2]; } /* Search region in southern hemisphere with multiple regions */ else if (dec1 < 0 && dec2 < 0) { ir1 = zreg1[iz2]; ir2 = zreg2[iz1]; } /* Search region spans equator */ else if (dec1 < 0 && dec2 >= 0) { ir1 = zreg1[12]; ir2 = zreg2[iz1]; jr1 = zreg1[0]; jr2 = zreg2[iz2]; nwrap = 2; } else { ir1 = 1; ir2 = 0; } nsrch = ir2 - ir1 + 1; if (verbose) fprintf (stderr,"GSCREG: searching %d regions: %d - %d\n",nsrch,ir1,ir2); if (jr1 > 0) { nsrch1 = jr2 - jr1 + 1; if (verbose) fprintf (stderr,"GSCREG: searching %d regions: %d - %d\n",nsrch1,jr1,jr2); } if (verbose) fprintf(stderr,"GSCREG: RA: %.5f - %.5f, Dec: %.5f - %.5f\n",ra1,ra2,dec1,dec2); if (nsrch < 1) { return (0); } nrgn = 0; for (iwrap = 0; iwrap < nwrap; iwrap++) { for (irow = ir1 - 1; irow < ir2; irow++) { /* Read next line of region table */ nbr = fitsrtline (ift,nbhead,ltab,table,irow,nchar,fitsline); if (nbr < nchar) { fprintf (stderr,"GSREG: %d / %d bytes read for row %d\n",nbr,nchar,irow); break; } /* Declination range of the gs region */ /* note: southern dechi and declow are reversed */ dechi = gscdec (fitsline, 10, 11, 12); declow = gscdec (fitsline, 7, 8, 9); if (dechi > declow) { decmin = declow; decmax = dechi; } else { decmax = declow; decmin = dechi; } if (decmax >= dec1 && decmin <= dec2) { /* Right ascension range of the Guide Star Catalog region */ ralow = gscra (fitsline, 1, 2, 3); rahi = gscra (fitsline, 4, 5, 6); if (rahi <= 0.0) rahi = 360.0; /* Check RA if 0h RA not between region RA limits */ if (ra1 < ra2) { if (ralow <= ra2 && rahi >= ra1) { /* Get region number from FITS table */ regnum = ftgeti4 (fitsline, &rkw[0]); if (verbose) fprintf (stderr,"GSCREG: Region %d added to search\n",regnum); /* Add this region to list, if there is space */ if (nrgn < nrmax) { rgns[nrgn] = regnum; nrgn = nrgn + 1; } } } /* Check RA if 0h RA is between region RA limits */ else { if (ralow > rahi) rahi = rahi + 360.0; if (ralow <= ra2 || rahi >= ra1) { /* Get region number from FITS table */ regnum = ftgeti4 (fitsline, &rkw[0]); if (verbose) fprintf (stderr,"GSCREG: Region %d added to search\n",regnum); /* Add this region to list, if there is space */ if (nrgn < nrmax) { rgns[nrgn] = regnum; nrgn = nrgn + 1; } } } } } /* Handle wrap-around through the equator */ ir1 = jr1; ir2 = jr2; jr1 = 0; jr2 = 0; } (void) close (ift); return (nrgn); } /* GSCRA -- returns right ascension in degrees from the GSC index table * This is converted from the hours, minutes, and seconds columns. */ double gscra (fitsline, hcol, mcol, scol) char *fitsline; /* Index table line */ int hcol; /* Column index for hours */ int mcol; /* Column index for minutes */ int scol; /* Column index for seconds */ { double ra; /* Right ascension in fractional degrees */ double hrs; /* Hours of right ascension */ double min; /* Minutes of right ascension */ double sec; /* Seconds of right ascension */ /* hours of right ascension */ hrs = ftgetr8 (fitsline, &rkw[hcol]); /* minutes of right ascension */ min = ftgetr8 (fitsline, &rkw[mcol]); /* seconds of right ascension */ sec = ftgetr8 (fitsline, &rkw[scol]); /* right ascension in degrees */ ra = hrs + (min / 60.0) + (sec / 3600.0); ra = ra * 15.0; return (ra); } /* GSCDEC -- returns the declination in degrees from the GSC index table. * This is converted from the sign, degrees, minutes, and seconds columns. */ double gscdec (fitsline, sgncol, dcol, mcol) char *fitsline; /* Index table line */ int sgncol; /* Column index for sign */ int dcol; /* Column index for degrees */ int mcol; /* Column index for minutes */ { double dec; /* Declination in fractional degrees */ char sgn[4]; /* Sign of declination */ double deg; /* Degrees of declination*/ double min; /* Minutes of declination */ /* Get declination sign from table */ (void) ftgetc (fitsline, &rkw[sgncol], sgn, 3); /* Get degrees of declination from table */ deg = ftgetr8 (fitsline, &rkw[dcol]); /* Get minutes of declination from table */ min = ftgetr8 (fitsline, &rkw[mcol]); dec = deg + (min / 60.0); /* Negative declination */ if (strchr (sgn, '-') != NULL) dec = -dec; return (dec); } /* GSCZONE -- find the zone number where a declination can be found */ int gsczone (dec) double dec; /* declination in degrees */ { int zone; /* gsc zone (returned) */ double zonesize; int ndeczones = 12; /* number of declination zones per hemisphere */ /* Width of declination zones */ zonesize = 90.0 / ndeczones; zone = ((int) (dec / zonesize)); if (dec < 0) zone = ndeczones - zone; return (zone); } /* GSCPATH -- Get HST Guide Star Catalog region FITS file pathname */ static void gscpath (refcat, regnum, path) int refcat; /* Catalog code (GSC or GSCACT) */ int regnum; /* Guide Star Catalog region number */ char *path; /* Pathname of GSC region FITS file */ { int zone = 0; /* Name of Guide Star Catalog zone directory */ int i; /* Get zone directory name given region number */ for (i = 0; i < 24; i++) { if (regnum >= zreg1[i] && regnum <= zreg2[i]) { zone = i; break; } } /* Set the pathname using the appropriate GSC CDROM directory */ /* Northern hemisphere disk (volume 1) */ if (regnum < zreg1[13]) { if (refcat == GSCACT) sprintf (path,"%s/%s/%04d.gsc", cdna, zdir[zone], regnum); else sprintf (path,"%s/gsc/%s/%04d.gsc", cdn, zdir[zone], regnum); } /* Southern hemisphere disk (volume 2) */ else { if (refcat == GSCACT) sprintf (path,"%s/%s/%04d.gsc", cdsa, zdir[zone], regnum); else sprintf (path,"%s/gsc/%s/%04d.gsc", cds, zdir[zone], regnum); } return; } /* Feb 16 1996 New program * May 14 1996 Change arguments to FITSRTOPEN * May 17 1996 Fix bug so equal magnitude limits accepts anything * May 17 1996 Use new FITSRTOPEN which internalizes FITS header * May 24 1996 Fix string decoding bug and region search * May 31 1996 Use stream I/O * Jun 10 1996 Remove unused variables after using lint * Jul 1 1996 Fix GSC pathname * Aug 6 1996 Minor changes after lint * Sep 17 1996 Fix bug causing incomplete region access; improve logging * Oct 17 1996 Change FITS table entry calls * 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 * Dec 12 1996 Make logging run on single line per region * Dec 17 1996 Add code to keep brightest or closest stars if too many found * Dec 18 1996 Add code to get star information by number * * Mar 20 1997 Remove unused variables and fixed logging in GSCRNUM after lint * Oct 10 1997 Use standard I/O instead of stream I/O when reading FITS files * Nov 6 1997 Do not print overrun unless logging * * May 13 1998 Do not set classd in gscread() * May 13 1998 Print all stars if classd is < -1 * May 27 1998 Include fitsio.h instead of fitshead.h * Jun 24 1998 Add string lengths to ra2str() and dec2str() calls * Aug 6 1998 Change fitsio.h to fitsfile.h * Sep 16 1998 Use limiting radius, if present * Sep 22 1998 Convert to desired output coordinate system * Oct 26 1998 Fix bug in region selection * Oct 29 1998 Correctly assign numbers when too many stars are found * * May 25 1999 Allocate table buffer only once * Jun 16 1999 Use SearchLim() * Aug 16 1999 Add RefLim() to get converted search coordinates right * Aug 16 1999 Fix bug to fix failure to search across 0:00 RA * Aug 25 1999 Return real number of stars from gscread() * Sep 10 1999 Set class selection with subroutine, not argument * 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 * Sep 22 1999 Rewrite table allocation so it works; make ltab static * Oct 21 1999 Include wcscat.h, unistd.h * * Mar 27 2000 Drop unused variables after lint * Mar 28 2000 Make default to read all classes * Jun 26 2000 Add coordinate system to SearchLim() arguments * Nov 29 2000 Add option to read cataog across the web * Dec 11 2000 Allow search engine URL in cdn[] * Dec 12 2000 Fix wrong web subroutine in gscrnum() * * Apr 24 2001 Add bandpass code and number of entries to object class * Apr 24 2001 Return individual entries if class is < -1 * May 23 2001 Add support for GSC-ACT * May 30 2001 Use GSC_NORTH and GSCACT_NORTH instead of *_PATH for consistency * Jun 27 2001 Print stars as found in gscread() if nstarmax < 1 * Jun 27 2001 Allocate distance array only if larger one is needed * Jun 27 2001 Add gscfree() to free table buffer and distance array * Sep 11 2001 Use single magnitude argument to gscread() and webread() * Sep 21 2001 Clean up web interface * Nov 20 2001 Change cos(degrad)) to cosdeg() * * Mar 28 2002 Change pathnames to /data/astrocat * Jul 31 2002 Always check for GSCACT_PATH, not GSCA_PATH * Oct 2 2002 Fix pass-through tab table header * * Mar 10 2003 Improve position limits * Apr 3 2003 Drop wrap in gscread(); add test for npos in gscread() * Apr 14 2003 Explicitly get revision date if nstarmax < 1 * Aug 22 2003 Add dradi for inner edge of annulus search * Sep 25 2003 Add gscbin() to fill an image with sources * Oct 6 2003 Update gscread() and gscbin() for improved RefLim() * Nov 18 2003 Initialize image size and bits/pixel from header in gscbin() * Dec 1 2003 Add missing tab to n=-1 header * Dec 12 2003 Fix bug in wcs2pix() call in gscbin() * * Aug 27 2004 Include math.h * * Jun 20 2006 Initialize uninitialized variables * Nov 16 2006 Fix binning * * Jan 10 2007 Add match=1 argument to webrnum() * Jan 10 2007 Rewrite web access in gscread() and gscrnum() to reduce code * * Sep 22 2009 Initialize lengths of FITS table columns * Sep 22 2009 Change region table keywords from DEC*LOW to DEC*LO */