/*++++++++++++++
.IDENTIFICATION pmm.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-A2.0 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   01-Mar-1997
.VERSION  1.2   02-May-1997: There was a problem when no radius provided...
.VERSION  1.3   03-Apr-1998: Problems with , between RA/DEC.
			Modified the SORT algorithm
.VERSION  1.4   27-May-1998: BUG in sort fixed !
.VERSION  2.0   20-Oct-1998: Version USNO-A2.0
.VERSION  2.1   12-Nov-2000: Swapping installed to run on Linux
.VERSION  2.2   20-Mar-2001: Accept -la and -ld alone
.VERSION  2.3   06-Mar-2002: Changed message with number of matches
.VERSION  2.4   05-Oct-2003: Use the program path to define the path to data
.COMMENTS       Access to PMM catalogue
	Remember that the environment variable PMMroot provides
	the root of PMM compressed catalogue; its default is
	/PMM2
---------------*/

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

#define RADIUS		(2.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*100)	/* Degree into 10mas */

/* Options for PMM queries.......	*/
int pmm_options = 0 ;	/* 1=verbose, 0x10 for Original Order*/
/* Linked records, to sort them...	*/
typedef struct linked_PMM {
    struct linked_PMM *prev ;
    PMMrec rec ;
} LPMMrec ;
static int mrec = 100 ;
static int irec, truncated ;
static int the_zone ;		/* Current zone	   */
static char input_id ;		/* Input is ID     */
static char whole_pmm ;		/* -whole option   */
LPMMrec *last, *arec;

