/*++++++++++++++
.IDENTIFICATION ucac1.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-UCAC1 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   18-May-2001
.VERSION  1.1   05-Dec-2002: Adapted message for truncated results
.VERSION  1.2   13-Jan-2004: Errors around 0h ?
.COMMENTS       Access to UCAC catalogue
	Remember that the environment variable UCACroot provides
	the root of UCAC compressed catalogue; its default is
	/UCAC1
---------------*/

#include <ucac.h>	/* Structure definitions */

#define RADIUS		(7.5/60.)	/* Default target radius */

#include <stdio.h>
#include <stdlib.h>	/* Malloc   */
#include <string.h>
#include <ctype.h>
#include <math.h>

#define ITEMS(a)	(sizeof(a)/sizeof(a[0]))
#define MIN(a,b)	((a)<(b) ? (a) : (b))
#define MAX(a,b)	((a)>(b) ? (a) : (b))
#define issign(c)	(((c)=='+')||((c)=='-'))

#define  sind(x)	sin(x*(M_PI/180.))
#define  cosd(x)	cos(x*(M_PI/180.))
#define  tand(x)	tan(x*(M_PI/180.))
#define asind(x)	(180./M_PI)*asin(x)
#define atand(x)	(180./M_PI)*atan(x)

#define adeg(x)		((x)*3600*1000)	/* Degree into mas */

/* Options for UCAC queries.......	*/
int ucac_options = 0 ;

/* Linked records, to sort them...	*/
typedef struct linked_UCAC {
    struct linked_UCAC *prev ;
    UCACrec rec ;
} LUCAC ;
static int mrec = 100 ;
static int irec, truncated ;
static int input_id ;		/* ID(1) Zone (2)  */
static char *prompt[] = {
   "Center", 	/* input_id=0 */
   "UCAC-ID", 	/* input_id=1 */
   "UCAC Zone",	/* input_id=2 */
} ;
static char whole_ucac ;	/* -whole option   */
LUCAC *last, *arec, *just_added;

/*  Definitions related to Center of Field */
static double radius = 0 ;
static double boxy[2] ;
static int (*digest_routine)() ;
static double localR[3][3] ;	/* Matrix, [0] = dir.cos. of center */
static double du2max = -1 ;

/* Other options */
static int opted = 6 ; 	/* Default edit Decimal, Epoch, ID */

/* Definitions of the record fields */
struct s_fields {
    char letter ;
    char selected ;
    short order ;
    double factor ;
    long lim[2] ;
};
static struct s_fields thefields[] =  {
    {'r',0,0,3.6e6},			/* MUST BE #0 */
    {'i',0,0,1.},			/* MUST BE #1 */
    {'x',0,0,(180./M_PI)*3.6e6},	/* MUST BE #2 */
    {'y',0,0,(180./M_PI)*3.6e6},	/* MUST BE #3 */
    {'a',0,0,3.6e6},			/* Must be #4 */
    {'d',0,0,3.6e6},			/* Must be #5 */
    {'m',0,0,100.},	/* Mag.	*/
    {'e',0,0,1000.},	/* Ep	*/
    {'o',0,0,1.},	/* Nobs	*/
    {'c',0,0,1.},	/* Ncat	*/
    {'p',0,0,10.},	/* pmt	*/
    {0}
} ;
static struct s_fields *compare_fields[16] ;
static struct s_fields *check_fields[16] ;
static double _factor_  = 1 ;
static long matched ;	/* Matching records 	      	*/
static int  stopid;	/* 1 = Stop after last ID read 	*/

static char usage[] = "\
Usage: ucac [-HELP|-v] [-R root_name] [-r[s] [min,]radius] [-b[s] x[,y]]\n\
            [-e edit_opt] [-f input_file] [-m max_records]\n\
	    [-c center | -i ID | -z zone] [-l! min,max] [-s !] [-whole]\n\
";
 
static char help[] = "\
  -HELP: display column explanations\n\
     -v: verbose option (details what the program does)\n\
     -b: target box in arcmin ;    -bs = target box in arcsec\n\
     -r: target radius in arcmin ; -rs = target radius in arcsec\n\
     -c: target center in decimal or sexagesimal (default in stdin)\n\
     -e: s=edit Sexa; m=edit mas; i=edit UCAC-ID, e=edit ObsDate, x=edit x,y\n\
     -f: specifies an input file (default stdin)\n\
    -id: query from an UCAC1 (8-digit) number\n\
     -m: max number of stars to retrieve\n\
    -l!: Set the limits (range) on one of the parameters (below)\n\
    -s!: Sort the result by the parameter ! (list below)\n\
  -zxxc: search on a zone (file name between 001 and 169)\n\
 -whole: search on the whole sky\n\
     -R: Root (directory) name where the UCAC files are located ($UCACroot)\n\
====The abbreviations of the parameters (symbolized !)are:\n\
      a=Alpha   c=Ncats    d=Delta     e=Epoch    i=UCAC1    m=mag   o=Nobs\n\
      p=pmotion r=distance s=sigmaPos  x=proj.E   y=proj.N   z=Zone\n\
";

