/*** File libwcs/uacread.c
 *** January 11, 2007
 *** By Jessica Mink, jmink@cfa.harvard.edu
 *** Harvard-Smithsonian Center for Astrophysics
 *** Copyright (C) 1996-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

 * Subroutines to read from the USNO A and SA catalogs
 */

#include <unistd.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "wcs.h"
#include "fitsfile.h"
#include "wcscat.h"

static int ucat=UA2;

/* USNO A-2.0 directory pathname; replaced by UA2_PATH environment variable
 * Use this if CDROMs have been transferred to a single hard disk
 * Otherwise set to null string ("") and use cdroot
 * This may also be a URL to a catalog search engine */
static char ua2path[64]="/data/astrocat/ua2";

/* Uncomment following line to use ESO USNO-A server for UA2
static char ua2path[64]="http://archive.eso.org/skycat/servers/usnoa-server";
 */

/* USNO SA-2.0 directory pathname; replaced by USA2_PATH environment variable
 * This may also be a URL to a catalog search engine */
static char usa2path[64]="/data/astrocat/usnosa20";

/* USNO SA-1.0 directory pathname; replaced by USA1_PATH environment variable
 * This may also be a URL to a catalog search engine */
static char usa1path[64]="/data/astrocat/usnosa10";

/* USNO A-1.0 directory pathname; replaced by UA1_PATH environment variable
 * Use this if CDROMs have been transferred to a single hard disk
 * Otherwise set to null string ("") and use cdroot
 * This may also be a URL to a catalog search engine */
static char ua1path[64]="/data/astrocat/ua1";

static char *uapath;

/* Root directory for CDROMs; replaced by UA_ROOT environment variable */
/* Ignored if uapath or UA*_PATH are set */
static char cdroot[32]="/cdrom";

/* Names of CDROM's for USNO A Catalogs */
static char cdname[11][8]={"ua001","ua002","ua003","ua004","ua005","ua006",
			"ua007","ua008","ua009","ua010","ua011"};

/* Disks for 24 zones of USNO A-1.0 Catalog */
static int zdisk1[24]={1,1,6,5,3,2,1,4,6,5,7,10,8,7,8,9,9,4,10,3,2,6,2,3};

/* Disks for 24 zones of USNO A-2.0 Catalog */
static int zdisk2[24]={1,1,9,7,5,4,3,2,1,6,7,10,9,8,8,11,10,11,6,4,2,3,3,2};

typedef struct {
    int rasec, decsec, magetc;
} UACstar;

static int nstars;	/* Number of stars in catalog */
static int cswap = 0;	/* Byte reverse catalog to Intel/DEC order if 1 */
static double *udist;	/* Array of distances to stars */
static int ndist = 0;

static FILE *fcat;
#define ABS(a) ((a) < 0 ? (-(a)) : (a))
#define NZONES 24

static double uacra();
static double uacdec();
static double uacmagr();
static double uacmagb();
int uacmagerr();
int uacgsc();
static int uacplate();

static int uaczones();
static int uaczone();
static int uacsra();
static int uacopen();
static int uacpath();
static int uacstar();
static void uacswap();

static int xplate = 0;	/* If nonzero, use objects only from this plate */
void setuplate (xplate0)
int xplate0;
{ xplate = xplate0; return; }
int getuplate ()
{ return (xplate); }


/* USACREAD -- Read USNO SA Catalog stars from CDROM */

