/***************************************************************************** ****************************************************************************** ** ** ** FILE: usno.c ** ** ** This file contains functions to read the USNOA2.0 catalog. ** ** ** ** ** ENVIRONMENT: ** C ** ****************************************************************************** ****************************************************************************** */ #include #include "taCatalogs.h" #define nChunks 96 /* Number of chunks per zone */ #define nMaxPerRead 100000 /* Maximum objects to read at once */ #define arrayIncr 100000 /* Number of objects per memory allocation */ /* Each zone is divided into 96 chunks, 3.75 degrees wide in ra, by the * accelerator file. */ #define zoneWidth 7.5 /* Width of each zone in dec (deg) */ #define chunkSize 3.75 /* Width of each chunk in ra (deg) */ #define stripeWidth 3.0 /* Generous width of each scan in nu (deg) */ /*============================================================================ ** ** ** ROUTINE: taUSNOMax ** ** ** Determine the maximum number of stars from USNOA2.0 that will intersect ** a strip, 3 degrees wide, along a section of a great circle (equinox ** J2000). ** This reads the accelerator files only, counting the objects within each ** chunk that intersects the scan, without reading the actual objects and ** determining which objects in the chunk don't intersect. **

** Returns the maximum number of intersecting objects, or -1 on error. ** ** ** **============================================================================ */ int taUSNOMax(double node /* IN: ascending node of GC (degrees) */, double incl /* IN: inclination of GC (degrees) */, double minMu /* IN: minimum mu scanned (degrees) */, double maxMu /* IN: maximum mu scanned (degrees) */, const char *catDir /* IN: directory containing USNOA2.0 files */) { int zone, ira, i; double ra, dec, mu, nu; int first; /* Have we read the accelerator file yet? */ char fileName[300]; int frec[nChunks], nrec[nChunks]; int nObjs = 0; FILE *fp; /* Each chunk is 3.75 degrees in RA by 7.5 degrees in dec, for a maximum * radius of 4.2 degrees. Each scan is about * 2.5 degrees wide, which we'll generously treat as 3.0 deg. Thus, * we will read a chunk if it's center lies within 3.0/2. + 4.2 of the * stripe great circle. */ double chunkRadius = sqrt(zoneWidth * zoneWidth / 4. + chunkSize * chunkSize / 4); double overlapRadius = stripeWidth / 2. + chunkRadius; /* degrees */ /* Loop over USNOA2.0 zones. The zones are each 7.5 deg wide in south * polar distance (SPD), starting with 0 < SPD < 7.5. */ for (zone = 0; zone < 24; zone++) { first = 1; /* For each 3.75 degree chunk (as divided up in the acc file) ... */ for (ira = 0; ira < nChunks; ira++) { /* Compute chunk center in great circle coordinates. */ ra = (ira + 0.5) * chunkSize; dec = (zone + 0.5) * zoneWidth - 90.; atEqToGC(ra, dec, &mu, &nu, node, incl); /* Skip chunk if it doesn't overlap the scan. */ if (fabs(nu) > overlapRadius) continue; if (mu < minMu - overlapRadius) mu = mu + 360.; if (mu > maxMu + overlapRadius) continue; /* It overlaps the scan. Read it in. */ /* If first chunk, read the accelerator file. */ if (first == 1) { /* Read the accelerator file */ sprintf(fileName, "%s/zone%04d.acc", catDir, zone * 75); if ((fp = fopen(fileName, "r")) == (FILE *)NULL) { shErrStackPush("taUSNOMax: unable to open acc file %s\n", fileName); return -1; } for (i = 0; i < nChunks; i++) if (fscanf(fp, "%*lf %d %d", &(frec[i]), &(nrec[i])) != 2) { fclose(fp); shErrStackPush("taUSNOMax: error reading acc file %s\n", fileName); return -1; } fclose (fp); first = 0; } nObjs += nrec[ira]; } } return nObjs; } /*============================================================================ ** ** ** ROUTINE: taUSNORead ** ** ** Read in that portion of USNOA2.0 that intersects a scan, as ** defined by a strip, 3 degrees wide, along a section of a great circle. ** ** ** **============================================================================ */ RET_CODE taUSNORead(double node /* IN: ascending node of GC (degrees) */, double incl /* IN: inclination of GC (degrees) */, double equinox /* IN: equinox of GC coords (years) */, double minMu /* IN: minimum mu scanned (degrees) */, double maxMu /* IN: maximum mu scanned (degrees) */, const char *catDir /* IN: directory containing USNOA2.0 files */, TA_USNO **array /* OUT: array of USNO objects */, int *nObj /* OUT: number of objects read */) { int zone; int arraySize; /* Equinox must be J2000 */ if (fabs(equinox - 2000.) > 0.00001) { shErrStackPush("taUSNORead: equinox not J2000"); return SH_GENERIC_ERROR; } /* Determine the maximum number of objects which intersect the scan, so * that we can make one memory allocation for that amount. */ if ((arraySize = taUSNOMax(node, incl, minMu, maxMu, catDir)) == -1) return SH_GENERIC_ERROR; /* Initialize the output array of USNO objects. */ if (*array != NULL) { shErrStackPush("taUSNORead: array must be NULL on entry"); return SH_GENERIC_ERROR; } if ((*array = (TA_USNO*) shMalloc(arraySize * sizeof(TA_USNO))) == NULL) { shErrStackPush("taUSNORead: couldn't allocate memory"); return SH_GENERIC_ERROR; } *nObj = 0; /* Loop over USNOA2.0 zones. The zones are each 7.5 deg wide in south * polar distance (SPD), starting with 0 < SPD < 7.5. */ for (zone = 0; zone < 24; zone++) { if (taUSNOZoneRead(node, incl, minMu, maxMu, catDir, zone, array, &arraySize, nObj) != SH_SUCCESS) { if (*array != NULL) shFree(*array); return SH_GENERIC_ERROR; }; } /* Resize arrays to final size */ if (*nObj < arraySize) { if ((*array = (TA_USNO*) shRealloc(*array, *nObj * sizeof(TA_USNO))) == NULL) { shErrStackPush("taUSNORead: couldn't reallocate memory"); return SH_GENERIC_ERROR; } } /* Return */ return SH_SUCCESS; } /*============================================================================ ** ** ** ROUTINE: taUSNOZoneRead ** ** ** Read in that portion of a USNOA2.0 zone that intersects a scan, as ** defined by a strip, 3 degrees wide, along a section of a great circle. ** The great circle definition is assumed to be for equinox J2000. ** ** ** **============================================================================ */ RET_CODE taUSNOZoneRead(double node /* IN: ascending node of GC (degrees) */, double incl /* IN: inclination of GC (degrees) */, double minMu /* IN: minimum mu scanned (degrees) */, double maxMu /* IN: maximum mu scanned (degrees) */, const char *catDir /* IN: directory containing USNO files */, int zone /* IN: USNOA2.0 SPD zone */, TA_USNO **array /* MOD: array of USNO objects */, int *arraySize /* MOD: number of objects allocated */, int *nObj /* MOD: number of objects read */) { int ira, i; double ra, dec, lambda, eta, mu, nu; int first = 1; /* Have we read the accelerator file yet? */ char fileName[300]; FILE *fp; int fd, nToGo, nToRead; int frec[nChunks], nrec[nChunks]; long pos; int buf[nMaxPerRead][3], mags; char *swapBuf, b; size_t nBytes; const int recSize = 12; /* Size of each object record (bytes) */ /* Each chunk is 3.75 degrees in RA by 7.5 degrees in dec, for a maximum * radius of 4.2 degrees. Each scan is about * 2.5 degrees wide, which we'll generously treat as 3.0 deg. Thus, * we will read a chunk if it's center lies within 3.0/2. + 4.2 of the * stripe great circle. */ double chunkRadius = sqrt(zoneWidth * zoneWidth / 4. + chunkSize * chunkSize / 4); double overlapRadius = stripeWidth / 2. + chunkRadius; /* degrees */ /* For each 3.75 degree chunk (as divided up in the accelerator file) ... */ for (ira = 0; ira < nChunks; ira++) { /* Compute chunk center in great circle coordinates. */ ra = (ira + 0.5) * chunkSize; dec = (zone + 0.5) * zoneWidth - 90.; atEqToGC(ra, dec, &mu, &nu, node, incl); /* Skip chunk if it doesn't overlap the scan. */ if (fabs(nu) > overlapRadius) continue; if (mu < minMu - overlapRadius) mu = mu + 360.; if (mu > maxMu + overlapRadius) continue; /* It overlaps the scan. Read it in. */ /* If first chunk, read the accelerator file and open the data file. */ if (first == 1) { /* Read the accelerator file */ sprintf(fileName, "%s/zone%04d.acc", catDir, zone * 75); if ((fp = fopen(fileName, "r")) == (FILE *)NULL) { shErrStackPush("taUSNOZoneRead: unable to open acc file %s\n", fileName); return SH_GENERIC_ERROR; } for (i = 0; i < nChunks; i++) if (fscanf(fp, "%*lf %d %d", &(frec[i]), &(nrec[i])) != 2) { fclose(fp); shErrStackPush("taUSNOZoneRead: error reading acc file %s\n", fileName); return SH_GENERIC_ERROR; } fclose (fp); /* Open the data file */ sprintf(fileName, "%s/zone%04d.cat", catDir, zone * 75); if ((fd = open(fileName, 0)) == -1) { shErrStackPush("taUSNOZoneRead: unable to open cat file %s\n", fileName); return SH_GENERIC_ERROR; } first = 0; } /* Go to the start of the chunk */ pos = (frec[ira] - 1) * recSize; /* Accelerator files are 1-indexed */ if (lseek(fd, pos, SEEK_SET) != pos) { shErrStackPush("taUSNOZoneRead: unable to lseek cat file %s\n", fileName); return SH_GENERIC_ERROR; } /* Read the objects, nMaxPerRead at a time */ nToGo = nrec[ira]; while (nToGo > 0) { nToRead = (nToGo <= nMaxPerRead ? nToGo : nMaxPerRead); nBytes = nToRead * recSize; if (read(fd, buf, nBytes) != nBytes) { shErrStackPush("taUSNOZoneRead: too few bytes read"); return SH_GENERIC_ERROR; } #ifdef SDSS_LITTLE_ENDIAN /* USNOA stored in BIG_ENDIAN order. Byte swap if necessary */ swapBuf = (char*) &(buf[0][0]); for (i = 0; i < nBytes; i += 4) { b = swapBuf[i]; swapBuf[i] = swapBuf[i+3]; swapBuf[i+3] = b; b = swapBuf[i+1]; swapBuf[i+1] = swapBuf[i+2]; swapBuf[i+2] = b; } #endif /* Copy the objects to the array */ for (i = 0; i < nToRead; i++) { /* Skip objects which are more than 1.5 degrees beyond stripe */ ra = buf[i][0] / 360000.; dec = buf[i][1] / 360000. - 90.; atEqToGC(ra, dec, &mu, &nu, node, incl); if (fabs(nu) > stripeWidth / 2.) continue; /* Reallocate memory if necessary */ if (*nObj == *arraySize) { *arraySize += arrayIncr; if ((*array = (TA_USNO*) shRealloc(*array,*arraySize*sizeof(TA_USNO))) == NULL) { shErrStackPush("taUSNOZoneRead: couldn't reallocate memory"); return SH_GENERIC_ERROR; } } /* Reformat mags bit-mask to add 1000 to all south fields, and * get rid of the sign and error flags */ mags = abs(buf[i][2]) % 1000000000; if (zone < 10) if (mags / 1000000 <= 606) mags += 1000000000; /* Add object to the array of potentially matching objects. */ atEqToSurvey(ra, dec, &lambda, &eta); (*array)[*nObj].ra = buf[i][0]; (*array)[*nObj].spd = buf[i][1]; (*array)[*nObj].mags = mags; (*array)[*nObj].lambda = (int) (lambda * 360000. + 0.5); (*nObj)++; } nToGo -= nToRead; } } /* Close up */ if (first == 0) close(fd); return SH_SUCCESS; }