static char HELP[] = "\
Access to UCAC1 (USNO CCD Astrograph Catalog) at CDS\n\
=============================================================================\n\
\n\
The fields have the following meaning:\n\
\n\
UCAC1     = Sequential number in original catalogue (in range 1-27425433)\n\
RA+Dec    = Position (Alpha, Delta) in ICRS (Hipparcos system) at Epoch\n\
            and their mean errors in mas (mean errors of RAcos(DE) and of DE)\n\
UCmag     = magnitude (bandpass 579-642nm, between V and R), accuracy 0.3mag\n\
No        = number of observations (images per star) used for the position\n\
Epoch     = Epoch (mean date, in Julian years) of the position\n\
pmRA,pmDE = Proper motions in mas/yr, Hipparcos system, and their mean errors\n\
Nc        = number of catalog positions used for proper motion determination\n\
r(\")     = distance from specified center, in arc seconds.\n\
=============================================================================\n\
For details concerning UCAC1, see  http://ad.usno.navy.mil/ucac/\n\
=============================================================================\n\
"; 

/*==================================================================
		Get the parameters
 *==================================================================*/

static struct s_fields *get_field(struct s_fields *f, char *p)
/*++++++++++++++++
.PURPOSE  Interpret the meaning of a parameter
.RETURNS  The relevant field / NULL
-----------------*/
{
    while (f->letter) {
	if (f->letter == *p) return(f) ;
	f++ ;
    }
    return((struct s_fields *)0) ;
}

/*============================================================================*/
static void tr_ou(double o[2], double u[3])
/*++++++++++++++++
.PURPOSE  Compute direction cosines
.RETURNS  ---
-----------------*/
{
    u[0] = u[1] = cosd(o[1]) ;
    u[0] *= cosd(o[0]) ;
    u[1] *= sind(o[0]) ;
    u[2] = sind(o[1]) ;
}

static void tr_oR ( o , R ) 
/*++++++++++++++++
.PURPOSE  Creates the rotation matrix R[3][3] defined as
	 R[0] (first row) = unit vector towards Zenith
	 R[1] (second row) = unit vector towards East
	 R[2] (third row) = unit vector towards North
.RETURNS ---
-----------------*/
  double o[2]; 		/* IN: original angles */
  double R[3][3];	/* OUT: rotation matrix */
{
    R[2][2] =  cosd(o[1]);
    R[0][2] =  sind(o[1]);
    R[1][1] =  cosd(o[0]);
    R[1][0] =  -sind(o[0]);
    R[1][2] =  0.e0;
    R[0][0] =  R[2][2] * R[1][1];  
    R[0][1] = -R[2][2] * R[1][0];
    R[2][0] = -R[0][2] * R[1][1];
    R[2][1] =  R[0][2] * R[1][0];
}

static void tr_uu( u1 , u2 , R ) 
/*++++++++++++++++
.PURPOSE  Rotates the unit vector u1 to u2, as  (u2) = (R) * (u1)
.RETURNS  ---
-----------------*/
  double u1[3]; 	/* IN: Unit vector */
  double u2[3]; 	/* OUT: Resulting unit vector after rotation */
  double R[3][3];	/* IN: rotation matrix (e.g. created by tr_oR)*/
{
  register int i;
  double u_stack[3];	/* allows same address for input/output      */

    for (i=0; i<3; i++) 
	u_stack[i] = R[i][0]*u1[0] +R[i][1]*u1[1] +  R[i][2]*u1[2] ;
    u2[0] = u_stack[0]; 	/* copies to output */
    u2[1] = u_stack[1]; 	/* copies to output */
    u2[2] = u_stack[2]; 	/* copies to output */
}

/*============================================================================*/

static int match_unum (char *string, double *value)
/*++++++++++++++++
.PURPOSE  Analyze the incoming and get unsigned number (no E notation)
.RETURNS  How many bytes matched
.REMARKS  Initial blanks return 0
-----------------*/
{
  char *p;
    for (p=string; isdigit(*p); p++) ;
    if (*p == '.') for (++p; isdigit(*p); p++) ;
    if (p == string) *value = 0 ;
    else *value = _factor_ * atof(string) ;
    return(p-string);
}

static int match_num (char *string, double *value)
/*++++++++++++++++
.PURPOSE  Analyze the incoming and get signed number (no E notation)
.RETURNS  How many bytes matched
.REMARKS  The internal _value_ is set
-----------------*/
{
  char *p, *s ;
    for (s=string; *s == ' '; s++) ;
    p = s ; 
    if ((*p == '+') || (*p == '-')) p++ ;
    p += match_unum(p, value) ;
    if (*s == '-') *value = -*value ;
    return(p-string);
}