int
usaread (cra,cdec,dra,ddec,drad,dradi,distsort,sysout,eqout,epout,mag1,mag2,
	 sortmag,nstarmax,unum,ura,udec,umag,uplate,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 */
			/* If < 1, print stars to stdout as found */
double	*unum;		/* Array of UA numbers (returned) */
double	*ura;		/* Array of right ascensions (returned) */
double	*udec;		/* Array of declinations (returned) */
double	**umag;		/* Array of red and blue magnitudes (returned) */
int	*uplate;	/* Array of plate numbers (returned) */
int	nlog;		/* Logging interval */
{
    int i;
    int uacread();

    ucat = USA2;
    i = uacread ("usac",distsort,cra,cdec,dra,ddec,drad,dradi,
		 sysout,eqout,epout,mag1,mag2,sortmag,nstarmax,
		 unum,ura,udec,umag,uplate,nlog);
    ucat = USA2;
    return (i);
}


/* UACREAD -- Read USNO A or SA Catalog stars from CDROM */

int
uacread (refcatname,distsort,cra,cdec,dra,ddec,drad,dradi,sysout,eqout,epout,
	 mag1,mag2,sortmag,nstarmax,unum,ura,udec,umag,uplate,nlog)

char	*refcatname;	/* Name of catalog (UAC, USAC, UAC2, USAC2) */
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) */
double	dradi;		/* Inner edge of annulus 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	sortmag;	/* Magnitude by which to sort (1 or 2) */
int	nstarmax;	/* Maximum number of stars to be returned */
double	*unum;		/* Array of UA numbers (returned) */
double	*ura;		/* Array of right ascensions (returned) */
double	*udec;		/* Array of declinations (returned) */
double	**umag;		/* Array of red and blue magnitudes (returned) */
int	*uplate;	/* Array of plate numbers (returned) */
int	nlog;		/* Logging interval */
{
    double ra1,ra2;	/* Limiting right ascensions of region in degrees */
    double dec1,dec2;	/* Limiting declinations of region in degrees */
    int nz;		/* Number of input UA zone files */
    int zlist[NZONES];	/* List of input UA zones */
    UACstar star;	/* UA catalog entry for one star */
    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 sysref=WCS_J2000;	/* Catalog coordinate system */
    double eqref=2000.0;	/* Catalog equinox */
    double epref=2000.0;	/* Catalog epoch */

    double rra1, rra2, rdec1, rdec2;
    double num;		/* UA numbers */
    int wrap, iwrap;
    int verbose;
    int znum, itot,iz, i;
    int jtable,jstar;
    int itable = 0;
    int nstar, nread;
    int uara1, uara2, uadec1, uadec2;
    double ra,dec, rdist, ddist;
    double mag, magb, magr;
    int istar, istar1, istar2, plate;
    int nzmax = NZONES;	/* Maximum number of declination zones */
/*    int isp;
    char ispc[2]; */
    int pass;
    int magsort;
    char *str;
    char cstr[32], rastr[32], numstr[32], decstr[32], catid[32];
    char *title;

    itot = 0;
    if (nlog > 0)
	verbose = 1;
    else
	verbose = 0;

    /* Set catalog code and path to catalog */
    if (strncmp (refcatname,"us",2)==0 ||
        strncmp (refcatname,"US",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("USA2_PATH")) != NULL)
		strcpy (usa2path,str);
	    ucat = USA2;
	    uapath = usa2path;
	    }
	else {
	    if ((str = getenv("USA1_PATH")) != NULL)
		strcpy (usa1path,str);
	    ucat = USA1;
	    uapath = usa1path;
	    }
	}
    else if (strncmp (refcatname,"ua",2)==0 ||
        strncmp (refcatname,"UA",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("UA2_PATH")) != NULL)
		strcpy (ua2path,str);
	    else if ((str = getenv("UA2_ROOT")) != NULL) {
		ua2path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA2;
	    uapath = ua2path;
	    }
	else {
	    if ((str = getenv("UA1_PATH")) != NULL)
		strcpy (ua1path,str);
	    else if ((str = getenv("UA1_ROOT")) != NULL) {
		ua1path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA1;
	    uapath = ua1path;
	    }
	}
    else {
	fprintf (stderr, "UACREAD:  %s not a USNO catalog\n", refcatname);
	return (0);
	}

    /* If root pathname is a URL, search and return */
    if (!strncmp (uapath, "http:",5)) {
	return (webread (uapath,refcatname,distsort,cra,cdec,dra,ddec,drad,
			 dradi,sysout,eqout,epout,mag1,mag2,sortmag,nstarmax,
			 unum,ura,udec,NULL,NULL,umag,uplate,nlog));
	}

    wcscstr (cstr, sysout, eqout, epout);

    SearchLim (cra,cdec,dra,ddec,sysout,&ra1,&ra2,&dec1,&dec2,verbose);

    /* mag1 is always the smallest magnitude */
    if (mag2 < mag1) {
	mag = mag2;
	mag2 = mag1;
	mag1 = mag;
	}
    if (sortmag == 1)
	magsort = 0;
    else
	magsort = 1;

    /* Find UA Star Catalog regions in which to search */
    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);
    nz = uaczones (rra1, rra2, rdec1, rdec2, nzmax, zlist, verbose);
    if (nz <= 0) {
	fprintf (stderr, "UACREAD:  no USNO A zones found\n");
	return (0);
	}

    /* Write header if printing star entries as found */
    if (nstarmax < 1) {
	char *revmessage;
	revmessage = getrevmsg();
	title = CatName (ucat, refcatname);
	printf ("catalog        %s\n", title);
	free ((char *)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, ucat);
	printf ("%s	ra          	dec         	", catid);
	printf ("magb 	magr  	arcmin\n");
	printf ("-------------	------------	------------	");
	printf ("-----	-----	------\n");
	}

    uara1 = (int) (rra1 * 360000.0 + 0.5);
    uara2 = (int) (rra2 * 360000.0 + 0.5);
    uadec1 = (int) ((rdec1 * 360000.0) + 32400000.5);
    uadec2 = (int) ((rdec2 * 360000.0) + 32400000.5);
    
    if (nstarmax > ndist) {
	if (ndist > 0)
	    free ((void *)udist);
	udist = (double *) malloc (nstarmax * sizeof (double));
	if (udist == NULL) {
	    fprintf (stderr,"UACREAD:  cannot allocate separation array\n");
	    return (0);
	    }
	ndist = nstarmax;
	}

    /* Loop through region list */
    nstar = 0;
    for (iz = 0; iz < nz; iz++) {

    /* Get path to zone catalog */
	znum = zlist[iz];
	if ((nstars = uacopen (znum)) != 0) {

	    jstar = 0;
	    jtable = 0;
	    for (iwrap = 0; iwrap <= wrap; iwrap++) {

	    /* Find first star based on RA */
		if (iwrap == 0 || wrap == 0)
		    istar1 = uacsra (rra1);
		else
		    istar1 = 1;

	    /* Find last star based on RA */
		if (iwrap == 1 || wrap == 0)
		    istar2 = uacsra (rra2);
		else
		    istar2 = nstars;

		if (istar1 == 0 || istar2 == 0)
		    break;

		nread = istar2 - istar1 + 1;
		itable = 0;

	    /* Loop through zone catalog for this region */
		for (istar = istar1; istar <= istar2; istar++) {
		    itable ++;
		    jtable ++;

		    if (uacstar (istar, &star)) {
			fprintf (stderr,"UACREAD: Cannot read star %d\n", istar);
			break;
			}

		/* Extract selected fields */
		    else {

		    /* Check position limits */
     			if ((star.decsec >= uadec1 && star.decsec <= uadec2) &&
			    ((wrap && (star.rasec>=uara1 || star.rasec<=uara2)) ||
			     (!wrap && (star.rasec>=uara1 && star.rasec<=uara2))
			    )){

			/* Check magnitude, distance, and plate number */
			    magb = uacmagb (star.magetc);
			    magr = uacmagr (star.magetc);
			    if (magsort == 1)
				mag = magr;
			    else
				mag = magb;

			    /* Check magnitude limits */
			    pass = 1;
			    if (mag1 != mag2 && (mag < mag1 || mag > mag2))
				pass = 0;

			    /* Check plate ID */
			    plate = uacplate (star.magetc);
			    if (xplate != 0 && plate != xplate)
				pass = 0;

			    /* Check position limits */
			    if (pass) {
				ra = uacra (star.rasec);
				dec = uacdec (star.decsec);
				wcscon (sysref,sysout,eqref,eqout,&ra,&dec,epout);

				/* Compute distance from search center */
				if (distsort || drad > 0)
				    dist = wcsdist (cra,cdec,ra,dec);
				else
				    dist = 0.0;
			    
				/* Check radial distance to search center */
				if (drad > 0.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) {
				num = (double) znum +
				      (0.00000001 * (double)istar);

			    /* Write star position and magnitudes to stdout */
				if (nstarmax < 1) {
				    CatNum (ucat, -13, 0, num, 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	%.2f	%.2f\n",
					magb, magr, dist / 60.0);
				    }

			    /* Save star position and magnitude in table */
				else if (nstar < nstarmax) {
				    unum[nstar] = num;
				    ura[nstar] = ra;
				    udec[nstar] = dec;
				    umag[0][nstar] = magb;
				    umag[1][nstar] = magr;
				    uplate[nstar] = plate;
				    udist[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) {
					unum[farstar] = num;
					ura[farstar] = ra;
					udec[farstar] = dec;
					umag[0][farstar] = magb;
					umag[1][farstar] = magr;
					uplate[farstar] = plate;
					udist[farstar] = dist;

				    /* Find new farthest star */
					maxdist = 0.0;
					for (i = 0; i < nstarmax; i++) {
					    if (udist[i] > maxdist) {
						maxdist = udist[i];
						farstar = i;
						}
					    }
					}
				    }

			    /* If too many stars, replace faintest star */
				else if (mag < faintmag) {
				    unum[faintstar] = num;
				    ura[faintstar] = ra;
				    udec[faintstar] = dec;
				    umag[0][faintstar] = magb;
				    umag[1][faintstar] = magr;
				    uplate[faintstar] = plate;
				    udist[faintstar] = dist;

			    /* Find new faintest star */
				    faintmag = 0.0;
				    for (i = 0; i < nstarmax; i++) {
					if (umag[magsort][i] > faintmag) {
					    faintmag = umag[magsort][i];
					    faintstar = i;
					    }
					}
				    }
				nstar++;
				jstar++;
				if (nlog == 1)
				    fprintf (stderr,"UACREAD: %04d.%08d: %9.5f %9.5f %s %5.2f %5.2f\n",
					znum,istar,ra,dec,cstr,magb,magr);

			    /* End of accepted star processing */
				}
			    }

		    /* End of individual star processing */
			}

		/* Log operation */
		    if (nlog > 0 && itable%nlog == 0)
			fprintf (stderr,"UACREAD: zone %d (%2d / %2d) %8d / %8d / %8d sources\r",
				znum, iz+1, nz, jstar, itable, nread);

		/* End of star loop */
		    }

		/* End of wrap loop */
		}

	/* Close zone input file */
	    (void) fclose (fcat);
	    itot = itot + itable;
	    if (nlog > 0)
		fprintf (stderr,"UACREAD: zone %d (%2d / %2d) %8d / %8d / %8d sources      \n",
			znum, iz+1, nz, jstar, jtable, nstars);

	/* End of zone processing */
	    }

    /* End of zone loop */
	}

/* Summarize search */
    if (nlog > 0) {
	if (nz > 1)
	    fprintf (stderr,"UACREAD: %d zones: %d / %d found\n",nz,nstar,itot);
	else
	    fprintf (stderr,"UACREAD: 1 zone: %d / %d found\n",nstar,itable);
	if (nstar > nstarmax)
	    fprintf (stderr,"UACREAD: %d stars found; only %d returned\n",
		     nstar,nstarmax);
	}
    return (nstar);
}