/*  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 ;
    char *alim ;
    double factor ;
    long lim[2] ;
};
static struct s_fields thefields[] =  {
    {'r',0,0,(char *)0,3.6e5},			/* MUST BE #0 */
    {'i',0,0,(char *)0,1.},			/* MUST BE #1 */
    {'x',0,0,(char *)0,(180./M_PI)*3.6e5},	/* MUST BE #2 */
    {'y',0,0,(char *)0,(180./M_PI)*3.6e5},	/* MUST BE #3 */
    {'a',0,0,(char *)0,3.6e5},			/* Must be #4 */
    {'d',0,0,(char *)0,3.6e5},			/* Must be #5 */
    {'B',0,0,(char *)0,10.},		
    {'R',0,0,(char *)0,10.},
    {'c',0,0,(char *)0,10.},
    {'o',0,0,(char *)0,1000.},
    {'p',0,0,(char *)0,1.},
    {'z',0,0,(char *)0,1.},
    {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 char usage[] = "\
Usage: pmm [-HELP] [-r[s] [min,]radius] [-b[s] x[,y]  [-c center]  \n\
           [-e edit_opt] [-f input_file] [-id zone-USNO_id] [-m max_records] \n\
	   [-l! min,max] [-s !] [-whole]\n\
";
 
static char help[] = "\
  -HELP: display column explanations\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 (defaut in stdin)\n\
     -e: s=edit RA/DE in Sexa; i=edit USNO-ID, e=edit ObsDate, x=edit x,y\n\
     -f: specifies an input file (default stdin)\n\
    -id: choice of a zone and USNO-A2.0 instead of target position\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\
 -whole: search on the whole sky\n\
====The abbreviations of the parameters (symbolized !)are:\n\
      a=Alpha   B=Bmag       c=Color(B-R)   d=Delta    i=USNO-ID  o=obsDate\n\
      p=Plate   r=distance   R=Rmag         x=proj.E   y=proj.N   z=Zone\n\
";

static char HELP[] = "\
Access to USNO-A2.0 (by Dave Monnet <dgm@nofs.navy.mil> Flagstaff Station)\n\
=============================================================================\n\
\n\
The fields have the following meaning:\n\
\n\
USNO2   = Designation in original catalogue, made of a zone number (4 digits\n\
          from 0000 to 1725 representing the distance in 0.1deg to the \n\
          South Pole, stepped by 75 = 7.5degrees), and a number in the \n\
          area (8 digits, order by IRCS (J2000) RA);\n\
	  ****NOTE that this number differs from the the USNO-A2.0****\n\
RA+Dec  = Position (Alpha, Delta) in ICRS (Hipparcos system, more accurate\n\
	  but compatible with J2000), mean epoch of plates.\n\
A*      = flags: A = star has an entry in the ACT Reference Catalog <I/246>\n\
		 * = magnitudes are probably wrong\n\
Bmag    = photographic blue magnitude (O in the North, J in the South)\n\
Rmag    = photographic  red magnitude (E in the North, F in the South)\n\
Epoch   = mean date (Julian years) of the blue and red plates\n\
Pl.     = plate number (positive for POSS, negative for ESO/SERC)\n\
r(\")    = distance from specified center, in arc seconds.\n\
=============================================================================\n\
For details concerning PMM and its products at Flagstaff Station, refer to \n\
        http://www.usno.navy.mil/pmm\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
-----------------*/
{
  int found, n ;
  char *p ;

    found = 0 ; p = string ; values[0] = values[1] = 0 ;
    n = match_num(p, values) ;
    if (n) found |= 1, p += n ;

    if (*p) {
	p++ ; 
	/* if ((*p == ',')||(*p == '/')) p++ ; */
	n = match_num(p, values+1) ;
	if (n) found |= 2 ;
    }
    return(found) ;
}

static int get_pmm (char *string, long *id)
/*++++++++++++++++
.PURPOSE  Interpret a string for the_zone + USNO2_id
.RETURNS  Number of bytes matched
-----------------*/
{
  double v ; char *p ;

    p = string ;
    p += match_unum(p, &v) ;
    if (*p) {	/* Zone is specified */
	the_zone = v ;
	++p;
	p += match_unum(p, &v) ;
    }
    *id = v ;
    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(PMMrec *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.6e5 ;
    	o[1] = rec->sd / 3.6e5  - 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 = 360000. * 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.6e5 ;
    	    o[1] = rec->sd / 3.6e5  - 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 'B':
	if (rec->mB < (*f)->lim[0]) return(0) ;
	if (rec->mB > (*f)->lim[1]) return(0) ;
	continue ;
      case 'c':
	if (rec->ci < (*f)->lim[0]) return(0) ;
	if (rec->ci > (*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 'o':
	if (rec->epoch < (*f)->lim[0]) return(0) ;
	if (rec->epoch > (*f)->lim[1]) return(0) ;
	continue ;
      case 'p':
	if (rec->plate < (*f)->lim[0]) return(0) ;
	if (rec->plate > (*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 ;
      case 'R':
	if (rec->mR < (*f)->lim[0]) return(0) ;
	if (rec->mR > (*f)->lim[1]) return(0) ;
	continue ;
      case 'z':
	if (rec->zone < (*f)->lim[0]) return(0) ;
	if (rec->zone > (*f)->lim[1]) return(0) ;
	continue ;
    }
    return(1) ;
}

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

static int compare(PMMrec *a, PMMrec *b)
{
  struct s_fields **f;
    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 'B':
	return((a->mB - b->mB) * (*f)->order) ;
      case 'c':
	return((a->ci - b->ci) * (*f)->order) ;
      case 'd':
	return((a->sd - b->sd) * (*f)->order) ;
      case 'j':
	return((a->epoch - b->epoch) * (*f)->order) ;
      case 'p':
	return((a->plate - b->plate) * (*f)->order) ;
      case 'r':
	return((a->rho - b->rho) * (*f)->order) ;
      case 'R':
	return((a->mR - b->mR) * (*f)->order) ;
      case 'z':
	return((a->zone - b->zone) * (*f)->order) ;
    }
    return(0) ;
}

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

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

    /* (0) Check if the new record has any use... */
    if ((irec == mrec) && (compare(new, &(last->rec)) >= 0)) {
        truncated++ ;
        return(0) ;
    }
 
    /* 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.
    */
    if (irec & (~255)) {
	curr = (LPMMrec *)0 ;
	for (k=irec; k; k >>= 1) {
	    dif = rand()%irec ;
	    /* printf("....Random #%d", dif) ; 	*/
	    next = &(arec[dif]) ;
	    dif = compare(new, &(next->rec)) ;
	    /* printf(", dif=%d\n", dif) ;	*/
	    if (dif > 0) continue ;
	    if (dif == 0) break ;
	    if (curr && (dif < pdif)) continue ; /* I'm away ! */
	    curr = next ; pdif = dif ;
	}
	if (!curr) curr=last, prec=(LPMMrec*)0;
	else if (dif == 0) {
	    prec = curr ;
	    curr = curr ?  curr->prev : last ;
	}
    }
    else curr=last, prec=(LPMMrec*)0;

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

    /* (2) If max attained, remove or return */
    if (irec == mrec) {
	truncated++ ;
	if (!prec) return(0) ;
	next = last ;
	if (prec == last) prec = (LPMMrec*)0;
	else last = last->prev ;
    }
    else next = arec + irec++ ;

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

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

static int digest(PMMrec *rec) 
/*++++++++++++++++
.PURPOSE  Digest routine: check and display if necessary
.RETURNS  0
.REMARKS  matched set to zero when unknown
-----------------*/
{
    if (!check(rec)) return(0) ;
    matched += 1 ;
    if (compare_fields[0]) 
    	add_rec(rec)  ;
    else {
	if (irec >= mrec) { truncated++; matched = 0; return(-1) ; }
	irec++ ; 
    	puts(pmm2a(rec, opted));
    }
    return(0) ;
}

long pmm_pos(long gra[2], long gsd[2])
/*++++++++++++++++
.PURPOSE  Launch Search on USNO 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 (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 += pmm_search(ra, sd, digest) ;
    }

    return(tested) ;
}

long pmm_center(double center[2])
/*++++++++++++++++
.PURPOSE  Search USNO 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.6e5 ;
    value = center[1] + dd ; sd[1] = (90.+value) * 3.6e5 ;
    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.6e5 ;
	if (ra[0] < 0) ra[0] += adeg(360) ;
	ra[1] = (center[0] + da)*3.6e5 ;
	ra[1] %= adeg(360) ;
    }

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

    /* Display the Constraints */
    if (pmm_options&1) {
	printf("#....pmm_center(%.5f, %.5f) maxrad=%.5f, limits_10mas:\n", 
	    center[0], center[1], maxrad) ;
	printf("#....       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(pmm_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 ;
  double center[2] ;
  double flims[2] ;
  char *the_center, argline[BUFSIZ] ;
  struct s_fields *f ;
  int goon = 1 ;
  long tested, along, id ;
  int i, sign ;
  int non_flagged_arg = 0 ;
  PMMrec rec ;
  FILE *input_file ;

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

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

    /* Define PMMroot from the program name */
    if (strchr(argv[0],'/') && (!getenv("PMMroot"))) {
	sprintf(line, "PMMroot=%s", *argv) ;
	for (p = line + strlen(line) -1 ; (p>line) && (*p != '/'); p--) ;
	if (p>line) for (--p; (p>line) && (*p != '/'); p--) ;
	if (strncmp(p, "/bin/", 5) == 0) {
	    *p = 0 ;
	    putenv(strdup(line)) ;
	}
    }

    /* 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':	/* 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 ;
	    }
	    else 	/* 2 limits ==> must check !! */
		thefields[0].selected = 1 ;
	    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 ;
		default:
		   fprintf(stderr, "****Bad argument: -e '%s'\n%s", --p, usage);
		    exit(1) ;
	    }
	    continue ;
	  case 'i':	/* Input ID numbers  ZZZZ-NNNNN	*/
	    input_id = 1 ;
	    if ((argc > 1) && (*argv[1] != '-')) {	/* ID as argument */
		++argv; --argc;
		the_center = strdup(*argv) ;
		/* It may happen that another separator is used that the -
		   for instance a blank or a dot. Just force the separator to -
		*/
		for (p=the_center; isdigit(*p); p++) ;
		if (*p) *p = '-' ; 
	    }
	    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->alim   = *argv ;
	    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-- ;
	    if (*p == '.') 	/* No Sort at All */
		continue ;
	    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	*/
	    pmm_options |= 1 ;
	    continue ;
	  case 'w':	/* Whole sky	*/
	    if (strcmp(p, "-whole") == 0) { whole_pmm = 1 ; 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) ;
    }

    /* 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 = (LPMMrec *)malloc(mrec*sizeof(LPMMrec)) ;

    while (goon) {
	if (whole_pmm) goon = 0 ;
	else if (the_center != line) 	/* -c was specified ! */
	    goon = 0 ;
	else if (thefields[4].selected && thefields[5].selected)
	    the_center = (char *)0, 
	    goon = 0 ;			/* -la and -ld specified */
	else if (the_center == line) {
    	    if (isatty(0)) printf("----Give %s: ", 
		input_id ? "Zone-USNO2_id" : "Center") ;
	    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_pmm) tested = pmm_pos((long *)0, (long *)0);	

	else if (input_id) {	/* Get from PMM-IDs */
	    get_pmm(the_center, &id) ;
	    printf("#PMM %s\n", the_center) ;
	    puts(pmm_head(opted)) ;
	    if (pmm_set(the_zone, id>0 ? id : 1) < 0) continue ;
	    pmm_read(&rec) ; tested++;
	    digest(&rec) ;
	    while (id == 0) {
		if (pmm_read(&rec) <= 0) break ; 
		tested++;
		if (digest(&rec) < 0) break ;
	    }
	}

	else if (!the_center) {		/* -la and -ld only  */
	    pmm_options |= 0x10 ;	
	    printf("#Limits: -la=%s -ld=%s\n", 
		thefields[4].alim, thefields[5].alim) ;
	    puts(pmm_head(opted)) ;
	    tested = pmm_pos(thefields[4].lim, thefields[5].lim) ;
	}

	else {			/* Get from Center   */
	    if (get_center(the_center, center) < 0) continue ;
	    printf("#Center: %s\n", the_center) ;
	    puts(pmm_head(opted)) ;
	    tested = pmm_center(center) ;
	}

	/*---- List saved records in order (replace backward by forward link) */
    	if (compare_fields[0]) { LPMMrec *gn, *gc, *gp;	
	    gn=(LPMMrec *)0, gc=last;
	    while (gc) gp=gc->prev, gc->prev=gn, gn=gc, gc=gp;
	    while (gn) puts(pmm2a(&(gn->rec), opted)), gn=gn->prev;
	    last = (LPMMrec *)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);
    }
    pmm_close() ;
    exit(0);
}