static int get_lim (char *string, double values[2])
/*++++++++++++++++
.PURPOSE  Find out the limits (2 numbers separated by ,)
.RETURNS  0 / 1 / 2 / 3
.REMARKS  Accept p=min,max=min+max ==> sqrt(min1^2+min2^2) / 
-----------------*/
{
  int found, round, n ;
  double val ;
  char *p ;

    found = 0 ; p = string ; values[0] = values[1] = 0 ;
    round = 0 ;
  NextRound:
    round ^= 1;
    n = match_num(p, &val) ;
    if (n) {
        found |= 1, p += n ;
	if (round) values[0] = val ;
	else       values[0] = sqrt(values[0]*values[0] + val*val) ;
    }

    if (*p) {
	p++ ; 
	n = match_num(p, &val) ;
	if (n) {
            found |= 2, p += n ;
	}
	if (round) values[1] = val ;
	else       values[1] = sqrt(values[1]*values[1] + val*val) ;
    }
    if (*p && round) goto NextRound;
    return(found) ;
}

static int get_ucac (char *string, long id[2])
/*++++++++++++++++
.PURPOSE  Interpret a string for UCAC Numbers (range)
.RETURNS  Number of bytes matched
-----------------*/
{
  double v ; char *p ;

    p = string ;
    p += match_unum(p, &v) ;
    id[0] = v ;
    if (*p) {	/* Zone is specified */
	++p;
	p += match_unum(p, &v) ;
	id[1] = v ;
    }
    else id[1] = id[0] ;
    return(p-string) ;
}

static int get_center (char *string, double center[2])
/*++++++++++++++++
.PURPOSE  Interpret a string for RA + DEC
.RETURNS  Number of bytes matched
-----------------*/
{
  double v, factor ;
  char *p, *s ;

    factor = _factor_ ;
    center[0] = center[1] = 0 ;
    p = string ; while (*p == ' ') p++ ;
    p += match_unum(p, center) ;
    while (*p == ' ') p++ ;
    if (*p == ',') 	/* For those HEASARC guys ! */
	for (++p; *p == ' '; p++) ;
    if ((*p == '+') || (*p == '-')) ;
    else {		/* Sexagesimal */
	while ((*p == ':') || isdigit(*p)) {
	    _factor_ /= 60.; 
	    if (*p == ':') p++ ;
	    p += match_unum(p, &v);
	    while (*p == ' ') p++ ;
	    center[0] += v ;
	}
	center[0] *= 15. ;	/* Right Ascension */
    }
    s = p ;
    if ((*s == '+') || (*s == '-')) p++ ;
    _factor_ = 1 ;
    while ((*p == ':') || isdigit(*p)) {
	if (*p == ':') p++ ;
	p += match_unum(p, &v);
	while (*p == ' ') p++ ;
	center[1] += v ;
	_factor_ /= 60. ;
    }
    if (*s == '-') center[1] = -center[1] ;

    _factor_ = factor ; 	/* Restore */
    return(p-string) ;
}

/*==================================================================
		Check routine (return 0 = don't keep, 1 = keep)
 *==================================================================*/