int
usarnum (nnum,sysout,eqout,epout,unum,ura,udec,umag,uplate,nlog)

int	nnum;		/* 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	*unum;		/* Array of UA numbers to find */
double	*ura;		/* Array of right ascensions (returned) */
double	*udec;		/* Array of declinations (returned) */
double	**umag;		/* Array of blue and red magnitudes (returned) */
int	*uplate;	/* Array of plate numbers (returned) */
int	nlog;		/* Logging interval */
{
    int i;
    int uacrnum();

    ucat = USA1;
    i = uacrnum ("USAC",nnum,sysout,eqout,epout,unum,ura,udec,umag,uplate,nlog);
    ucat = UA1;
    return (i);
}


int
uacrnum (refcatname,nnum,sysout,eqout,epout,unum,ura,udec,umag,uplate,nlog)

char	*refcatname;	/* Name of catalog (UAC, USAC, UAC2, USAC2) */
int	nnum;		/* 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	*unum;		/* Array of UA numbers to find */
double	*ura;		/* Array of right ascensions (returned) */
double	*udec;		/* Array of declinations (returned) */
double	**umag;		/* Array of blue and red magnitudes (returned) */
int	*uplate;	/* Array of plate numbers (returned) */
int	nlog;		/* Logging interval */
{
    UACstar star;	/* UA catalog entry for one star */
    int sysref=WCS_J2000;	/* Catalog coordinate system */
    double eqref=2000.0;	/* Catalog equinox */
/*    int isp;
    char ispc[2]; */

    int znum;
    int jnum;
    int nzone;
    int nfound = 0;
    double ra,dec;
    double magr, magb;
    double dstar;
    int istar, plate;
    char *str;

    /* Set catalog code and path to catalog */
    if (strncmp (refcatname,"us",2)==0 ||
        strncmp (refcatname,"US",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("USA2_PATH")) != NULL)
		strcpy (usa2path,str);
	    ucat = USA2;
	    uapath = usa2path;
	    }
	else {
	    if ((str = getenv("USA1_PATH")) != NULL)
		strcpy (usa1path,str);
	    ucat = USA1;
	    uapath = usa1path;
	    }
	}
    else if (strncmp (refcatname,"ua",2)==0 ||
        strncmp (refcatname,"UA",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("UA2_PATH")) != NULL)
		strcpy (ua2path,str);
	    else if ((str = getenv("UA2_ROOT")) != NULL) {
		ua2path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA2;
	    uapath = ua2path;
	    }
	else {
	    if ((str = getenv("UA1_PATH")) != NULL)
		strcpy (ua1path,str);
	    else if ((str = getenv("UA1_ROOT")) != NULL) {
		ua1path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA1;
	    uapath = ua1path;
	    }
	}
    else {
	fprintf (stderr, "UACREAD:  %s not a USNO catalog\n", refcatname);
	return (0);
	}

    /* If root pathname is a URL, search and return */
    if (!strncmp (uapath, "http:",5)) {
	return (webrnum (uapath,refcatname,nnum,sysout,eqout,epout,1,
			 unum,ura,udec,NULL,NULL,umag,uplate,nlog));
	}


