**
** ROUTINE: taUSNOMax
**
**
** Determine the maximum number of stars from USNO-B 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) */,
double minNu /* IN: minimum nu scanned (degrees) */,
double maxNu /* IN: maximum nu scanned (degrees) */,
const char *catDir /* IN: directory containing USNO 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;
int zone1, zone2, raDec, chunk1, chunk2, chunk;
FILE *fp;
double chunkRadius, overlapRadius;
/* Are limits in ra/dec, as indicated by incl=0 (speeds things if they are)*/
if (fabs(incl) < 0.0000001)
raDec = 1;
else
raDec = 0;
/* Each chunk is 3.75 degrees in RA by 0.1 degrees in dec, for a maximum
* radius of 1.9 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. + 1.9 of the stripe great circle. */
chunkRadius = sqrt(zoneWidth * zoneWidth / 4. + chunkSize * chunkSize / 4);
overlapRadius = stripeWidth / 2. + chunkRadius; /* degrees */
/* Loop over USNO zones. The zones are each 0.1 deg wide in south
* polar distance (SPD), starting with 0 < SPD < 0.1. */
if (raDec)
{
zone1 = (int) ((minNu + 90.) / zoneWidth);
zone2 = (int) ((maxNu + 90.) / zoneWidth);
chunk1 = (int) (minMu / chunkSize);
chunk2 = (int) (maxMu / chunkSize);
if (chunk2 < chunk1) chunk2 += nChunks;
}
else
{
zone1 = 0;
zone2 = nZones - 1;
chunk1 = 0;
chunk2 = nChunks - 1;
}
for (zone = zone1; zone <= zone2; zone++)
{
first = 1;
/* For each 3.75 degree chunk (as divided up in the acc file) ... */
for (chunk = chunk1; chunk <= chunk2; chunk++)
{
ira = chunk;
if (ira >= nChunks) ira -= nChunks;
if (! raDec)
{
/* 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 (mu + overlapRadius < minMu) mu = mu + 360.;
if (mu - overlapRadius > maxMu) continue;
if (maxNu > minNu)
{
if (nu + overlapRadius < minNu) continue;
if (nu - overlapRadius > maxNu) continue;
}
else
{
if (fabs(nu) > 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/%03d/b%04d.acc", catDir, zone / 10, zone);
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;
}
/*============================================================================
**