static int check(UCACrec *rec)
{
  struct s_fields **f;
  double o[2], u[3], du, dr2 ;
  int i, flag ;

    /* Verify first the Position */
    flag = 0 ;
    rec->xy[0] = rec->xy[1] = NULLxy ;
    if (du2max >= 0) {
    	o[0] = rec->ra / 3.6e6 ;
    	o[1] = rec->sd / 3.6e6  - 90. ;
    	tr_ou(o, u); flag |= 1 ;
    	du = u[0] - localR[0][0] ; dr2  = du * du ;
    	du = u[1] - localR[0][1] ; dr2 += du * du ;
    	du = u[2] - localR[0][2] ; dr2 += du * du ;
    	if (dr2 > du2max) return(0) ;
    	rec->rho = thefields[0].factor * 2.*asind(0.5*sqrt(dr2)) ;
	if (opted&8) {	/* Compute x,y */
      	    tr_uu(u, u, localR) ;
      	    rec->xy[0] = thefields[2].factor*u[1]/u[0] ;
      	    rec->xy[1] = thefields[3].factor*u[2]/u[0] ;
      	    flag |= 2 ;
	}
    }
    else rec->rho = -1 ;

    /* printf("#id=%4d,%8d, rho=%8d, mags=%4d,%4d,%4d\n", 
	rec->zone,rec->id,rec->rho, rec->mB, rec->mR, rec->ci) ;
    */

    for (f=check_fields ; *f; f++) switch((*f)->letter) {
      case_xy:
      	if (!flag) {
    	    o[0] = rec->ra / 3.6e6 ;
    	    o[1] = rec->sd / 3.6e6  - 90. ;
    	    tr_ou(o, u); flag |= 1 ;
      	}
      	if (!(flag&2)) {
      	    tr_uu(u, u, localR) ;
      	    rec->xy[0] = thefields[2].factor*u[1]/u[0] ;
      	    rec->xy[1] = thefields[3].factor*u[2]/u[0] ;
      	    flag |= 2 ;
      	}
	if (rec->xy[i] < (*f)->lim[0]) return(0) ;
	if (rec->xy[i] > (*f)->lim[1]) return(0) ;
      	break ;
      case 'x': i = 0 ; goto case_xy ;
      case 'y': i = 1 ; goto case_xy ;
      case 'a': 
	if ((*f)->lim[0] <= (*f)->lim[1]) {
	    if (rec->ra < (*f)->lim[0]) return(0) ;
	    if (rec->ra > (*f)->lim[1]) return(0) ;
	}
	else {
	    if ((rec->ra < (*f)->lim[0]) &&
	        (rec->ra > (*f)->lim[1])) return(0) ;
	}
	continue ;
      case 'm':
	if (rec->mag < (*f)->lim[0]) return(0) ;
	if (rec->mag > (*f)->lim[1]) return(0) ;
	continue ;
      case 'c':
	if (rec->nc < (*f)->lim[0]) return(0) ;
	if (rec->nc > (*f)->lim[1]) return(0) ;
	continue ;
      case 'd':
	if (rec->sd < (*f)->lim[0]) return(0) ;
	if (rec->sd > (*f)->lim[1]) return(0) ;
	continue ;
      case 'e':
	if (rec->epoch < (*f)->lim[0]) return(0) ;
	if (rec->epoch > (*f)->lim[1]) return(0) ;
	continue ;
      case 'i':
        if (rec->id    < (*f)->lim[0]) return(0) ;
        if (rec->id    > (*f)->lim[1]) return(stopid) ;
	continue ;
      case 'p':
	if (rec->pmtot < 0) {
	    u[0] = rec->pmra ;
	    u[1] = rec->pmsd ;
	    rec->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ;
	}
	if (rec->pmtot < (*f)->lim[0]) return(0) ;
	if (rec->pmtot > (*f)->lim[1]) return(0) ;
	continue ;
      case 'r':
	if (rec->rho < (*f)->lim[0]) return(0) ;
	if (rec->rho > (*f)->lim[1]) return(0) ;
	continue ;
    }
    return(1) ;
}

/*==================================================================
		Comparison routines
 *==================================================================*/

static int compare(UCACrec *a, UCACrec *b)
{
  struct s_fields **f;
  double u[3] ;
    for (f=compare_fields ; *f; f++) switch((*f)->letter) {
      case 'x':
	return((a->xy[0] - b->xy[0]) * (*f)->order) ;
      case 'y':
	return((a->xy[1] - b->xy[1]) * (*f)->order) ;
      case 'a': 
	return((a->ra - b->ra) * (*f)->order) ;
      case 'c':
	return((a->nc - b->nc) * (*f)->order) ;
      case 'i':
	return((a->id - b->id) * (*f)->order) ;
      case 'o':
	return((a->no - b->no) * (*f)->order) ;
      case 'd':
	return((a->sd - b->sd) * (*f)->order) ;
      case 'm':
	return((a->mag - b->mag) * (*f)->order) ;
      case 'p':
	if (a->pmtot < 0) {
	    u[0] = a->pmra ; u[1] = a->pmsd ;
	    a->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ;
	}
	if (b->pmtot < 0) {
	    u[0] = b->pmra ; u[1] = b->pmsd ;
	    b->pmtot = 0.5 + sqrt(u[0]*u[0] + u[1]*u[1]) ;
	}
	return((a->pmtot - b->pmtot) * (*f)->order) ;
      case 'e':
	return((a->epoch - b->epoch) * (*f)->order) ;
      case 'r':
	return((a->rho - b->rho) * (*f)->order) ;
    }
    return(0) ;
}

/*==================================================================
		Sorting the Records
 *==================================================================*/