/* Loop through star list */
    for (jnum = 0; jnum < nnum; jnum++) {

    /* Get path to zone catalog */
	znum = (int) unum[jnum];
	if ((nzone = uacopen (znum)) != 0) {
	    dstar = (unum[jnum] - znum) * 100000000.0;
	    istar = (int) (dstar + 0.5);
	    if (istar > nzone) {
		fprintf (stderr,"UACRNUM: Star %d > max. in zone %d\n",
			 istar,nzone);
		break;
		}

	    if (uacstar (istar, &star)) {
		fprintf (stderr,"UACRNUM: Cannot read star %d\n", istar);
		break;
		}

	    /* Extract selected fields */
	    else {
		ra = uacra (star.rasec); /* Right ascension in degrees */
		dec = uacdec (star.decsec); /* Declination in degrees */
		magb = uacmagb (star.magetc); /* Blue magnitude */
		magr = uacmagr (star.magetc); /* Red magnitude */
		plate = uacplate (star.magetc);	/* Plate number */
		wcscon (sysref, sysout, eqref, eqout, &ra, &dec, epout);

		/* Save star position and magnitude in table */
		ura[nfound] = ra;
		udec[nfound] = dec;
		umag[0][nfound] = magb;
		umag[1][nfound] = magr;
		/* br2sp (NULL, magb, magr, ispc);
		isp = (1000 * (int)ispc[0]) + (int)ispc[1]; */
		uplate[nfound] = plate;

		nfound++;
		if (nlog == 1)
		    fprintf (stderr,"UACRNUM: %04d.%08d: %9.5f %9.5f %5.2f %5.2f\n",
			     znum,istar,ra,dec,magb,magr);

		/* Log operation */
		if (nlog > 0 && jnum%nlog == 0)
		    fprintf (stderr,"UACRNUM: %4d.%8d  %8d / %8d sources\r",
			     znum, istar, jnum, nnum);

		(void) fclose (fcat);
		/* End of star processing */
		}

	    /* End of star */
	    }

	/* End of star loop */
	}

    /* Summarize search */
    if (nlog > 0)
	fprintf (stderr,"UACRNUM:  %d / %d found\n",nfound,nnum);

    return (nfound);
}


/* USACBIN -- Fill a FITS WCS image with USNO SA Catalog stars */

int
usabin (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;		/* Logging interval */
{
    int i;
    int uacbin();

    ucat = USA2;
    i = uacbin ("usac", wcs, header, image, mag1, mag2, sortmag, magscale, nlog);
    ucat = USA2;
    return (i);
}


/* UACBIN -- Fill a FITS WCS image with USNO A or SA Catalog stars */

int
uacbin (refcatname, wcs, header, image, mag1, mag2, sortmag, magscale, nlog)

char	*refcatname;	/* Name of catalog (UAC, USAC, UAC2, USAC2) */
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 or 2) */
double	magscale;	/* Scaling factor for magnitude to pixel flux
			 * (number of catalog objects per bin if 0) */
