/*** File libwcs/actread.c *** January 10, 2007 *** By Jessica Mink, jmink@cfa.harvard.edu *** Copyright (C) 1999-2007 *** 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 #include "fitsfile.h" #include "wcs.h" #include "wcscat.h" /* pathname of ACT CDROM or catalog search engine URL */ char actcd[64]="/data/act"; #define MAXREG 100 static double *gdist; /* Array of distances to stars */ static int ndist = 0; static int actreg(); struct StarCat *actopen(); void actclose(); static int actstar(); static int actsize(); static int actsra(); /* ACTREAD -- Read USNO ACT Star Catalog stars from CDROM */ int actread (cra,cdec,dra,ddec,drad,dradi,distsort,sysout,eqout,epout,mag1,mag2, sortmag,nstarmax,gnum,gra,gdec,gpra,gpdec,gmag,gtype,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 */ double drad; /* Limiting separation in degrees (ignore if 0) */ double dradi; /* Inner edge of annulus in degrees (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 sortmag; /* Magnitude by which to sort (1 to nmag) */ 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 *gpra; /* Array of right ascension proper motions (returned) */ double *gpdec; /* Array of declination proper motions (returned) */ double **gmag; /* Array of visual 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 nreg = 0; /* Number of ACT regions in search */ int rlist[MAXREG]; /* List of input region files */ int sysref=WCS_J2000; /* Catalog coordinate system */ double eqref=2000.0; /* Catalog equinox */ double epref=2000.0; /* Catalog epoch */ double secmarg = 60.0; /* Arcsec/century margin for proper motion */ struct StarCat *starcat; struct Star *star; int verbose; int wrap; int pass; int magsort; int rnum, ireg; int jstar, iw; int nrmax,nstar,i, ntot; int istar, istar1, istar2; double num, ra, dec, rapm, decpm, mag, magb, magv; double rra1, rra2, rra2a, rdec1, rdec2; double rdist, ddist; char *str; char cstr[32], decstr[32], rastr[32]; int nbytes; ntot = 0; if (nlog == 1) verbose = 1; else verbose = 0; /* Set path to ACT Catalog */ if ((str = getenv("ACT_PATH")) == NULL ) str = actcd; /* If pathname is a URL, search and return */ if (!strncmp (str, "http:",5)) { return (webread (str,"act",distsort,cra,cdec,dra,ddec,drad,dradi, sysout,eqout,epout,mag1,mag2,sortmag,nstarmax, gnum,gra,gdec,gpra,gpdec,gmag,gtype,nlog)); } if (sortmag > 0 && sortmag < 3) magsort = sortmag - 1; else magsort = 1; 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; } /* Allocate table for distances of stars from search center */ if (nstarmax > 10) nbytes = nstarmax * sizeof (double); else nbytes = 10 * sizeof (double); if (nstarmax > ndist) { if (ndist > 0) free ((void *) gdist); gdist = (double *) malloc (nbytes); if (gdist == NULL) { fprintf (stderr,"ACTREAD: cannot allocate separation array\n"); return (0); } ndist = nstarmax; } /* Allocate catalog entry buffer */ star = (struct Star *) calloc (1, sizeof (struct Star)); star->num = 0.0; nstar = 0; jstar = 0; rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,secmarg, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); if (wrap) { rra2a = rra2; rra2 = 360.0; } else rra2a = 0.0; nrmax = MAXREG; /* Write header if printing star entries as found */ if (nstarmax < 1) { char *revmessage; revmessage = getrevmsg(); printf ("catalog ACT\n"); ra2str (rastr, 31, cra, 3); printf ("ra %s\n", rastr); dec2str (decstr, 31, cdec, 2); printf ("dec %s\n", decstr); printf ("rpmunit mas/year\n"); printf ("dpmunit mas/year\n"); 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); printf ("act_id ra dec "); printf ("magb magv ura udec arcmin\n"); printf ("---------- ------------ ------------ "); printf ("----- ----- ------ ------ ------\n"); } /* If searching through RA = 0:00, split search in two */ for (iw = 0; iw <= wrap; iw++) { /* Find ACT Star Catalog regions in which to search */ nreg = actreg (rra1,rra2,rdec1,rdec2,nrmax,rlist,verbose); if (nreg <= 0) { fprintf (stderr,"ACTREAD: no ACT regions found\n"); free ((void *)star); return (0); } /* Loop through region list */ for (ireg = 0; ireg < nreg; ireg++) { /* Open catalog file for this region */ starcat = actopen (rlist[ireg]); rnum = rlist[ireg]; /* Set first and last stars to check */ istar1 = actsra (starcat, star, rra1); istar2 = actsra (starcat, star, rra2); if (verbose) fprintf (stderr,"ACTREAD: Searching stars %d.%d through %d.%d\n", rnum,istar1,rnum,istar2); /* Loop through catalog for this region */ for (istar = istar1; istar <= istar2; istar++) { if (actstar (starcat, star, istar)) { fprintf (stderr,"ACTREAD: Cannot read star %d\n", istar); break; } /* Magnitude */ magv = star->xmag[0]; magb = star->xmag[1]; mag = star->xmag[magsort]; /* Check magnitude limits */ pass = 1; if (mag1 != mag2 && (mag < mag1 || mag > mag2)) pass = 0; /* Get position in output coordinate system */ if (pass) { rapm = star->rapm; decpm = star->decpm; ra = star->ra; dec = star->dec; wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); /* 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; } } if (pass) { /* ID number */ num = (double) rlist[ireg] + (star->num / 100000.0); /* Write star position and magnitudes to stdout */ if (nstarmax < 1) { ra2str (rastr, 31, ra, 3); dec2str (decstr, 31, dec, 2); dist = wcsdist (cra,cdec,ra,dec) * 60.0; printf ("%010.5f %s %s", num,rastr,decstr); printf (" %.2f %.2f %6.1f %6.1f %.2f\n", magb, magv, gpra[nstar] * 3600000.0 * cosdeg(dec), gpdec[nstar] * 3600000.0, dist / 60.0); } /* Save star position and magnitude in table */ if (nstar < nstarmax) { gnum[nstar] = num; gra[nstar] = ra; gdec[nstar] = dec; gpra[nstar] = rapm; gpdec[nstar] = decpm; gmag[0][nstar] = magb; gmag[1][nstar] = magv; 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 farthest star */ else if (distsort) { if (dist < maxdist) { gnum[farstar] = num; gra[farstar] = ra; gdec[farstar] = dec; gpra[farstar] = rapm; gpdec[farstar] = decpm; gmag[0][farstar] = magb; gmag[1][farstar] = magv; gdist[farstar] = dist; /* Find new farthest star */ maxdist = 0.0; for (i = 0; i < nstarmax; i++) { if (gdist[i] > maxdist) { maxdist = gdist[i]; farstar = i; } } } } /* Else if too many stars, replace faintest star */ else if (mag < faintmag) { gnum[faintstar] = num; gra[faintstar] = ra; gdec[faintstar] = dec; gpra[farstar] = rapm; gpdec[farstar] = decpm; gmag[0][faintstar] = magb; gmag[1][faintstar] = magv; gdist[faintstar] = dist; faintmag = 0.0; /* Find new faintest star */ for (i = 0; i < nstarmax; i++) { if (gmag[magsort][i] > faintmag) { faintmag = gmag[magsort][i]; faintstar = i; } } } nstar++; if (nlog == 1) fprintf (stderr,"ACTREAD: %11.6f: %9.5f %9.5f %5.2f %5.2f\n", num,ra,dec,magb,mag); /* End of accepted star processing */ } /* Log operation */ jstar++; if (nlog > 0 && istar%nlog == 0) fprintf (stderr,"ACTREAD: %5d / %5d / %5d sources\r", nstar,jstar,starcat->nstars); /* End of star loop */ } ntot = ntot + starcat->nstars; if (nlog > 0) fprintf (stderr,"ACTREAD: %4d / %4d: %5d / %5d / %5d sources from region %4d \n", ireg+1,nreg,nstar,jstar,starcat->nstars,rlist[ireg]); /* Close region input file */ actclose (starcat); } rra1 = 0.0; rra2 = rra2a; } /* close output file and summarize transfer */ if (nlog > 0) { if (nreg > 1) fprintf (stderr,"ACTREAD: %d regions: %d / %d found\n",nreg,nstar,ntot); else fprintf (stderr,"ACTREAD: 1 region: %d / %d found\n",nstar,ntot); if (nstar > nstarmax) fprintf (stderr,"ACTREAD: %d stars found; only %d returned\n", nstar,nstarmax); } free ((void *)star); return (nstar); } /* ACTRNUM -- Read HST Guide Star Catalog stars from CDROM */ int actrnum (nstars,sysout,eqout,epout, gnum,gra,gdec,gpra,gpdec,gmag,gtype,nlog) 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 *gpra; /* Array of right ascension proper motions (returned) */ double *gpdec; /* Array of declination proper motions (returned) */ double **gmag; /* Array of V, B magnitudes (returned) */ int *gtype; /* Array of object types (returned) */ int nlog; /* 1 for diagnostics */ { int sysref; /* Catalog coordinate system */ double eqref; /* Catalog equinox */ double epref; /* Catalog epoch */ struct StarCat *starcat = NULL; struct Star *star; int rnum; int jstar; int istar, nstar, snum; double num, ra, dec, rapm, decpm, mag, magb; char *str; /* Set path to ACT Catalog */ if ((str = getenv("ACT_PATH")) == NULL ) str = actcd; /* If pathname is a URL, search and return */ if (!strncmp (str, "http:",5)) { return (webrnum (str,"act",nstars, sysout,eqout,epout, 1, gnum,gra,gdec,gpra,gpdec,gmag,gtype,nlog)); } /* Allocate catalog entry buffer */ star = (struct Star *) calloc (1, sizeof (struct Star)); star->num = 0.0; nstar = 0; /* Loop through star list */ for (jstar = 0; jstar < nstars; jstar++) { rnum = (int) gnum[jstar]; snum = (int) (((gnum[jstar] - (double)rnum) * 100000.0) + 0.01); /* Open file for this region of ACT catalog */ starcat = actopen (rnum); if (starcat == NULL) { free ((void*) star); return (0); } sysref = starcat->coorsys; eqref = starcat->equinox; epref = starcat->epoch; /* Find star in catalog */ istar = snum; if (actstar (starcat, star, istar)) { fprintf (stderr,"ACTRNUM: Cannot read star %d\n", istar); gra[nstar] = 0.0; gdec[nstar] = 0.0; gmag[0][nstar] = 0.0; gmag[1][nstar] = 0.0; gtype[nstar] = 0; continue; } /* If star has been found in catalog */ /* ID number */ num = star->num; /* Position in degrees at designated epoch */ ra = star->ra; dec = star->dec; rapm = star->rapm; decpm = star->decpm; wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); /* Magnitude */ mag = star->xmag[0]; magb = star->xmag[1]; /* Save star position and magnitude in table */ gra[nstar] = ra; gdec[nstar] = dec; gpra[nstar] = rapm; gpdec[nstar] = decpm; gmag[0][nstar] = magb; gmag[1][nstar] = mag; nstar++; if (nlog == 1) fprintf (stderr,"ACTRNUM: %11.6f: %9.5f %9.5f %5.2f %5.2f \n", num, ra, dec, magb, mag); /* End of star loop */ } /* Summarize search */ if (nlog > 0) fprintf (stderr,"ACTRNUM: %d / %d found\n",nstar,nstars); if (starcat != NULL) actclose (starcat); free ((void*) star); return (nstar); } /* ACTBIN -- Fill FITS WCS image with USNO ACT Star Catalog stars */ int actbin (wcs, header, image, mag1, mag2, sortmag, magscale, nlog) 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; /* 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; /* 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 */ int nreg = 0; /* Number of ACT regions in search */ int rlist[MAXREG]; /* List of input region files */ int sysref=WCS_J2000; /* Catalog coordinate system */ double eqref=2000.0; /* Catalog equinox */ double epref=2000.0; /* Catalog epoch */ double secmarg = 60.0; /* Arcsec/century margin for proper motion */ struct StarCat *starcat; struct Star *star; int verbose; int wrap; int pass; int magsort; int rnum, ireg; int ix, iy; int jstar, iw; int nrmax,nstar, ntot; int istar, istar1, istar2; double num, ra, dec, rapm, decpm, mag, magb, magv; double rra1, rra2, rra2a, rdec1, rdec2; double rdist, ddist; char *str; char cstr[32]; double xpix, ypix, flux; int offscl; int bitpix, w, h; /* Image bits/pixel and pixel width and height */ double logt = log(10.0); ntot = 0; if (nlog == 1) verbose = 1; else verbose = 0; /* Set path to ACT Catalog */ str = getenv("ACT_PATH"); if (sortmag > 0 && sortmag < 3) magsort = sortmag - 1; else magsort = 1; /* 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); /* If RA range includes zero, split it in two */ wrap = 0; if (ra1 > ra2) wrap = 1; else wrap = 0; /* make mag1 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; nstar = 0; jstar = 0; rra1 = ra1; rra2 = ra2; rdec1 = dec1; rdec2 = dec2; RefLim (cra,cdec,dra,ddec,sysout,sysref,eqout,eqref,epout,epref,secmarg, &rra1, &rra2, &rdec1, &rdec2, &wrap, verbose); if (wrap) { rra2a = rra2; rra2 = 360.0; } else rra2a = 0.0; nrmax = MAXREG; /* If searching through RA = 0:00, split search in two */ for (iw = 0; iw <= wrap; iw++) { /* Find ACT Star Catalog regions in which to search */ nreg = actreg (rra1,rra2,rdec1,rdec2,nrmax,rlist,verbose); if (nreg <= 0) { fprintf (stderr,"ACTBIN: no ACT regions found\n"); free ((void *)star); return (0); } /* Loop through region list */ for (ireg = 0; ireg < nreg; ireg++) { /* Open catalog file for this region */ starcat = actopen (rlist[ireg]); rnum = rlist[ireg]; /* Set first and last stars to check */ istar1 = actsra (starcat, star, rra1); istar2 = actsra (starcat, star, rra2); if (verbose) fprintf (stderr,"ACTBIN: Searching stars %d.%d through %d.%d\n", rnum,istar1,rnum,istar2); /* Loop through catalog for this region */ for (istar = istar1; istar <= istar2; istar++) { if (actstar (starcat, star, istar)) { fprintf (stderr,"ACTBIN: Cannot read star %d\n", istar); break; } /* ID number */ num = star->num; /* Magnitude */ magv = star->xmag[0]; magb = star->xmag[1]; mag = star->xmag[magsort]; /* Check magnitude limits */ pass = 1; if (mag1 != mag2 && (mag < mag1 || mag > mag2)) pass = 0; /* Get position in output coordinate system */ if (pass) { rapm = star->rapm; decpm = star->decpm; ra = star->ra; dec = star->dec; wcsconp (sysref, sysout, eqref, eqout, epref, epout, &ra, &dec, &rapm, &decpm); /* 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 * (-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,"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 */ jstar++; if (nlog > 0 && istar%nlog == 0) fprintf (stderr,"ACTBIN: %5d / %5d / %5d sources\r", nstar,jstar,starcat->nstars); /* End of star loop */ } ntot = ntot + starcat->nstars; if (nlog > 0) fprintf (stderr,"ACTBIN: %4d / %4d: %5d / %5d / %5d sources from region %4d \n", ireg+1,nreg,nstar,jstar,starcat->nstars,rlist[ireg]); /* Close region input file */ actclose (starcat); } rra1 = 0.0; rra2 = rra2a; } /* close output file and summarize transfer */ if (nlog > 0) { if (nreg > 1) fprintf (stderr,"ACTBIN: %d regions: %d / %d found\n",nreg,nstar,ntot); else fprintf (stderr,"ACTBIN: 1 region: %d / %d found\n",nstar,ntot); } free ((void *)star); return (nstar); } /* ACTREG -- from RA and Dec ranges, figure out which ACT files to search * Build a list containing the numeric part of the CDROM file names. */ static int regions[48]={ 0, 30, 100, 130, 200, 230, 300, 330, 400, 430, 500, 530, 600, 630, 700, 730, 800, 830, 900, 930, 1000,1030,1100,1130,1200,1230,1300,1330,1400,1430, 1500,1530,1600,1630,1700,1730,1800,1830,1900,1930, 2000,2030,2100,2130,2200,2230,2300,2330}; static double reghour[49]={ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5,14.0,14.5, 15.0,15.5,16.0,16.5,17.0,17.5,18.0,18.5,19.0,19.5, 20.0,20.5,21.0,21.5,22.0,22.5,23.0,23.5,24.0}; static int actreg (ra1, ra2, dec1, dec2, nrmax, rgns, verbose) double ra1, ra2; /* Right ascension limits in degrees */ double dec1, dec2; /* Declination limits in degrees */ int nrmax; /* Maximum number of regions to find */ int *rgns; /* Region numbers (returned)*/ int verbose; /* 1 for diagnostics */ { int nsrch; /* Number of regions found (returned) */ int i, ir, irx; int ir1 = 0; int ir2 = 0; /* Zero out regions to be searched */ for (i = 0; i < nrmax; i++) rgns[i] = 0; /* Find region range to search based on declination */ /* Find RA regions to search */ ra1 = ra1 / 15.0; ra2 = ra2 / 15.0; irx = 0; /* Find first region to search */ for (ir = 1; ir < 49; ir++) { if (ra1 >= reghour[ir-1] && ra1 <= reghour[ir]) { ir1 = ir - 1; break; } } /* Find last region to search */ for (ir = 1; ir < 49; ir++) { if (ra2 >= reghour[ir-1] && ra2 <= reghour[ir]) { ir2 = ir - 1; break; } } if (ir2 >= ir1) { for (ir = ir1; ir <= ir2; ir++) { if (irx < nrmax) rgns[irx++] = regions[ir]; } } else if (ir2 < ir1) { for (ir = ir1; ir < 48; ir++) { if (irx < nrmax) rgns[irx++] = regions[ir]; } for (ir = 0; ir <= ir2; ir++) { if (irx < nrmax) rgns[irx++] = regions[ir]; } } nsrch = irx; if (verbose) { fprintf (stderr,"ACTREG: RA: %.5f - %.5f, Dec: %.5f - %.5f\n", ra1,ra2,dec1,dec2); fprintf (stderr,"ACTREG: searching %d regions:",nsrch); for (ir = 0; ir < nsrch; ir++) fprintf (stderr," %04d",rgns[ir]); fprintf (stderr,"\n"); } return (nsrch); } /* ACTOPEN -- Open ACT catalog region file, returning number of entries */ struct StarCat * actopen (regnum) int regnum; /* ACT Catalog region number */ { FILE *fcat; struct StarCat *sc; int lfile, lpath; char *actfile; char *path; /* Full pathname for catalog file */ char *cdpath; /* Set the pathname using the appropriate ACT CDROM directory */ if ((cdpath = getenv("ACT_PATH")) == NULL ) cdpath = actcd; lpath = strlen (cdpath) + 32; path = (char *) calloc (lpath, 1); /* Declination zoned regions */ if (regnum > 0 && regnum < 5) sprintf (path,"%s/data2/act%1d.dat", cdpath, regnum); /* Right ascension zoned regions */ else sprintf (path,"%s/data1/act%04d.dat", cdpath, regnum); /* Find length of ACT catalog region file */ lfile = actsize (path); /* Check for existence of catalog */ if (lfile < 2) { fprintf (stderr,"ACTOPEN: Binary catalog %s has no entries\n", path); free (path); return (0); } /* Open ACT region file */ if (!(fcat = fopen (path, "r"))) { fprintf (stderr,"ACTOPEN: ACT region file %s cannot be read\n",path); free (path); return (0); } /* Set ACT catalog header information */ sc = (struct StarCat *) calloc (1, sizeof (struct StarCat)); sc->byteswapped = 0; sc->nbent = 161; sc->nstars = lfile / sc->nbent; /* Separate filename from pathname and save in structure */ actfile = strrchr (path,'/'); if (actfile) actfile = actfile + 1; else actfile = path; if (strlen (actfile) < 24) strcpy (sc->isfil, actfile); else strncpy (sc->isfil, actfile, 23); /* Set other catalog information in structure */ sc->inform = 'J'; sc->coorsys = WCS_J2000; sc->epoch = 2000.0; sc->equinox = 2000.0; sc->ifcat = fcat; sc->sptype = 2; /* ACT region files are all RA-sorted */ sc->rasorted = 1; return (sc); } void actclose (sc) struct StarCat *sc; /* Star catalog descriptor */ { fclose (sc->ifcat); free ((void *)sc); return; } /* ACTSRA -- Find star closest to given RA in ACT catalog file */ static int actsra (sc, st, dra) struct StarCat *sc; /* Star catalog descriptor */ struct Star *st; /* Current star entry */ double dra; /* Right ascension in degrees */ { char rastr[32], raxstr[32], ramins[32], ramaxs[32]; int istar0, istarx, nrep, ismin, ismax; double rax, ramin, ramax; int verbose = 0; /* Keep RA between 0 and 360 degrees */ if (dra > 360.0) rax = dra - 360.0; else rax = dra; ismin = 1; if (actstar (sc, st, ismin)) { fprintf (stderr,"ACTSRA: Cannot read star %d\n", ismin); return (0); } else ramin = st->ra; ismax = sc->nstars; if (actstar (sc, st, ismax)) { fprintf (stderr,"ACTSRA: Cannot read star %d\n", ismax); return (0); } else ramax = st->ra; istarx = sc->nstars / 2; for (nrep = 0; nrep < 32; nrep++) { if (actstar (sc, st, istarx)) { fprintf (stderr,"ACTSRA: Cannot read star %d\n", istarx); return (0); } /* Find next catalog number to read */ if (st->ra < rax) { ismin = istarx; ramin = st->ra; istar0 = istarx; if (ismax - istarx > 1) istarx = istarx + (ismax - istarx) / 2; else if (ismax - istarx > 0) istarx = istarx + 1; } else if (st->ra > rax) { ismax = istarx; ramax = st->ra; istar0 = istarx; if (istarx - ismin > 1) istarx = istarx - ((istarx - ismin) / 2); else if (istarx - ismin > 0) istarx = istarx - 1; } else break; if (verbose) { ra2str (rastr, 16, st->ra, 3); ra2str (raxstr, 16, rax, 3); ra2str (ramins, 16, ramin, 3); ra2str (ramaxs, 16, ramax, 3); fprintf (stderr,"%9d: %s -> %s %9d: %s %9d: %s\n", istarx, rastr, raxstr, ismin,ramins,ismax,ramaxs); } if (istarx == istar0) break; } /* Make sure final star is real */ if (actstar (sc, st, istarx)) { fprintf (stderr,"ACTSRA: Cannot read star %d\n", istarx); return (0); } else return (istarx); } /* ACTSTAR -- Get ACT catalog entry for one star; return 0 if successful */ static int actstar (sc, st, istar) struct StarCat *sc; /* Star catalog descriptor */ struct Star *st; /* Current star entry */ int istar; /* Star sequence number in ACT catalog region file */ { int nbr; long offset; char dsgn; char line[256]; int irh,irm,idd,idm; double rs, ds; /* double bvmag; */ /* Drop out if catalog pointer is not set */ if (sc == NULL) return (1); /* Drop out if catalog is not open */ if (sc->ifcat == NULL) return (2); /* Drop out if star number is too large */ if (istar > sc->nstars) { fprintf (stderr, "ACTSTAR: %d > %d is not in catalog\n", istar, sc->nstars); return (3); } /* Move file pointer to start of correct star entry */ if (istar > 0) { offset = (istar - 1) * sc->nbent; if (fseek (sc->ifcat, offset, SEEK_SET)) return (4); } /* Read catalog entry */ if ((nbr = fread (line, sc->nbent, 1, sc->ifcat)) > sc->nbent) { fprintf (stderr, "ACTSTAR: %d / %d bytes read from %s\n", nbr, sc->nbent, sc->isfil); return (5); } st->num = (double) istar; /* Read position for this star */ irh = atoi (line); irm = atoi (line+3); rs = atof (line+6); dsgn = line[14]; idd = atoi (line+15); idm = atoi (line+18); ds = atof (line+21); /* Convert position to degrees */ st->ra = hrdeg ((double)irh + ((double)irm)/60.0 + rs / 3600.0); st->dec = (double) idd + ((double)idm) / 60.0 + ds / 3600.0; if (dsgn == '-') st->dec = -st->dec; /* Read proper motion and convert it to to degrees/year */ st->rapm = hrdeg (atof (line+28)) / 3600.0; st->decpm = atof (line+36) / 3600.0; /* Set V, B, B-V magnitudes */ st->xmag[0] = atof (line+75); st->xmag[1] = atof (line+68); st->xmag[2] = atof (line+82); st->isp[0] = (char) 0; st->isp[1] = (char) 0; return (0); } /* ACTSIZE -- return size of one ACT catalog file in bytes */ static int actsize (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); } /* Feb 11 1999 New program * Apr 13 1999 Fix bugs which caused failure on crossing 0:00 h * May 12 1999 Fix bug for all searches * May 21 1999 Fix bug with proper motion so it is in deg/yr, not sec/yr * Jun 16 1999 Use SearchLim() * Aug 16 1999 Add RefLim() to get converted search coordinates right * Aug 25 1999 Return real number of stars from actread() * 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 Delete unused varaiables after lint * * Jan 5 2000 Add 2 to spectral type string so there can be proper termination * Mar 15 2000 Add proper motions to returns from actread() and actrnum() * May 31 2000 Get spectral type from bv2sp() * Jun 2 2000 Free all allocated data structures * Jun 9 2000 Fix bug which caused memory overflow if limiting number * Jun 26 2000 Add coordinate system to SearchLim() arguments * Sep 25 2000 Set sc->sptype to 2 to indicate presence of spectral type * Nov 29 2000 Add option to read catalog using HTTP * Dec 11 2000 Allow catalog search engine URL in actcd[] * * Jan 11 2001 All printing is to stderr * Jun 14 2001 Drop spectral type approximation * Sep 11 2001 Change to single magnitude argeument * Sep 11 2001 Add sort magnitude argument to actread() * * Apr 8 2002 Fix extraneous declaration of actsize() * Oct 2 2002 Print stars as found in actread() if nstarmax < 1 * * Mar 11 2003 Fix position limit testing * Apr 3 2003 Drop unused type variables in actstar() and actread() * Apr 14 2003 Explicitly get revision date if nstarmax < 1 * Jun 2 2003 Print proper motion as mas/year * Aug 22 2003 Add radi argument for inner edge of search annulus * Aug 22 2003 Fix bug in actrnum: mags were not in one array * Sep 25 2003 Add actbin() to fill an image with sources * Oct 6 2003 Update actread() and actbin() for improved RefLim() * Nov 18 2003 Initialize image size and bits/pixel from header in actbin() * Dec 1 2003 Add missing tab to n=-1 header * Dec 12 2003 Fix bug in wcs2pix() call in actbin() * * Aug 30 2004 Include math.h * * Jun 20 2006 Initialize uninitialized variables * Sep 26 2006 Increase length of rastr and destr from 16 to 32 * Nov 16 2006 Fix binning * * Jan 10 2007 Add match=1 argument to webrnum() * Jan 10 2007 Rewrite web access in actread() and actrnum() to reduce code */