static int add_rec(UCACrec *new)
/*++++++++++++++++
.PURPOSE  Add a new record, link it.
.RETURNS  0 (not kept) / 1
-----------------*/
{
  LUCAC *curr, *prec, *next ;
  int k, n, dif, pdif ;

    /* (0) Check if the new record has any use... */
    if ((irec == mrec) && (compare(new, &(last->rec)) >= 0)) {
        truncated++ ;
        return(0) ;
    }

    curr = prec = (LUCAC *)0 ;
    /* Compare with record added just before -- can be a good starting point */
    if (just_added) {
        pdif = compare(new, &(just_added->rec));
	if (pdif <= 0) curr = just_added ;
	if ((pdif == 0) && just_added->prev)
	    prec = just_added, curr = just_added->prev; 
    }
    else pdif = 9999 ;
 
    /* When there is a large number of records,
       use random checks to find out where the new record is
       to be inserted. The number of random checks is log2(Nrecs).
       This choice of the starting point (curr) requires that the 
       compare function returns a value which is a scalable difference
       between two records.
       In this test, prec is set ONLY when when two records are identical.
       This means that no further displacement in the linked list is needed.
    */
    if (pdif && (irec & (~255))) for (k=(irec>>3); k; k >>= 1) {
	n = rand()%irec ;
	next = &(arec[n]) ;
	dif = compare(new, &(next->rec)) ;
	/* printf("....Random #%d, dif=%d\n", n, dif) ; */
	if (dif  > 0) continue ;
	if (dif == 0) { prec = next; curr = next->prev; break; }
	if (curr && (dif <= pdif)) continue ; /* I'm away ! */
	curr = next ; pdif = dif ;
    }
    if (!curr) curr=last ;

    /* (1) Find where in the list we've to insert the new record */
    n = 0 ;
    if (!prec) while (curr && (compare(new, &(curr->rec)) < 0)) 
	prec=curr, curr=curr->prev, n++ ;
    /* printf("....(tested %d links)\n", n) ;	*/

    /* (2) If max attained, put the new in place of last */
    if (irec == mrec) {
	truncated++ ;
	next = last ;	
	if (prec == last) prec = (LUCAC*)0;
	else last = last->prev ;
    }
    else next = arec + irec++ ;

    /* (3) Insert the new record, and set the links */
    just_added = next ;
    *(&(next->rec)) = *new ;
    next->prev = curr;
    if (prec) prec->prev = next;
    else last = next ;

    if (ucac_options&1) {
#if 0
	fprintf(stderr, "----irec  = %3d %08d%7d\n", 
	    irec, just_added->rec.id, just_added->rec.rho);
	for (prec=last; prec; prec=prec->prev) fprintf(stderr, 
	    "%08d%7d,", prec->rec.id,  prec->rec.rho) ;
	fprintf(stderr, "\n") ;
#endif
    }
    return(1) ;
}

/*==================================================================
		Search in a Circle
 *==================================================================*/

static int digest(UCACrec *rec) 
/*++++++++++++++++
.PURPOSE  Digest routine: check and display if necessary
.RETURNS  0
-----------------*/
{
  int st ;
    st = check(rec) ;
    if (st <= 0) return(st) ;
    matched += 1 ;
    if (compare_fields[0]) 
    	add_rec(rec)  ;
    else {
	if (irec >= mrec) { truncated++; matched=0; return(-1) ; }
	irec++ ; 
    	puts(ucac2a(rec, opted));
    }
    return(0) ;
}

long ucac_loop()
/*++++++++++++++++
.PURPOSE  Loop on read & test of UCAC records
.RETURNS  Number of tested records.
-----------------*/
{
  long tested = 0 ;
  UCACrec rec ;
    while(1) {
        if (ucac_read(&rec) <= 0) break ;
	tested++;
	if (digest(&rec) < 0) break ;
    }
    return(tested) ;
}

long ucac_pos(long gra[2], long gsd[2])
/*++++++++++++++++
.PURPOSE  Launch Search on UCAC stars from limits in RA + DE
.RETURNS  Number of tested records.
.REMARKS  Merge with specified RA / DE limits
-----------------*/
{
  long ra[2], sd[2], sra[4], tra[4], tested ;
  int i, j ;

    tested = 0 ;
    if ((!thefields[5].selected) && (!thefields[4].selected))
       return(ucac_search(gra, gsd, digest)) ;

    if (gsd) sd[0] = gsd[0], sd[1] = gsd[1] ; 
    else     sd[0] = 0,      sd[1] = adeg(180) ;
    if (thefields[5].selected) {	/* Limits in DE */
	sd[0] = MAX(sd[0], thefields[5].lim[0]) ;
	sd[1] = MIN(sd[1], thefields[5].lim[1]) ;
	if (sd[0] > sd[1]) return(0);	/* Mismatch Dec */
    }

    sra[0] = tra[0] = 0;
    sra[1] = tra[1] = adeg(360)-1 ;	/* Selected RA	*/
    sra[2] = tra[2] = sra[3] = tra[3] = -1 ;
    if (gra) {
	tra[0] = gra[0];
	if (gra[0] <= gra[1]) 
	     tra[1] = gra[1] ;
	else tra[3] = gra[1], tra[2] = 0 ;
    }

    if (thefields[4].selected) {	/* Limits in RA */
	sra[0] = thefields[4].lim[0] ;
	if (thefields[4].lim[0] <= thefields[4].lim[1]) 
	     sra[1] = thefields[4].lim[1] ;
	else sra[3] = thefields[4].lim[1], sra[2] = 0 ;
    }

    for (i=0; i<4; i += 2) for (j=0; j<4; j += 2) {
	if (sra[j] < 0) continue ;
	if (tra[i] < 0) continue ;
	ra[0] = MAX(tra[i], sra[j]) ;
	ra[1] = MIN(tra[i+1], sra[j+1]) ;
	if (ra[0] > ra[1]) continue ;
	tested += ucac_search(ra, sd, digest) ;
    }

    return(tested) ;
}