int	nlog;		/* Logging interval */
{
    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 nz;		/* Number of input UA zone files */
    int zlist[NZONES];	/* List of input UA zones */
    UACstar star;	/* UA catalog entry for one star */
    int sysref=WCS_J2000;	/* Catalog coordinate system */
    double eqref=2000.0;	/* Catalog equinox */
    double epref=2000.0;	/* Catalog epoch */

    double rra1, rra2, rdec1, rdec2;
    int wrap, iwrap;
    int verbose;
    int znum, itot,iz;
    int jtable,jstar;
    int itable = 0;
    int nstar, nread;
    int uara1, uara2, uadec1, uadec2;
    double ra,dec, rdist, ddist;
    double mag, magb, magr;
    int istar, istar1, istar2, plate;
    int nzmax = NZONES;	/* Maximum number of declination zones */
/*    int isp;
    char ispc[2]; */
    int pass;
    int ix, iy;
    int magsort;
    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);

    itot = 0;
    if (nlog > 0)
	verbose = 1;
    else
	verbose = 0;

    /* 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 code and path to catalog */
    if (strncmp (refcatname,"us",2)==0 ||
        strncmp (refcatname,"US",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("USA2_PATH")) != NULL)
		strcpy (usa2path,str);
	    ucat = USA2;
	    uapath = usa2path;
	    }
	else {
	    if ((str = getenv("USA1_PATH")) != NULL)
		strcpy (usa1path,str);
	    ucat = USA1;
	    uapath = usa1path;
	    }
	}
    else if (strncmp (refcatname,"ua",2)==0 ||
        strncmp (refcatname,"UA",2)==0) {
	if (strchr (refcatname, '2') != NULL) {
	    if ((str = getenv("UA2_PATH")) != NULL)
		strcpy (ua2path,str);
	    else if ((str = getenv("UA2_ROOT")) != NULL) {
		ua2path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA2;
	    uapath = ua2path;
	    }
	else {
	    if ((str = getenv("UA1_PATH")) != NULL)
		strcpy (ua1path,str);
	    else if ((str = getenv("UA1_ROOT")) != NULL) {
		ua1path[0] = 0;
		strcpy (cdroot,str);
		}
	    ucat = UA1;
	    uapath = ua1path;
	    }
	}
    else {
	fprintf (stderr, "UACBIN:  %s not a USNO catalog\n", refcatname);
	return (0);
	}

    /* 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);

    /* mag1 is always the smallest magnitude */
    if (mag2 < mag1) {
	mag = mag2;
	mag2 = mag1;
	mag1 = mag;
	}
    if (sortmag == 1)
	magsort = 0;
    else
	magsort = 1;

    /* Find UA Star Catalog regions in which to search */
    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);
    nz = uaczones (rra1, rra2, rdec1, rdec2, nzmax, zlist, verbose);
    if (nz <= 0) {
	fprintf (stderr, "UACBIN:  no USNO A zones found\n");
	return (0);
	}

    uara1 = (int) (rra1 * 360000.0 + 0.5);
    uara2 = (int) (rra2 * 360000.0 + 0.5);
    uadec1 = (int) ((rdec1 * 360000.0) + 32400000.5);
    uadec2 = (int) ((rdec2 * 360000.0) + 32400000.5);
    
    /* Loop through region list */
    nstar = 0;
    for (iz = 0; iz < nz; iz++) {

    /* Get path to zone catalog */
	znum = zlist[iz];
	if ((nstars = uacopen (znum)) != 0) {

	    jstar = 0;
	    jtable = 0;
	    for (iwrap = 0; iwrap <= wrap; iwrap++) {

	    /* Find first star based on RA */
		if (iwrap == 0 || wrap == 0)
		    istar1 = uacsra (rra1);
		else
		    istar1 = 1;

	    /* Find last star based on RA */
		if (iwrap == 1 || wrap == 0)
		    istar2 = uacsra (rra2);
		else
		    istar2 = nstars;

		if (istar1 == 0 || istar2 == 0)
		    break;

		nread = istar2 - istar1 + 1;
		itable = 0;

	    /* Loop through zone catalog for this region */
		for (istar = istar1; istar <= istar2; istar++) {
		    itable ++;
		    jtable ++;

		    if (uacstar (istar, &star)) {
			fprintf (stderr,"UACBIN: Cannot read star %d\n", istar);
			break;
			}

		/* Extract selected fields */
		    else {

		    /* Check position limits */
     			if ((star.decsec >= uadec1 && star.decsec <= uadec2) &&
			    ((wrap && (star.rasec>=uara1 || star.rasec<=uara2)) ||
			     (!wrap && (star.rasec>=uara1 && star.rasec<=uara2))
			    )){

			/* Check magnitude, distance, and plate number */
			    magb = uacmagb (star.magetc);
			    magr = uacmagr (star.magetc);
			    if (magsort == 1)
				mag = magr;
			    else
				mag = magb;

			    /* Check magnitude limits */
			    pass = 1;
			    if (mag1 != mag2 && (mag < mag1 || mag > mag2))
				pass = 0;

			    /* Check plate ID */
			    plate = uacplate (star.magetc);
			    if (xplate != 0 && plate != xplate)
				pass = 0;

			    /* Check position limits */
			    if (pass) {
				ra = uacra (star.rasec);
				dec = uacdec (star.decsec);
				wcscon (sysref,sysout,eqref,eqout,&ra,&dec,epout);

				/* 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,"UACBIN: %04d.%08d: %9.5f %9.5f %s",
					     znum,istar,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 */
				}
			    }

		    /* End of individual star processing */
			}

		/* Log operation */
		    if (nlog > 0 && itable%nlog == 0)
			fprintf (stderr,"UACBIN: zone %d (%2d / %2d) %8d / %8d / %8d sources\r",
				znum, iz+1, nz, jstar, itable, nread);

		/* End of star loop */
		    }

		/* End of wrap loop */
		}

	/* Close zone input file */
	    (void) fclose (fcat);
	    itot = itot + itable;
	    if (nlog > 0)
		fprintf (stderr,"UACBIN: zone %d (%2d / %2d) %8d / %8d / %8d sources      \n",
			znum, iz+1, nz, jstar, jtable, nstars);

	/* End of zone processing */
	    }

    /* End of zone loop */
	}