long ucac_center(double center[2])
/*++++++++++++++++
.PURPOSE  Search UCAC stars from a central position:
	- either within circle of radius (degrees)
	- or within box of half-dimensions boxy
.RETURNS  Number of tested records.
-----------------*/
{
  long ra[2], sd[2];
  unsigned char *p ;
  double value, sr, cosdec, da, dd, maxrad ;

    ra[0] = 0, ra[1] = adeg(360) - 1 ;
    dd = 180. ;

    if ((radius==0) && (boxy[1]==0)) radius = RADIUS ;
    if (radius > 0) dd = maxrad = radius ;
    if (boxy[1] > 0) 		/* Search in a Box */
	dd = maxrad = sqrt(boxy[0]*boxy[0] + boxy[1]*boxy[1]) ;

    /* Don't accept too large searches... */
    if (dd > 45.) {
	fprintf(stderr, "****Radius [%.5f]deg or Box [%.5fx%.5f] too large\n",
	    maxrad, boxy[0], boxy[1]) ;
	exit(1) ;
    }
 
    /* Derive first the limits in Dec and RA */
    value = center[1] - dd ; sd[0] = (90.+value) * 3.6e6 ;
    value = center[1] + dd ; sd[1] = (90.+value) * 3.6e6 ;
    if (sd[0] < 0) sd[0] = 0 ;
    if (sd[1] > adeg(180)) sd[1] = adeg(180) ;

    cosdec = cosd(center[1]) ; 
    sr = sind(maxrad) ;
    du2max = 2. * sind(maxrad/2.) ;

    /* Compute limits on RA */
    if (sr < cosdec)  {		/* Pole not included in circle */
    	da = asind(sr/cosdec) ;
	ra[0] = (center[0] - da)*3.6e6 ;
	if (ra[0] < 0) ra[0] += adeg(360) ;
	ra[1] = (center[0] + da)*3.6e6 ;
	ra[1] %= adeg(360) ;
    }

    /* Define the constants */
    tr_oR(center, localR) ;
    du2max *= du2max ;

    /* Display the Constraints */
    if (ucac_options&1) {
	printf("#....ucac_center(%.5f, %.5f) maxrad=%.5f, limits_mas:\n", 
	    center[0], center[1], maxrad) ;
	printf("#.... (mas) ra=[%ld,%ld]", ra[0], ra[1]) ;
	printf(" sd=[%ld,%ld]\n", sd[0], sd[1]) ;
    }

    /* Merge the -la & -ld limits, and Launch Search */
    return(ucac_pos(ra, sd)) ;
}

/*==================================================================
		Main Program
 *==================================================================*/

/* Signal: if interrupted, return the signal number */
#include <signal.h>
void OnSignal(signo)
{
    if (isatty(2)) 
	fprintf(stderr, "\n****Signal #%d received, STOP\n", signo) ;
    exit(signo) ;
}