/* Summarize search */
    if (nlog > 0) {
	if (nz > 1)
	    fprintf (stderr,"UACBIN: %d zones: %d / %d found\n",nz,nstar,itot);
	else
	    fprintf (stderr,"UACBIN: 1 zone: %d / %d found\n",nstar,itable);
	}
    return (nstar);
}


/* Declination zone numbers */
int azone[NZONES]={0,75,150,225,300,375,450,525,600,675,750,825,900,
	      975,1050,1125,1200,1275,1350,1425,1500,1575,1650,1725};

/* UACZONES -- figure out which UA zones will need to be searched */

static int
uaczones (ra1, ra2, dec1, dec2, nzmax, zones, verbose)

double	ra1, ra2;	/* Right ascension limits in degrees */
double	dec1, dec2; 	/* Declination limits in degrees */
int	nzmax;		/* Maximum number of zones to find */
int	*zones;		/* Region numbers (returned)*/
int	verbose;	/* 1 for diagnostics */

{
    int nrgn;		/* Number of zones found (returned) */
    int iz,iz1,iz2,i;

    for (i = 0; i < nzmax; i++)
	zones[i] = 0;

    nrgn = 0;

/* Find zone range to search based on declination */
    iz1 = uaczone (dec1);
    iz2 = uaczone (dec2);

/* Tabulate zones to search */
    i = 0;
    if (iz2 >= iz1) {
	for (iz = iz1; iz <= iz2; iz++)
	    zones[i++] = azone[iz];
	}
    else {
	for (iz = iz2; iz <= iz1; iz++)
	    zones[i++] = azone[iz];
	}

    nrgn = i;
    if (verbose) {
	fprintf(stderr,"UACZONES:  %d zones: %d - %d\n",nrgn,zones[0],zones[i-1]);
	fprintf(stderr,"UACZONES: RA: %.5f - %.5f, Dec: %.5f - %.5f\n",ra1,ra2,dec1,dec2);
	}

    return (nrgn);
}


/* UACRA -- returns right ascension in degrees from the UA star structure */

static double
uacra (rasec)

int rasec;	/* RA in 100ths of arcseconds from UA catalog entry */
{
    return ((double) (rasec) / 360000.0);
}


/* UACDEC -- returns the declination in degrees from the UA star structure */

static double
uacdec (decsec)

int decsec;	/* Declination in 100ths of arcseconds from UA catalog entry */
{
    return ((double) (decsec - 32400000) / 360000.0);
}


/* UACMAGR -- returns the red magnitude from the UA star structure */

static double
uacmagr (magetc)

int magetc;	/* Quality, plate, and magnitude from UA catalog entry */
{
    if (magetc < 0)
	return ((double) (-magetc % 1000) * 0.1);
    else
	return ((double) (magetc % 1000) * 0.1);
}


/* UACMAGB -- returns the blue magnitude from the UA star structure */

static double
uacmagb (magetc)

int magetc;	/* Quality, plate, and magnitude from UA catalog entry */
{
    if (magetc < 0)
	return ((double) ((-magetc / 1000) % 1000) * 0.1);
    else
	return ((double) ((magetc / 1000) % 1000) * 0.1);
}


/* UACMAGERR -- returns 1 if magnitude is uncertain from UA star structure */

int
uacmagerr (magetc)

int magetc;	/* Quality, plate, and magnitude from UA catalog entry */
{
    if (magetc < 0)
	return ((-magetc / 1000000000) % 10);
    else
	return ((magetc / 1000000000) % 10);
}


/* UACGSC -- returns 1 if UA star is in the HST Guide Star Catalog */

int
uacgsc (magetc)

int magetc;	/* Quality, plate, and magnitude from UA catalog entry */
{
    if (magetc < 0)
	return (1);
    else
	return (0);
}


/* UACPLATE -- returns the plate number from the UA star structure */

static int
uacplate (magetc)

int magetc;	/* Quality, plate, and magnitude from UA catalog entry */
{
    if (magetc < 0)
	return ((-magetc / 1000000) % 1000);
    else
	return ((magetc / 1000000) % 1000);
}


/* UACZONE -- find the UA zone number where a declination can be found */

static int
uaczone (dec)

double dec;	/* declination in degrees */
{
    double zonesize = 7.5;	/* number of degrees per declination zone */
    int zone;

    zone = (int) ((dec + 90.0) / zonesize);
    if (zone > 23)
	zone = 23;
    else if (zone < 0)
	zone = 0;
    return (zone);
}


/* UACSRA -- Find UA star closest to specified right ascension */

static int
uacsra (rax0)