main (argc, argv) int argc; char **argv;
{
  char line[BUFSIZ], *p, *a ;
  double center[2] ;
  double flims[2] ;
  char *the_center, *pgm, argline[BUFSIZ] ;
  struct s_fields *f ;
  int goon = 1 ;
  long tested, along, id ;
  int i, sign ;
  int non_flagged_arg = 0 ;
  UCACrec rec ;
  FILE *input_file ;

    if (argc < 2) {
	fprintf(stderr, "%s%s", usage, help) ;
	exit(1) ;
    }

    for (i=1; i<35; i++) signal(i, OnSignal) ;

    /* Keep program name to define the default UCAC root name */
    pgm = argv[0] ;
    /* fprintf(stderr, "....pgm=%s\n", pgm ? pgm : "(nil)") ; */

    /* NOTE: The options can also be  RA-DEC Radius */
    the_center = line ;
    input_file = stdin ;
    while (--argc > 0) {
	p = *++argv;
	if (*p == '-') switch(p[1]) {
	  case 'H':	/* HELP */
	    printf("%s", HELP);
	    exit(0) ;
	  case 'h':	/* Help */
	    fprintf(stderr, "%s%s", usage, help) ;
	    exit(0) ;
	  case 'f':	/* Next parameter = input file */
	    argc--, argv++;
	    input_file = fopen(*argv, "r") ;
	    if (!input_file) { perror(*argv) ; exit(1) ; }
	    continue ;
	  case 'b':	/* Next parameter = Box limits */
	    argc--, argv++;
	    switch(p[2]) {
	      case 'd': _factor_ = 1; break ;
	      case 's': _factor_ = 1./3600. ; break ;
	      default:  _factor_ = 1./60. ; break ;
	    }
	    if (get_lim(*argv, flims) < 2) flims[1] = flims[0] ;
	    _factor_ = 1.; 
	    flims[0] /= 2.0 ; flims[1] /= 2.0 ; 
	    boxy[0] = flims[0] ; boxy[1] = flims[1] ;
	    thefields[2].lim[1] = thefields[2].factor*tand(flims[0]) ;
	    thefields[3].lim[1] = thefields[3].factor*tand(flims[1]) ;
	    thefields[2].lim[0] = -thefields[2].lim[1] ;
	    thefields[3].lim[0] = -thefields[3].lim[1] ;
	    thefields[2].selected = thefields[3].selected = 1 ;
	    continue ;
	  case 'R':	/* Root Name of Catalogue  */
	    if (p[2])  p += 2;
	    else p = *++argv, argc-- ;
	    a = malloc(12+strlen(p)) ;
	    sprintf(a, "UCACroot=%s", p) ;
	    putenv(a) ;
	    pgm = (char *)0 ;	/* Forget about implied from program name */
	    continue ;
	  case 'r':	/* Next parameter = Radius */
	    argc--, argv++;
	    switch(p[2]) {
	      case 'd': _factor_ = 1; break ;
	      case 's': _factor_ = 1./3600. ; break ;
	      default:  _factor_ = 1./60. ; break ;
	    }
	    if (get_lim(*argv, flims) < 2) {
		flims[1] = flims[0] ;
		flims[0] = 0 ;
	    }
	    radius = flims[1] ;
	    _factor_ = 1.; 
	    thefields[0].lim[0] = thefields[0].factor*flims[0] ;
	    thefields[0].lim[1] = thefields[0].factor*flims[1] ;
	    continue ;
	  case 'm':	/* Max number of records */
	    argc--, argv++;
	    mrec = atoi(*argv) ;
	    if (mrec < 1) mrec = 1 ;
	    continue ;
	  case 'e':	/* Edit option	*/
	    p += 2;	/* Is the argument of sort glued with -e ? */
	    if (!*p) p = *++argv, argc-- ;
	    opted = 0 ;
	    if (isdigit(*p)) opted = atoi(p) ;
	    else while(*p) switch(*(p++)) {
		case '.': case ' ': case 'd': continue ;
		case 's': opted |= 1 ; continue ;
		case 'i': opted |= 2 ; continue ;
		case 'e': opted |= 4 ; continue ;
		case 'x': opted |= 8 ; continue ;
		case 'm': opted |= 16; continue ;
		default:
		   fprintf(stderr, "****Bad argument: -e '%s'\n%s", --p, usage);
		    exit(1) ;
	    }
	    continue ;
	  case 'i':	/* Input ID */
	    input_id = 1 ;
	    thefields[1].selected = 1;
	    if (argc > 1) {	/* ID as argument */
		stopid = -1 ;	/* Indicates to stop after last ID read */
		argc--, argv++;
		the_center = strdup(*argv) ;
	    }
	    continue ;
	  case 'z':	/* A Zone */
	    input_id = 2 ;
	    if (argc > 1) {	/* Zone as argument */
		ucac_stop(1) ;	/* Must stop at EOF */
		argc--, argv++;
		the_center = strdup(*argv) ;
	    }
	    continue ;
	  case 'l':	/* Set up the limits */
	    argc--, argv++;
	    i = get_lim(*argv, flims);
	    if (i==1) flims[1] = flims[0] ;
	    if (i==2) flims[0] = -0x7fffffff ;
	    if (!(f = get_field(thefields, p+2))) {
		fprintf(stderr, "****Unknown field in %s\n%s", p, usage);
		exit(1) ;
	    }
	    f->selected = 1 ;
	    f->lim[0] = flims[0] * f->factor ;
	    f->lim[1] = flims[1] * f->factor ;
	    if (p[2] == 'd') 	/* Limits in Declination */
		f->lim[0] += adeg(90), f->lim[1] += adeg(90) ;
	    if ((f->lim[0] > f->lim[1]) && (p[2] != 'a'))
		along = f->lim[0], f->lim[0] = f->lim[1], f->lim[1] = along ;
	    continue ;
	  case 's':	/* Sort Order	*/
	    p += 2;	/* Is the argument of sort glued with -s ? */
	    if (!*p) p = *++argv, argc-- ;
	    for (i = 0 ; *p && (i < ITEMS(compare_fields)); i++, p++) {
		sign = 1 ;
		if (*p == '-') sign = -1, p++ ;
		else if (*p == '+') p++ ;
		f = get_field(thefields, p) ;
		if (!f) {
		    fprintf(stderr, "****Unknown field in -s %s\n%s", 
		    *argv, usage) ;
		    exit(1) ;
		}
	    	f->order    = sign ;
		compare_fields[i] = f ;
	    }
	    continue ;
	  case 'c':	/* Get Center	*/
	    argc--, argv++;
	    if (argc < 1) {
		fprintf(stderr, "****Unspecified center in -c argument\n%s",
		    usage) ;
		exit(1) ;
	    }
	    strcpy(argline, *argv);
	    the_center = argline ;
	    while (argc > 1) {		/* Position as several arguments ? */
		if (issign(argv[1][0])) {
		    if (! isdigit(argv[1][1])) break ;
		}
		else if (! isdigit(argv[1][0])) break ;
		argc--, argv++ ;
		strcat(argline, " ") ;
		strcat(argline, *argv) ;
	    }
	    continue ;
	  case 'v':	/* Verbose	*/
	    ucac_options |= 1 ;
	    continue ;
	  case 'w':	/* Whole sky	*/
	    if (strcmp(p, "-whole") == 0) { 
	        whole_ucac = 1 ; 
	        mrec = 999999999 ;
		continue ; 
	    }
	    /* NO BREAK */
	  default:
	    fprintf(stderr, "****Unknown argument: %s\n%s", p, usage);
	    exit(1) ;
	}
	/* Non-Option: if  */
	if (isdigit(*p) && (non_flagged_arg < 2)) {
	    non_flagged_arg++ ;
	    if (the_center == line) { the_center = *argv; continue ; }
	    get_lim(*argv, flims);
	    radius = flims[0]/60. ;
	    continue ;
	}
	fprintf(stderr, "****Unknown argument: %s\n%s", p, usage);
	exit(1) ;
    }

    /* Define UCAC from the program name */
    if (pgm && (!getenv("UCACroot"))) {
	a = malloc(12+strlen(pgm)) ;
	sprintf(a, "UCACroot=%s", pgm) ;
	for (p = a + strlen(a) -1 ; (p>a) && (*p != '/'); p--) ;
	if (p>a) for (--p; (p>a) && (*p != '/'); p--) ;
	if (strncmp(p, "/bin/", 5) == 0) {
	    strcpy(p, "/bindat") ;
	    putenv(a) ;
	}
	else free(a) ;
    }

    /* Set up the check_fields */
    for (f=thefields, i=0; f->letter; f++)
	if (f->selected) check_fields[i++] = f ;

    /* Set up the array of results */
    if (compare_fields[0]) arec = (LUCAC *)malloc(mrec*sizeof(LUCAC)) ;

    while (goon) {
	if (whole_ucac) goon = 0 ;
	else if (the_center == line) {
    	    if (isatty(0)) printf("----Give %s: ",  prompt[input_id]) ;
	    if (!fgets(line, sizeof(line), input_file)) break ;
	    p = line + strlen(line) ;
	    while ((p>line) && (iscntrl(p[-1]))) p-- ;
	    *p = 0 ;
	    if (strncasecmp(line, "quit", 4) == 0) break ;
	    if (ispunct(line[0])) { puts(line); continue ; }
	}
	else goon = 0 ;
    	irec = truncated = tested = matched = 0 ;

	if (whole_ucac) {
	    printf("#UCAC %s\n", "(whole)") ;
	    puts(ucac_head(opted)) ;
	    tested = ucac_pos((long *)0, (long *)0);	
	}

	else switch(input_id) {

	  case 1:		/* Get from UCAC-IDs */
	    ucac_options &= ~0x10 ;
	    get_ucac(the_center, thefields[1].lim);
	    printf("#UCAC %s\n", the_center) ;
	    puts(ucac_head(opted)) ;
	    if (ucac_set(thefields[1].lim[0]) == 0)
	        tested = ucac_loop() ;
	    break ;
	    
	  case 2:		/* Get in Zone */
	    ucac_options &= ~0x10 ;
	    printf("#UCAC Zone %s\n", the_center) ;
	    puts(ucac_head(opted)) ;
	    if (ucac_zopen(atoi(the_center)) == 0)
	       tested = ucac_loop() ;
	    break ;

	  case 0:		/* Get from Center */
	    ucac_options &= ~0x10 ;
	    if (compare_fields[0]) {
		if (((compare_fields[0])->letter == 'r') 
		  &&((compare_fields[0])->order > 0)) 
		  ucac_options |= 0x10 ;
	    }
	    if (get_center(the_center, center) < 0) continue ;
	    printf("#Center: %s\n", the_center) ;
	    puts(ucac_head(opted)) ;
	    tested = ucac_center(center) ;
	}

	/*---- List saved records in order (replace backward by forward link) */
    	if (compare_fields[0]) { LUCAC *gn, *gc, *gp;	
	    gn=(LUCAC *)0, gc=last;
	    while (gc) gp=gc->prev, gc->prev=gn, gn=gc, gc=gp;
	    while (gn) puts(ucac2a(&(gn->rec), opted)), gn=gn->prev;
	    last = just_added = (LUCAC *)0;
	}
	if (truncated)  {
            printf("#+++Truncated to %d out of ", mrec) ;
            if (matched)  printf("%ld", matched) ;
            else printf("(not-computed)") ;
            printf(" matches (%ld tested)\n", tested) ;
        }
	else printf("#--- %ld matches (%ld tested)\n", matched, tested);
    }
    ucac_close() ;
    exit(0);
}