double	rax0;		/* Right ascension in degrees for which to search */
{
    int istar, istar1, istar2, nrep;
    double rax, ra1, ra, rdiff, rdiff1, rdiff2, sdiff;
    UACstar star;	/* UA catalog entry for one star */
    char rastrx[32];
    int debug = 0;

    rax = rax0;
    if (debug)
	ra2str (rastrx, 31, rax, 3);
    istar1 = 1;
    if (uacstar (istar1, &star))
	return (0);
    ra1 = uacra (star.rasec);
    istar = nstars;
    nrep = 0;
    while (istar != istar1 && nrep < 20) {
	if (uacstar (istar, &star))
	    break;
	else {
	    ra = uacra (star.rasec);
	    if (ra == ra1)
		break;
	    if (debug) {
		char rastr[32];
		ra2str (rastr, 31, ra, 3);
		fprintf (stderr,"UACSRA %d %d: %s (%s)\n",
			 nrep,istar,rastr,rastrx);
		}
	    rdiff = ra1 - ra;
	    rdiff1 = ra1 - rax;
	    rdiff2 = ra - rax;
	    if (nrep > 20 && ABS(rdiff2) > ABS(rdiff1)) {
		istar = istar1;
		break;
		}
	    nrep++;
	    sdiff = (double)(istar - istar1) * rdiff1 / rdiff;
	    istar2 = istar1 + (int) (sdiff + 0.5);
	    ra1 = ra;
	    istar1 = istar;
	    istar = istar2;
	    if (debug) {
		fprintf (stderr," ra1=    %.5f ra=     %.5f rax=    %.5f\n",
			 ra1,ra,rax);
		fprintf (stderr," rdiff=  %.5f rdiff1= %.5f rdiff2= %.5f\n",
			 rdiff,rdiff1,rdiff2);
		fprintf (stderr," istar1= %d istar= %d istar1= %d\n",
			 istar1,istar,istar2);
		}
	    if (istar < 1)
		istar = 1;
	    if (istar > nstars)
		istar = nstars;
	    if (istar == istar1)
		break;
	    }
	}
    return (istar);
}

/* UACOPEN -- Open UA Catalog zone catalog, returning number of entries */

static int
uacopen (znum)

int znum;	/* UA Catalog zone */
{
    char zonepath[64];	/* Pathname for input UA zone file */
    UACstar star;	/* UA catalog entry for one star */
    int lfile;
    
/* Get path to zone catalog */
    if (uacpath (znum, zonepath)) {
	fprintf (stderr, "UACOPEN: Cannot find zone catalog for %d\n", znum);
	return (0);
	}

/* Find number of stars in zone catalog by its length */
    lfile = getfilesize (zonepath);
    if (lfile < 2) {
	fprintf (stderr,"UA zone catalog %s has no entries\n",zonepath);
	return (0);
	}
    else
	nstars = lfile / 12;

/* Open zone catalog */
    if (!(fcat = fopen (zonepath, "rb"))) {
	fprintf (stderr,"UA zone catalog %s cannot be read\n",zonepath);
	return (0);
	}

/* Check to see if byte-swapping is necessary */
    cswap = 0;
    if (uacstar (1, &star)) {
	fprintf (stderr,"UACOPEN: cannot read star 1 from UA zone catalog %s\n",
		 zonepath);
	return (0);
	}
    else {
	if (star.rasec > 360 * 360000 || star.rasec < 0) {
	    cswap = 1;
	    /* fprintf (stderr,"UACOPEN: swapping bytes in UA zone catalog %s\n",
		     zonepath); */
	    }
	else if (star.decsec > 180 * 360000 || star.decsec < 0) {
	    cswap = 1;
	    /* fprintf (stderr,"UACOPEN: swapping bytes in UA zone catalog %s\n",
		     zonepath); */
	    }
	else
	    cswap = 0;
	}

    return (nstars);
}


/* UACPATH -- Get UA Catalog region file pathname */

static int
uacpath (zn, path)

int zn;		/* UA zone number */
char *path;	/* Pathname of UA zone file */

{
    int iz;		/* Zone index (0000 = 0, 0075 = 1, ...) */
    int icd;		/* CDROM number if multiple CDROMs used */

    /* Return error code and null path if zone is out of range */
    if (zn < 0 || zn > 1725) {
	fprintf (stderr, "UACPATH: zone %d out of range 0-1725\n",zn);
	path[0] = 0;
	return (-1);
	}

    /* Set path for USNO SA zone catalog */
    if (ucat == USA1 || ucat == USA2)
	sprintf (path,"%s/zone%04d.cat", uapath, zn);

    /* Set zone catalog path when USNO A is in a single directory */
    else if (strlen (uapath) > 0)
	sprintf (path,"%s/zone%04d.cat", uapath, zn);

    /* Set zone catalog path when USNO A is read from CDROMs */
    else {
	iz = zn / 75;
	if (ucat == UA1)
	    icd = zdisk1[iz];
	else
	    icd = zdisk2[iz];
	sprintf (path,"%s/%s/zone%04d.cat", cdroot, cdname[icd-1], zn);
	}

    return (0);
}


/* UACSTAR -- Get UA catalog entry for one star; return 0 if successful */

static int
uacstar (istar, star)

int istar;	/* Star sequence number in UA zone catalog */
UACstar *star;	/* UA catalog entry for one star */
{
    int nbs, nbr, nbskip;

    if (istar < 1 || istar > nstars) {
	fprintf (stderr, "UACstar %d is not in catalog\n",istar);
	return (-1);
	}
    nbskip = 12 * (istar - 1);
    if (fseek (fcat,nbskip,SEEK_SET))
	return (-1);
    nbs = sizeof (UACstar);
    nbr = fread (star, nbs, 1, fcat) * nbs;
    if (nbr < nbs) {
	fprintf (stderr, "UACstar %d / %d bytes read\n",nbr, nbs);
	return (-2);
	}
    if (cswap)
	uacswap ((char *)star);
    return (0);
}


/* UACSWAP -- Reverse bytes of UA Catalog entry */

static void
uacswap (string)

char *string;	/* Start of vector of 4-byte ints */

{
char *sbyte, *slast;
char temp0, temp1, temp2, temp3;
int nbytes = 12; /* Number of bytes to reverse */

    slast = string + nbytes;
    sbyte = string;
    while (sbyte < slast) {
	temp3 = sbyte[0];
	temp2 = sbyte[1];
	temp1 = sbyte[2];
	temp0 = sbyte[3];
	sbyte[0] = temp0;
	sbyte[1] = temp1;
	sbyte[2] = temp2;
	sbyte[3] = temp3;
	sbyte = sbyte + 4;
	}
    return;
}

/* Nov 15 1996	New subroutine
 * Dec 11 1996	Set ra<0 to ra+360 and ra>360 to ra-360
 * Dec 16 1996	Add code to read a specific star
 * Dec 17 1996	Keep closest stars, not brightest, if searching within radius
 *
 * Mar 12 1997	Set paths for SA-1.0 and multiple CDROM A-1.0
 * Mar 20 1997	Clean up UACRNUM after lint
 * Apr 23 1997	Fix bug which rejected stars in HST Guide Star Catalog
 * Nov  6 1997	Don't print star overrun unless logging
 *
 * Feb 20 1998	Speed up processing by searching in arcseconds, not degrees
 * Feb 20 1998	Speed up processing by searching RA and Dec, then rest
 * Apr 20 1998	Fix bug so stars within radius can be found
 * Jun 24 1998	Add string lengths to ra2str() and dec2str() calls
 * Jun 24 1998	Initialize byte-swapping flag in UACOPEN()
 * Sep 15 1998	Fix bug setting A 1.0 region path on CDROM found by Naoki Yasuda
 * Sep 22 1998	Convert coordinates in these subroutines
 * Oct  8 1998	Fix bug in uacread call in usaread() and uacrnum call in usarnum()
 * Oct 26 1998	Fix bug in search algorith for non-J2000 searches
 * Oct 29 1998	Correctly assign numbers when too many stars are found
 * Nov 20 1998	Add support for USNO A-2.0 and SA-2.0 catalogs
 * Nov 24 1998	Fix bug reading SA-2.0 catalog
 *
 * Feb  9 1999	Improve documentation
 * 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 uacread()
 * Sep 10 1999	Set plate 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
 * Oct 20 1999	Include wcscat.h
 * Oct 21 1999	Clean up code after lint
 *
 * Jun  9 2000	Fix bug detecting swapped files on Alphas and PCs if RA=0
 * Jun 26 2000	Add coordinate system to SearchLim() arguments
 * Sep 22 2000	Return approximate spectral type instead of plate number
 * Nov 28 2000	Add option to read catalog using HTTP
 * Dec  1 2000	Return plate number, not bad spectral type
 * Dec 11 2000	Path may be set to catalog search engine URL
 * Dec 15 2000	Do away with separate usapath variable; use uapath
 *
 * Jun 14 2001	Make sure star number is correct in uacrnum()
 * Jun 27 2001	Use RefCat codes for ucat
 * Jun 27 2001	Print stars as found in uacread() if nstarmax < 1
 * Jun 27 2001	Allocate udist only when larger array is needed
 * Sep 11 2001	Change to single magnitude argeument
 * Sep 11 2001	Add sort magnitude argument to uacread()
 * Sep 19 2001	Drop fitshead.h; it is in wcs.h
 * Sep 21 2001	Add commented-out URL of ESO web catalog server
 * Nov 20 2001	Change cos(degrad()) to cosdeg()
 *
 * Apr 10 2002	Simplify use of magsort
 * Jul 31 2002	Drop extra magb argument in uacrnum()
 * Oct  2 2002	Print current scat revision message
 *
 * Jan 21 2003	Print arcminute radial distance, not arcsecond for instant out
 * Feb  4 2003	Open catalog file rb instead of r (Martin Ploner, Bern)
 * Mar 10 2003	Improve test for position
 * Apr  3 2003	Drop unused variables after lint
 * Apr 14 2003	Explicitly get revision date if nstarmax < 1
 * May 27 2003	Use getfilesize() to get file size
 * Aug 22 2003	Add radi argument for inner edge of search annulus
 * Sep 25 2003	Add usabin() and uacbin() to fill an image with sources
 * Oct  6 2003	Update uacread() and uacbin() for improved RefLim()
 * Nov 18 2003	Initialize image size and bits/pixel from header in uacbin()
 * Dec  1 2003	Add missing tab to n=-1 header
 * Dec 12 2003	Fix bug in wcs2pix() call in uacbin(); fix usacbin() subroutine
 *
 * Aug 30 2004	Include fitsfile.h and math.h
 *
 * Sep 13 2006	Change default paths to /data/astrocat
 * Sep 26 2006	Increase length of rastr and destr from 16 to 32
 * Nov 15 2006	Fix binning
 *
 * Jan 10 2007	Add match=1 argument to webrnum()
 * Jan 10 2007	Add dradi arguemnt to uacread() call in usacread()
 */