/*++++++++++++++
.IDENTIFICATION nomad_sub.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO Catalogues
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   02-Aug-2006: Adapted from USNO
.VERSION  1.1   17-Feb-2007: Problem with RA close to limit.
.VERSION  1.2   31-Jul-2007: edition position only
.VERSION  1.3   03-Mar-2008: Bug for RA selection!!
.VERSION  1.31  05-May-2011: HIP number can be > 1000000 !!!
.COMMENTS       Edition of NOMAD data.

---------------*/

#include <nomad1def.h>	/* Structure definitions */
#include <nomad1.h>	/* Structure definitions */

#include <stdio.h>
#include <stdlib.h>	/* Malloc   */
#include <unistd.h>
#include <string.h>
#include <ctype.h>
#include <fcntl.h>	/* O_BINARY */
#include <math.h>
#define  sind(x)	sin(x*(M_PI/180.))
#define  cosd(x)	cos(x*(M_PI/180.))
#define  asind(x)	(180./M_PI)*asin(x)

#ifndef O_BINARY	/* Required for MS-DOS	*/
#define O_BINARY 0
#endif
#include <stdio.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))

typedef int (*SORT_FCT)(const void *, const void *) ;

/* Options for NOMAD queries.......	*/
int nomad_options ;		/* 0x10 for OrigOrder	*/

static char *oval_name[3] = {
  "pm=" /* 0 */, "mag="   /* 1 */, 
  "Ep=" /* 2 */
};

/* Sources defined in the catalog: USNO-B, 2MASS, YB6, UCAC2, Tycho, Hip 
 *  	Oplate, Eplate
 * */
static char source_letter[] = " BMYCTH.oe"; 

static struct {
    char *root;			/* Root name of NOMAD	*/
    char *name ;		/* Name of opened file	*/
    int fno;
    char line[Lheader] ;	/* First line of file	*/
    short SDz;			/* Zone opened 		*/
    char reclen;		/* 36 or 37		*/
    short lename;		/* Length of name buffer*/
    short nchunks;		/* Actual # chunks	*/
    int mbuf ;			/* Allocated bytes mbuf	*/
    unsigned char *abuf;	/* Loaded data		*/
    int4 oval[4];		/* Offsets in oval_name */
    int4 chunks[2*Nchunks];	/* List (filepos, ID)	*/
    int4 *echunks;		/* (Last+1) valid chunk	*/
    CHUNK dat;			/* All details of Chunk	*/
} NOMADcat ;

static int stopat_eof ;		/* 1=stop at EOF, 2=stop at EndOfChunk */
static int swapping ;		/* Set to 1 or 2 (swap)	*/
static int the_status;		/* Last Status Read     */

/* The first byte of a data record is related to the record length;
    the implied length of the variable part is defined here.
*/
static unsigned char varlen[256] = {
     0, 2, 2, 4, 2, 4, 4, 6, 2, 4, 4, 6, 4, 6, 6, 8,
     2, 4, 4, 6, 4, 6, 6, 8, 4, 6, 6, 8, 6, 8, 8,10,
     2, 4, 4, 6, 4, 6, 6, 8, 4, 6, 6, 8, 6, 8, 8,10,
     4, 6, 6, 8, 6, 8, 8,10, 6, 8, 8,10, 8,10,10,12,
     4, 6, 6, 8, 6, 8, 8,10, 6, 8, 8,10, 8,10,10,12,
     6, 8, 8,10, 8,10,10,12, 8,10,10,12,10,12,12,14,
     6, 8, 8,10, 8,10,10,12, 8,10,10,12,10,12,12,14,
     8,10,10,12,10,12,12,14,10,12,12,14,12,14,14,16,
     4, 6, 6, 8, 6, 8, 8,10, 6, 8, 8,10, 8,10,10,12,
     6, 8, 8,10, 8,10,10,12, 8,10,10,12,10,12,12,14,
     6, 8, 8,10, 8,10,10,12, 8,10,10,12,10,12,12,14,
     8,10,10,12,10,12,12,14,10,12,12,14,12,14,14,16,
     8,10,10,12,10,12,12,14,10,12,12,14,12,14,14,16,
    10,12,12,14,12,14,14,16,12,14,14,16,14,16,16,18,
    10,12,12,14,12,14,14,16,12,14,14,16,14,16,16,18,
    12,14,14,16,14,16,16,18,14,16,16,18,16,18,18,20
} ;

/*==================================================================
		Internal Utilities
 *==================================================================*/

static int strloc(char *text, int c)
/*++++++++++++++++
.PURPOSE  Locate specified character
.RETURNS  Index of located char
-----------------*/
{
  char *s;
    for ( s = text; *s; s++)      if (*s == c)    break;
    return(s-text);
}

static int get(unsigned char *a, int b, int len)
/*++++++++++++++++
.PURPOSE  Get a value made of 'len' bits
          starting from byte position 'a', bit offset 'b'.
.RETURNS  The value
.REMARKS  Independent of the architecture.
-----------------*/
{
   static unsigned int mask[] = {
    0x00, 0x01, 0x03, 0x07, 0x0f, 0x1f, 0x3f, 0x7f,
    0x00ff, 0x01ff, 0x03ff, 0x07ff, 0x0fff, 0x1fff, 0x3fff, 0x7fff,
    0x00ffff, 0x01ffff, 0x03ffff, 0x07ffff, 0x0fffff,
              0x1fffff, 0x3fffff, 0x7fffff,
    0x00ffffff, 0x01ffffff, 0x03ffffff, 0x07ffffff, 0x0fffffff,
                0x1fffffff, 0x3fffffff, 0x7fffffff,
    0xffffffff } ;
  int value;
  unsigned char *ac;
  int nb ;      /* remaining bits to get */

    ac = a + (b>>3);    /* Byte position  */
    value = *(ac++);    /* Initialisation */
    nb = len;           /* Bits to read.  */
    nb -= 8 - (b&7);    /* Useful bits    */

    /* We don't care about the leftmost bits,
       these will be removed at the end
       (we assume that len <= 32...)
    */

    while (nb >= 8) { value = (value<<8) | *(ac++); nb -= 8; }
    if (nb < 0)         /* I've read too much */
        value >>= (-nb);
    else if (nb)        /* Remainder bits  */
        value = (value<<nb) | (*ac >> (8-nb));

    return(value&mask[len]);
}

static void swap2 (short *array, int nshort)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of shorts if necessary
.RETURNS  0/1/2 (type of swap)
.REMARKS  Useful for big-endian machines
-----------------*/
{
  register int m, n ;
  register short *p, *e ;
    for (p=array, e=p+nshort; p<e; p++) {
        n = *p ;
	m =  (n>>8) & 0xff ;
	m |= (n<<8) ;
	*p = m ;
    }
}

#define swap4(a,n)	swap((int *)(a), n)
static void swap (int *array, int nint)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of integers if necessary
.RETURNS  0/1/2 (type of swap)
.REMARKS  Useful for big-endian machines
-----------------*/
{
  char *p, *e ;
  int n ;
  
    p = (char *)array ; e = p + 4*nint ;
    if (swapping == 1) swap2((short *)array, 2*nint) ;
    else if (swapping == 2) while (p < e) {
        n = p[0] ; p[0] = p[3] ; p[3] = n ;
        n = p[1] ; p[1] = p[2] ; p[2] = n ;
        p += 4 ;
    }
}


static int4 *binloc(int4 *start, int4 *end, int4 value, int step)
/*++++++++++++++++
.PURPOSE  Find index in table such that    start[i] <= value < start[i+step]
.RETURNS  NULL if pointer < *start / the pointer otherwise
-----------------*/
{
  int4 *low, *high, *med ;
  int i;

    if (value < *start) return((int4 *)0) ;
    /* Compute the highest item */
    i = (end-start)%step;	
    if (i == 0) i = step;
    low=start, high=end-i;
    if (value >= *high) return(high) ;
    while ((high-low) > step) {
    	med = low + step*((high-low)/2/step) ;
    	if (value < *med) high = med;
    	else low = med ;
    }
    return(low) ;
}

static int4 ra4(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  int4 ra ;
    ra = (pc[4]<<24) | (pc[5]<<16) | (pc[6]<<8) | pc[7] ;
    return(ra) ;
}

static int4 ra3(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  int4 ra ;
    ra = (pc[4]<<16) | (pc[5]<<8) | pc[6];
    return(ra) ;
}

#if 0
static int binlocf(char *start, int bytes, int reclen, int4 (*f)(), int4 value)
/*++++++++++++++++
.PURPOSE  Find index in table such that   f(p) <= value < f(p+reclen)
.RETURNS  Byte position in start / -1 / bytes
-----------------*/
{
  char *low, *high, *med ;
    if (value < (*f)(start)) return(-1) ;
    low=start, high=start+bytes-reclen;
    if (value >= (*f)(high)) return(bytes) ;
    while ((high-low) > reclen) {
    	med = low + reclen*((high-low)/reclen/2) ;
    	if (value < (*f)(med)) high = med;
    	else low = med ;
    }
    return(low - start) ;
}
#endif

/*==================================================================
		NOMAD File Manipulation
 *==================================================================*/
static void init(char *root)
/*++++++++++++++++
.PURPOSE  Initialisation of NOMAD access.
.RETURNS  ---
-----------------*/
{
  static int value ;
  char *v ;
  int n ;

    if (NOMADcat.root) return ; 		/* Already opened... */

    /* Verify the Swapping ! */
  
    value = 0x010203 ;
    v = (char *)(&value) ;
    if ((v[0] == 0) && (v[1] == 1) && (v[2] == 2)) 
        swapping = 0 ; 	/* No swap necessary */
    else if ((v[0] == 1) && (v[1] == 0) && (v[2] == 3)) 
        swapping = 1 ; 	/* Half-word swap */
    else if ((v[0] == 3) && (v[1] == 2) && (v[2] == 1)) 
        swapping = 2 ; 	/* Full-word swap */
    else {
        fprintf(stderr, "****Irrationnal Byte Swap %02x%02x%02X%02x\n", 
            v[0]&0xff, v[1]&0xff, v[2]&0xff, v[3]&0xff) ;
	exit(2) ;
    }
    if (swapping && (nomad_options&1)) fprintf(stderr, 
        "#...swapping type=%d\n", swapping) ;

    /* Choose the root name */
    if (root)  NOMADcat.root = root ;
    if (!NOMADcat.root) NOMADcat.root = getenv("NOMADroot") ;
    if (!NOMADcat.root) NOMADcat.root = "/NOMAD" ;
    n = 20 + strlen(NOMADcat.root);	/* Size of file name */
    if (n <  48) n = 128 ;
    NOMADcat.name = malloc(NOMADcat.lename = n) ;
    NOMADcat.fno = NOMADcat.SDz  = NOMADcat.dat.RAz = -1 ;
    memset(&NOMADcat.dat, 0, sizeof(CHUNK)) ;
}

static int close_chunk()
/*++++++++++++++++
.PURPOSE  Close the currently opened chunk
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    memset(&NOMADcat.dat, 0, sizeof(CHUNK));
    NOMADcat.dat.ppos= -1;	/* Unknown */
    NOMADcat.dat.RAz = -1;
    return(0) ;
}

int nomad_close()
/*++++++++++++++++
.PURPOSE  Close the currently opened file
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    if (NOMADcat.fno < 0) 	/* Was never opened! */
	return(0);
    close_chunk();
    /* Remove Information on existing chunks */
    NOMADcat.nchunks = 0; NOMADcat.echunks = NOMADcat.chunks;
    memset(NOMADcat.chunks, 0, sizeof(NOMADcat.chunks)) ;
    if (close(NOMADcat.fno) < 0) 
        fprintf(stderr, "****"), perror(NOMADcat.name) ;
    NOMADcat.fno = -1 ;
    NOMADcat.SDz = -1 ;
    NOMADcat.reclen  = 0 ;
    NOMADcat.name[0] = 0 ;
    NOMADcat.line[0] = 0 ;
    return(0) ;
}

static int nomad_fopen(char *name)
/*++++++++++++++++
.PURPOSE  Open a name-specified file
.RETURNS  0 = OK / -1 = error
-----------------*/
{
  char buffer[sizeof(NOMADcat.line)];
  int file, len, i;
  char *p;

    if (nomad_options&1) printf("#....nomad_fopen(%s)\n", name) ;
    if (!NOMADcat.root) init((char *)0) ;
    if (NOMADcat.fno >= 0) {
	fprintf(stderr, "++++nomad_fopen: Closing %s\n", NOMADcat.name);
	nomad_close();
    }
    file = open(name, O_BINARY);
    if (file < 0) {  fprintf(stderr, "****"); perror(name);  return (-1); }
    len = read(file, buffer, sizeof(buffer)) ;
    if (len <= 0)  {  fprintf(stderr, "****"); perror(name);  return (-1); }

    if (NOMADcat.name != name) strncpy(NOMADcat.name, name, NOMADcat.lename);
    NOMADcat.fno  = file;
    memset(NOMADcat.line,   0, sizeof(NOMADcat.line));
    memset(NOMADcat.chunks, 0, sizeof(NOMADcat.chunks));
    memset(&NOMADcat.dat,   0, sizeof(CHUNK));

    /* The end of the first line indicates the starting point
       of the List of Chunks 
    */
    len = strloc(buffer, '\n');
    if (len >= sizeof(NOMADcat.line)) {
	fprintf(stderr,  "****File %s: too long first line (%d>=%d)\n",
	   name, len, (int)sizeof(NOMADcat.line)) ;
	return(-1);
    }
    memcpy(NOMADcat.line, buffer, len);

    /* Example of a header line:
NOMAD-1.0(36) 111/N1111.bin 0000001-0642643  pm=-15000 mag=-5000 Ep=+19010 
UCAC2=+35000000 xtra=7,1
    */

    /* Verify magic */
    i = strloc(header_text, '('); 
    if (strncmp(buffer, header_text, i+1) != 0) {
    	fprintf(stderr, "****File %s: bad magic, is not NOMAD-1.0\n", name);
    	return (-1);
    }
    NOMADcat.reclen = atoi(buffer+i+1);

    p = strchr(buffer+i+2, 'N'); 	/* Find the zone */
    NOMADcat.SDz = atoi(p+1);
    NOMADcat.dat.RAz = -1;

    /* Find the offset parameters (in oval_name) */
    for (i=0; i<ITEMS(oval_name); i++) {
	p = strstr(buffer, oval_name[i]);
	if (!p) {
	    fprintf(stderr,  "****File %s: missing parameter '%s'\n", 
		name, oval_name[i]);
	    return(-1);
	}
	p += strloc(p, '=') +1;
	NOMADcat.oval[i] = atoi(p);
    }

    /* Read the Nchunks*(Offset, ID) designating Chunks. */
    if (read(file, NOMADcat.chunks, sizeof(NOMADcat.chunks)) <= 0) {
	fprintf(stderr, "****"); perror(NOMADcat.name); 
	return(-1);
    }
    if (swapping) swap4(NOMADcat.chunks, 
	    sizeof(NOMADcat.chunks)/sizeof(int));

    /* Find the actual end of Chunks -- there may be zeroes */
    for (i=2; (i<ITEMS(NOMADcat.chunks)) && (NOMADcat.chunks[i]); i += 2) ;
    NOMADcat.nchunks = i>>1;		/* 2 numbers (o,ID) per chunk   */
    NOMADcat.nchunks--;			/* Actually, last indicates EOF */
    NOMADcat.echunks = &(NOMADcat.chunks[i]);

    return(0) ;
}

int nomad_zopen(int zone)
/*++++++++++++++++
.PURPOSE  Open a speficied zone (0..1799)
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    if (nomad_options&1) printf("#...nomad_zopen(N%04d)\n", zone);
    if (!NOMADcat.root) init((char *)0);
    if (NOMADcat.SDz == zone) {
        close_chunk();
	return(0);
    }
    nomad_close();		/* Need another file! */
    sprintf(NOMADcat.name, 
      "%s/%03d/N%04d.bin", NOMADcat.root, zone/10, zone);
    return(nomad_fopen(NOMADcat.name));
}

int nomad_stop(int flag)
/*++++++++++++++++
.PURPOSE  Specify if reading should stop at EOF (1), endofChunck (2), never (0)
.RETURNS  Previous status
-----------------*/
{
  int ostat ;
    ostat = stopat_eof ;
    stopat_eof = flag ;
    return(ostat) ;
}

static int read_chunk(int chunkno)
/*++++++++++++++++
.PURPOSE  Load the specified chunk
.RETURNS  0 = OK / -1 = error
.REMARKS  The largest chunck is ~10Mbytes
	Structure of one chunck on the file:
	+00 (int4)   : Oa = length of preface = offset of accelerator
	+04 (3*int4) : values id_min RA_min SD_min
	+16 (3*int4) : values id_max RA_max SD_max
	+28 (int2)   : N2 = number of short 'extras'
	+30 (int2)   : N4 = number of long  'extras'
	+32 (N4*int4): list of long  'extras'
	+.. (N2*int2): list of short 'extras'
	+Oa (3N*int4): Accelerator, a list of 3N values containing
	               (offset, id, ra)
-----------------*/
{
  int4 o, len, *aint4 ;
  int  ino;
  short *ashort;

    /* The file must be started and loaded */
    if (nomad_options&1) printf("#...read_chunk(%d) in zone N%04d\n",
        chunkno, NOMADcat.SDz) ;
    if (NOMADcat.SDz < 0) return(-1) ;	/* No zone loaded  */

    /* The table of Offsets gives the (location,ID) in the File */
    if (chunkno >= Nchunks-1)	/* Too large chunk number! */
	return(-1);

    /* Accept any chunk number of single-chunk file (pole) */
    if (NOMADcat.reclen > NOMAD1reclen) chunkno = 0;
    ino = chunkno<<1;
      o = NOMADcat.chunks[ino] ;
    len = NOMADcat.chunks[ino+2] - o ;
    if (len <= 0) return(-1) ;

    /* Verify chunck not already loaded -- not necessary to reload ! */
    if (NOMADcat.dat.RAz != chunkno) {

        if (lseek(NOMADcat.fno, o, 0) != o) {
            fprintf(stderr, "****Can't move to %d, file: %s\n",
                o, NOMADcat.name) ;
            return(-1) ;
        }

        /* Need to read the whole chunck -- enough memory ? */
        if (NOMADcat.mbuf < len) {
            if (NOMADcat.abuf) free(NOMADcat.abuf) ;
            NOMADcat.mbuf = 1+(len|0xffff) ;    /* Multiple of 64Kbytes */
            NOMADcat.abuf = malloc(NOMADcat.mbuf) ;
        }
        NOMADcat.dat.len = read(NOMADcat.fno, NOMADcat.abuf, len) ;
        if (NOMADcat.dat.len != len) {
            fprintf(stderr, "****Got %d/%d bytes, ", NOMADcat.dat.len, len) ;
            if (NOMADcat.dat.len < 0) perror(NOMADcat.name);
            else fprintf(stderr, "file: %s\n", NOMADcat.name) ;
            return(-1) ;
        }
	/* NOTE: Buffer was rounded -- the actual size may be 
	   up to 3 bytes less. We therefore diminish the chunk length,
	   but will test more carefully for the end-of-chunk.
	   Actual size given in the last element of the Accelerator
	   (see below)
	NOMADcat.dat.len -= 3;
	*/
    
        aint4 = (int4 *)NOMADcat.abuf;	/* Starts by 7 integers */
        if (swapping) swap4(aint4, 7);	/* Len.preface, minmax	*/
        len = aint4[0];			/* Preface length...	*/

	/* Extra values: aint4[7] contains the numbers of N2 and N4 */
        ashort = (short *)(aint4+7);
        if (swapping) swap2(ashort, 2);
	NOMADcat.dat.ntra2 = ashort[0];
	NOMADcat.dat.ntra4 = ashort[1];
        NOMADcat.dat.xtra4 = aint4+8;
        NOMADcat.dat.xtra2 = (short *)(NOMADcat.dat.xtra4 + ashort[1]);
	if (swapping) {
	    swap2(NOMADcat.dat.xtra2, ashort[0]);
	    swap4(NOMADcat.dat.xtra4, ashort[1]);
	}
    
	NOMADcat.dat.RAz  = chunkno;
	NOMADcat.dat.ppos = -1;
        NOMADcat.dat.id0  = aint4[1];
        NOMADcat.dat.ra0  = aint4[2];	/* = chunkno<<24; */
        NOMADcat.dat.sd0  = aint4[3];   /* = NOMADcat.SDz*360000 */
        NOMADcat.dat.id1  = aint4[4];
        NOMADcat.dat.ra1  = aint4[5];   /* = NOMADcat.dat.ra0 | 0xffffff; */
        NOMADcat.dat.sd1  = aint4[6];   /* = NOMADcat.dat.sd0 + 360000 */
	/* if (NOMADcat.reclen > NOMAD1reclen) 
	    NOMADcat.dat.ra1 = 360*3600*1000 -1;
	*/
    
        /* The Accelerator starts just after the Preface.
	   The number of elements in this index is computed from
	   the position of the data in index[0].
	*/
        NOMADcat.dat.index = (int4 *)(NOMADcat.abuf + len) ;
	memcpy(&o, NOMADcat.abuf + len, sizeof(int));
	if (swapping) swap4(&o, 1);
        NOMADcat.dat.endex = (int4 *)(NOMADcat.abuf + o) ;
        if (swapping) swap4(NOMADcat.dat.index, 
	    NOMADcat.dat.endex - NOMADcat.dat.index);

    	/* Set the actual length from the accelerator */
	aint4 = NOMADcat.dat.endex - 3;
	NOMADcat.dat.len = aint4[0];
    }

    /* Initialize the CHUNK parameters */
    NOMADcat.dat.pos = NOMADcat.dat.index[0];
    NOMADcat.dat.id  = NOMADcat.dat.id0;
    NOMADcat.dat.ppos = -1;

    return(0);
}

static int nomad_next()
/*++++++++++++++++
.PURPOSE  Move to the next redord
.RETURNS  1 = OK / 0 = EOF / -1 = Error
-----------------*/
{
  int len, status;

    if (NOMADcat.dat.pos >= NOMADcat.dat.len) return(0);
    status = NOMADcat.abuf[NOMADcat.dat.pos];
    len = NOMADcat.reclen + varlen[status];
    NOMADcat.dat.ppos = NOMADcat.dat.pos;
    NOMADcat.dat.pos += len;
    NOMADcat.dat.id  += 1  ;
    if (NOMADcat.dat.pos < NOMADcat.dat.len) return(1);
    if (NOMADcat.dat.pos <= NOMADcat.dat.len+3) /* Leave room for align */
	return(0);
    NOMADcat.dat.ppos = -1;	/* We've lost position of preceding record */
    NOMADcat.dat.pos -= len;
    NOMADcat.dat.id -= 1;
    fprintf(stderr, "****nomad_next %04d#%03d, id=%d: %d(%d) > %d\n",
	NOMADcat.SDz, NOMADcat.dat.RAz, NOMADcat.dat.id, 
	NOMADcat.dat.pos, len, NOMADcat.dat.len);
    return(-1);
}

/*==================================================================
		NOMAD File Seek, by Position or Identification
 *==================================================================*/

static int4 aseek(int4 value)
/*++++++++++++++++
.PURPOSE  Position the file to specified RA value (in mas)
.RETURNS  The NOMAD1 number in that zone / 0
-----------------*/
{
  int a0, a1, da;
  int4 *aval;
  int4 (*fra)();	/* Function to get the RA   */

    /* Verify -- maybe we're in the correct chunk ? */
    if ((NOMADcat.dat.RAz >= 0) && 
       (value >= NOMADcat.dat.ra0) && (value <= NOMADcat.dat.ra1)) ;
    else {
	if (read_chunk(value>>24)) return(-1) ;
    }

    if ((value < NOMADcat.dat.ra0) || (value >  NOMADcat.dat.ra1)) 
        return(-1) ;
    da = value - NOMADcat.dat.ra0 ;
    fra = NOMADcat.reclen == NOMAD1reclen ? ra3 : ra4 ;

    /* We look whether we're not just at the correct location ??? */
    if ((NOMADcat.dat.pos>=0) && (NOMADcat.dat.ppos>=0)) {
	a0 = (*fra)(NOMADcat.abuf + NOMADcat.dat.ppos);
	a1 = (*fra)(NOMADcat.abuf + NOMADcat.dat.pos);
	if ((da > a0) && (da <= a1)) 
	    return(NOMADcat.dat.id);
	/* Go back to the beginning of the chunk (done in read_chunk) */
	if (read_chunk(NOMADcat.dat.RAz)) return(-1);
    }

    /* Find the RA position in the index  (o, ID, RA)
       RA(n) <= value < RA(n+1)) 
    */
    aval = binloc(NOMADcat.dat.index+2, NOMADcat.dat.endex, value, 3) -2;
    if (!aval) aval = NOMADcat.dat.index;

    /* Move in the chunk */
    while (NOMADcat.dat.pos < NOMADcat.dat.len) {
	a1 = (*fra)(NOMADcat.abuf + NOMADcat.dat.pos);
	if (a1 >= da) return(NOMADcat.dat.id);
	if (nomad_next()<=0) break;
    }
    return(0);		/* Impossible... */
}

int4 nomad_seek(int4 ra, int4 sd)
/*++++++++++++++++
.PURPOSE  Seek to specified position (in mas)
.RETURNS  The corresponding NOMAD number (relative to SD zone)
-----------------*/
{
  int4 id;
  int z;

    if (!NOMADcat.root) init((char *)0) ;
    z = sd/Dstep ;
    if ((z<0) || (z >= 1800))
        return(-1) ; 		/* Outside Boundaries */

    /* Maye we're at the right position ? */
    if (NOMADcat.SDz == z) {
	id = aseek(ra);
	if (id > 0) return(id);
    }

    /* We can't find in the loaded zone. Reload. */
    if (nomad_zopen(z)) return(-1) ;

    /* Seek now until the specified RA found */
    return(aseek(ra));
}

static int set_id(int4 id)
/*++++++++++++++++
.PURPOSE  Set the current record to id 
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
.REMARKS  We assume we're in the correct SDzone.
-----------------*/
{
  int4 *aval, di;
  int i; 

    /* Maybe we're already at the correct location ? */
    if (NOMADcat.dat.id == id) return(0);

    /* A zero (negative) ID sets to beginning of file */
    if (id <= 0) id = 1;

    /* Locate first which chunk is concerned -- i.e. search in listofChunks */
    aval = binloc(NOMADcat.chunks+1, NOMADcat.echunks, id, 2);
    if (!aval) return(-1);		/* e.g. id<0 */
    aval--;				/* aval -> offset */
    i = (aval-NOMADcat.chunks)>>1;	/* 2 numbers (o,ID) per chunk */
    if (i >= NOMADcat.nchunks) return(-1) ;

    /* Maybe we're already in the correct chunk -- taken care by read_chunk */
    if (i == NOMADcat.dat.RAz) {
	di = id - NOMADcat.dat.id;
	/* When new ID close to an index boundary, better to use the Index */
	if ((di>(Dacc/3)) && ((id-NOMADcat.dat.id0)%Dacc) < (Dacc/3))
	  di = -1;
    }
    else di = -1;

    /* When asking for just a bit further, use successive read's */
    if (di<0) {
        if (read_chunk(i) < 0) return(-1);
        /* Locate in the index of the Chunk */
        aval = binloc(NOMADcat.dat.index+1, NOMADcat.dat.endex, id, 3);
        NOMADcat.dat.id  = *aval;
        NOMADcat.dat.pos = *--aval;
        NOMADcat.dat.ppos= -1;		/* Not known... */
    }

    /* Loop the in <~1000 (Dacc) records */
    while (NOMADcat.dat.id < id) {
	if (nomad_next()<= 0) return(-1);
    }
    return(0) ;
}

int nomad_set(int zone, int4 id)
/*++++++++++++++++
.PURPOSE  Position in specified zone + id
.RETURNS  0 / -1
-----------------*/
{
    if (!NOMADcat.root) init((char *)0) ;
    /* Go to specified zone */
    if (zone == NOMADcat.SDz) ;
    else if (nomad_zopen(zone)) return(-1);

    return(set_id(id));
}

/*==================================================================
		Convert the Input Record(s)
 *==================================================================*/
static int extra_value(int value, int max, int max2)
/*++++++++++++++++
.PURPOSE  Retrieve a value from Index
.RETURNS  The Value
-----------------*/
{
    if (value <= max) return(value);
    return( value > max2 ?
        NOMADcat.dat.xtra4[value-max2-1] :
    	NOMADcat.dat.xtra2[value-max -1]
    );
}

static int ed_rec(NOMADrec *prec)
/*++++++++++++++++
.PURPOSE  Convert the current compressed NOMAD record into its standard 
	  structure
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
-----------------*/
{
  unsigned char *pc, presence;
  int status, m, i;

    memset(prec, 0, sizeof(NOMADrec));

    /* Convert the compacted record */
    pc = (NOMADcat.abuf + NOMADcat.dat.pos) ;
    presence = *pc;			/* presence of mags and IDs */
    prec->zone = NOMADcat.SDz;
    prec->Uzone= 0;
    prec->id   = NOMADcat.dat.id;

    /* USNO-B Name */
    prec->Uid = (pc[1]<<16)|(pc[2]<<8)|pc[3];
    if (prec->Uid) {
	prec->Uzone = prec->zone;
	if (prec->Uid&0x800000) prec->Uzone -= 1;
	if (prec->Uid&0x400000) prec->Uzone += 1;
	prec->Uid &= 0x3fffff;
    }

    prec->ra   = (pc[4]<<16) | (pc[5]<<8) | pc[6];
    if (NOMADcat.reclen > NOMAD1reclen) {
	++pc;
	prec->ra <<= 8;
	prec->ra |= pc[6]; 
    }
    prec->ra  += NOMADcat.dat.ra0;

    prec->e_ra = get(pc+7, 0, 10);
    prec->epra = get(pc+8, 2, 10) + NOMADcat.oval[2];

    prec->sd   = get(pc+9, 4, 20);
    prec->sd  += NOMADcat.dat.sd0;

    prec->e_sd = get(pc+12, 0, 10);
    prec->epsd = get(pc+13, 2, 10) + NOMADcat.oval[2];

    prec->pmra = extra_value(get(pc+14, 4, 14), EXTRA_14_2, EXTRA_14_4)
	       + NOMADcat.oval[0];
    prec->pmsd = extra_value(get(pc+16, 2, 14), EXTRA_14_2, EXTRA_14_4)
	       + NOMADcat.oval[0];
    prec->e_pmra = extra_value(get(pc+18, 0, 13), EXTRA_13_2, EXTRA_13_4);
    prec->e_pmsd = extra_value(get(pc+19, 5, 13), EXTRA_13_2, EXTRA_13_4);

    status = ((pc[21]&0x3f)<<24) | (pc[22]<<16) | (pc[23]<<8) | pc[24];
    the_status = status;  /* printf("=%02X", status); */
    prec->flags= status>>12;

    pc += NOMAD1reclen;		/* Fixed-length part done. */

    if (presence&0x80) {	/* UCAC2 */
	prec->UCAC2 = (pc[0]<<24)|(pc[1]<<16)|(pc[2]<<8) | pc[3];
	pc += 4;
    }

    if (presence&0x40) {	/* Tycho */
	prec->TYCfar = pc[0]>>7;
	prec->TYC3 = get(pc, 29, 3);
	if (prec->TYC3 == 0) 
	    prec->HIP = (get(pc, 1, 28))%1000000;
	else {
	    prec->TYC1 = get(pc, 1, 14);
	    prec->TYC2 = get(pc, 15, 14);
	}
	pc += 4;
    }

    /* Magnitudes */
    for (i=0, m=0x20; i<6; i++, m >>= 1) {
	if (presence&m) {
	    prec->mag[i] = ((pc[0]<<8)|pc[1]) + NOMADcat.oval[1];
	    pc += 2;
	}
	else prec->mag[i] = 0x8000 ;	/* -32768 */
    }

    /* Sources */
    prec->abvr[0] = status&7;
    prec->abvr[1] = (status>>3)&7;	/* Blue photometry */
    prec->abvr[2] = (status>>6)&7;	/* V    photometry */
    prec->abvr[3] = (status>>9)&7;	/* Red  photometry */
    if (prec->flags&NOMAD_OMAGBIT) prec->abvr[1] = 8;
    if (prec->flags&NOMAD_EMAGBIT) prec->abvr[3] = 9;

    prec->pmtot= -1;	/* Computed if required only */
    prec->rho  = -1;

    return(0) ;
}

int nomad_read(NOMADrec *rec)
/*++++++++++++++++
.PURPOSE  Read next record in NOMAD
.RETURNS  0(End) / -1 (Err) / 1(OK)
.REMARKS  When end if a RA file is found, action depends on stopat_eof
	When rec is NULL, jusr set to to next.
-----------------*/
{
  int i; 
    if (!NOMADcat.root) init((char *)0) ;
    if ((!NOMADcat.abuf) || (NOMADcat.SDz < 0)) {
	i = NOMADcat.SDz;
	if (i < 0) i=0, fprintf(stderr, 
	    "++++nomad_read: no zone specified, read from South Pole\n");
	if (nomad_set(i, 0) < 0) return(-1);
    }
    if (NOMADcat.dat.pos >= NOMADcat.dat.len) {	/* End of Chunk ? */
	if (stopat_eof == 2) return(0) ;
	i = NOMADcat.dat.RAz + 1;		/* Next Chunk */
	if (i >= NOMADcat.nchunks) {		/* Need to read a new file */
	    if (stopat_eof) return(0);
	    i = NOMADcat.SDz+1;			/* New Zone...*/
	    if (i >= 1800) return(0);		/* FINAL END  */
	    if (nomad_zopen(i)) return(-1);
	    i = 0;
	}
	if (read_chunk(i)) return(-1);
    }
    if (rec) ed_rec(rec) ;
    nomad_next();
    return(1) ;
}

int nomad_get(int zone, int4 id, NOMADrec *rec)
/*++++++++++++++++
.PURPOSE  Read specified record in NOMAD
.RETURNS  0 / -1
-----------------*/
{
    if (nomad_set(zone, id) < 0) return(-1) ;
    return(nomad_read(rec)) ;
}

static int edra(char *buf, int4 ra, int opt /* 0=deg, 1=sexa 16=mas */)
/*++++++++++++++++
.PURPOSE  Edit a RA (unit=1mas) into HH:MM:SS.SSSS
.RETURNS  Number of bytes (13 / 11 / 10)
-----------------*/
{
  int i ; 
  int4 value;
    if (opt&16) {	/* Edit in mas */
	i = 10 ;
	buf[--i] = '0' + (ra%10);
	for (value=ra/10; value; value /= 10) buf[--i] = '0' + (value%10);
	while (i>0) buf[--i] = ' ';
	return(10);
    }
    if (opt) {		/* Sexagesimal -- round it */
    	value = ra;
	value = value - (value/3);	/* 2/3 for number in 0.1ms */
    	for (i=12; i>8; i--) buf[i] = '0' + value%10, value /= 10 ;
    	buf[8] = '.' ;
    	buf[7] = '0' + value%10, value /= 10 ;
    	buf[6] = '0' + value%6,  value /= 6 ;
    	buf[5] = ':' ;
    	buf[4] = '0' + value%10, value /= 10 ;
    	buf[3] = '0' + value%6,  value /= 6 ;
    	buf[2] = ':' ;
    	buf[1] = '0' + value%10, value /= 10 ;
    	buf[0] = '0' + value ;
    	return(13) ;
    }
    else {		/* Caution with 32-bit limit !! */
	value = 3*(ra%3600000);		/* Degree fraction */
	value -= ((2*value+13)/27);	/* 10^7/3600000 = 3 * (1 - 2/27) */
    	for (i=10; i>3; i--) buf[i] = '0' + value%10, value /= 10 ;
    	buf[3] = '.' ;
	/* Degrees */
	value = ra/3600000;
    	buf[2] = '0' + value%10, value /= 10 ;
    	buf[1] = '0' + value%10, value /= 10 ;
    	buf[0] = '0' + value ;
    	return(11) ;
    }
}

static int edsd(char *buf, int4 sd, int opt /* 0=deg, 1=sexa 16=mas */)
/*++++++++++++++++
.PURPOSE  Edit a SPD (unit=1mas) into +DD:MM:SS.SSS
.RETURNS  Number of bytes (13 / 11 / 10)
-----------------*/
{
  int i ; 
  int4 value ;
    if (opt&16) {	/* Edit in mas */
	i = 10;
	buf[--i] = '0' + (sd%10);
	for (value=sd/10; value; value /= 10) buf[--i] = '0' + (value%10);
	while (i>0) buf[--i] = ' ';
	return(10);
    }
    value = sd - 90*3600000;
    if (value > 0) buf[0] = '+' ;
    else buf[0] = '-', value = -value ;
    if (opt) {		/* Sexagesimal */
	buf[12] = '0' + value%10, value /= 10 ;
    	buf[11] = '0' + value%10, value /= 10 ;
    	buf[10] = '0' + value%10, value /= 10 ;
    	buf[9] = '.' ;
    	buf[8] = '0' + value%10, value /= 10 ;
    	buf[7] = '0' + value%6,  value /= 6 ;
    	buf[6] = ':' ;
    	buf[5] = '0' + value%10, value /= 10 ;
    	buf[4] = '0' + value%6,  value /= 6 ;
    	buf[3] = ':' ;
    	buf[2] = '0' + value%10, value /= 10 ;
    	buf[1] = '0' + value ;
    	return(13) ;
    }
    else {
    	value = 3*value ; 		/* 10^6/360000 = 3 * (1 - 2/27) */
    	value -= ((2*value+13)/27) ; 
    	for (i=10; i>3; i--) buf[i] = '0' + value%10, value /= 10 ;
    	buf[3] = '.' ;
    	buf[2] = '0' + value%10, value /= 10 ;
    	buf[1] = '0' + value ;
    	return(11) ;
    }
}

static int edmag(char *buf, int mag) 
/*++++++++++++++++
.PURPOSE  Edit a Magnitude in mmag
.RETURNS  Number of bytes (6)
.REMARKS  Value -32768 = unknown
-----------------*/
{
  int value ;
  value = mag ;
    buf[0] = buf[1] = buf[2] = buf[3] = buf[4] = buf[5] = ' ' ;
    if (value <= -32768) 
	buf[1] = buf[2] = buf[3] = '-'; 
    else {
        if (value < 0) value = -value, buf[0] = '-';
	buf[5] = '0' + value%10, value /= 10 ;
        buf[4] = '0' + value%10, value /= 10 ;
        buf[3] = '0' + value%10, value /= 10 ;
        buf[2] = '.' ;
        buf[1] = '0' + value%10, value /= 10 ;
        if (value) buf[0] = '0' + value ;
    }
    return(6) ;
}

static int edmas(char *buf, int4 value)
/*++++++++++++++++
.PURPOSE  Edit a value in units of mas as arcsec
.RETURNS  Number of bytes written. Round to 0.01arcsec
-----------------*/
{
  static char aval[16] ;
  char s, *p, *b ;  int4 x ;
    if (value < 0) s = '-', x = (5-value)/10;
    else           s = 0  , x = (6+value)/10;
    p = aval + sizeof(aval) ;
    *--p = 0 ;
    *--p = '0' + x%10 ; x /= 10 ;
    *--p = '0' + x%10 ; x /= 10 ;
    *--p = '.' ;
    do { *--p = '0' + x%10, x /= 10 ; }
    while (x) ;
    if (s) *--p = s ;

    for (b = buf, x = 9 - strlen(p) ; --x >= 0; b++) *b = ' ' ;
    while (*p) *(b++) = *(p++) ;
    *b = 0 ;
    return(b-buf) ;
}

/*==================================================================
		Edit a Record
 *==================================================================*/

char *nomad_head(int opt 
	/*0=deg, 1=sexa, 2=ID, 4=Ep, 8=x,y, 16=mas 32=xy 128=pos.only */)
/*++++++++++++++++
.PURPOSE  Get the Header for one edition
.RETURNS  The Header, starting by #
-----------------*/
{
  static char buf[500] ;
  char *p; int i, b ;

    p = buf ; 

  if (opt&128) ; else {		/* V1.2 */

    if (opt&2) {	/* Edit the NOMAD-ID */
    	sprintf(p, "#%-11s|", "NOMAD1.0");
    	p += strlen(p) ;
    }
    /* else putchar('#'); */

    /* USNO-B Number */
    if (opt&8) {
	sprintf(p, "%-12s ", "USNO-B1.0");
    	p += strlen(p) ;
    }
    else if (p == buf) *(p++) = '#';
    else *(p++) = 'B';

    /* UCAC2 Number */
    if (opt&8) {
	sprintf(p, p==buf ? "#%7s ":"%8s ", "UCAC2");
	p += strlen(p);
    }
    else if (p == buf) *(p++) = '#';
    else *(p++) = 'C';

    /* Tycho Number */
    if (opt&8) {
	sprintf(p, "%12s ", "Tycho-2/HIP");
	p += strlen(p);
    }
    else *(p++) = 'T';
    *(p++) = 'Y';
    *(p++) = 'M';
    *(p++) = '|';
  }

    /* Edit the Position */
    i = edra(p, 0, opt&17) + edsd(p, 0, opt&17);
    b = (i-15)/2;
    if (p == buf)   { *(p++) = '#'; b--; }	/* V1.2 */
    while(--b >= 0) { *(p++) = ' ', i--; }
    if (opt&16)
	 strcpy(p, "RA   (mas)  SPD");
    else strcpy(p, "RA  (J2000) Dec");
    i -= 15, p += 15;
    while(--i >= 0) *(p++) = ' ';
    if (opt&128) {		/* That's all for position-only */
	*p = 0;
        return(buf) ;
    }
    *(p++) = ' ';
    *(p++) = 'r';
    *(p++) = ' ';

    /* Edit the sigma's */
    strcpy(p, " sRA  sDE|"); p += strlen(p);

    /* Epoch */
    if (opt&4) {
	sprintf(p, "%6s %6s|", "Ep.RA", "EpDec");
        p += strlen(p);
    }

    /* Proper Motions + Sigmas */
    sprintf(p, "%8s %8s %5s %5s|", "pmRA", "pmDE", "spRA", "spDE");
    p += strlen(p);

    /* Magnitudes */
    for (i=0; i<6; i++) {
	*(p++) = ' ';
	*(p++) = NOMADmag[i];
	*(p++) = 'm';
	*(p++) = 'a';
	*(p++) = 'g';
	if (i<3) { *(p++) = '.'; *(p++) = 'r'; }
	else *(p++) = ' ';
	*(p++) = (i%3) == 2 ? '|' : ' ';
    }

    /* Flags */
    if (opt&8) {
	sprintf(p, "%5s", "Flags"); p += strlen(p);
    }
    *(p++) = 'R';
    *(p++) = '|';

    /* Distance in arcsec */
    if (opt&32) sprintf(p, " ;%9s %9s", "x(\")", "y(\")");
    else sprintf(p, " ;%9s", "r(\")");
    p += strlen(p);

    *p = 0;
    return(buf) ;
}

char *nomad2a(NOMADrec *pr, 
   int opt /*0=deg, 1=sexa, 2=ID, 4=Epoch, 8=allIDs 16=mas 32=xy 128=pos */)
/*++++++++++++++++
.PURPOSE  Edit (in a static buffer) one record.
.RETURNS  The edited record
-----------------*/
{
  static char buf[500] ;
  char *p ; int i ;
  int4 value ;

    p = buf ;

  if (opt&128) ; else {		/* V1.2 */
    if (opt&2) {	/* Edit the NOMAD-ID */
    	sprintf(p, "%04d-%07d|", pr->zone, pr->id) ;
    	p += strlen(p) ;
    }
    /* else putchar(' '); */

    /* USNO-B Number */
    if (opt&8) {
	if (pr->Uid ||  pr->Uzone)
             sprintf(p, "%04d-%07d ", pr->Uzone, pr->Uid);
	else sprintf(p, "%4s-%-7s ",  "-", "-");
	p += strlen(p);
    }
    else *(p++) = pr->Uid ? 'B' : '.';

    /* UCAC2 Number */
    if (opt&8) {
	if (pr->UCAC2) 
	     sprintf(p, "%08d ", pr->UCAC2);
	else sprintf(p, "%8d ", 0);
	p += strlen(p);
    }
    else *(p++) = pr->UCAC2 ? 'C' : '.';

    /* Tycho Number */
    if (opt&8) {
	if (pr->TYC3) sprintf(p, 
		"%04d-%05d-%d", pr->TYC1, pr->TYC2, pr->TYC3);
	else if (pr->HIP) 
	     sprintf(p, "    %06d  ", pr->HIP);
	else sprintf(p, "%10d  ", 0);
        p += strlen(p);
	*(p++) = pr->TYCfar>=2 ? 	/* Large distance */
	    '*' : ' ';
    }
    else if (pr->TYC3) *(p++) = 'T';
    else if (pr->HIP)  *(p++) = 'H';
    else *(p++) = '.';

    /* Flags 2MASS, YB6 */
    *(p++) = pr->flags&NOMAD_YB6   ? 'Y' : '.';
    *(p++) = pr->flags&NOMAD_2MASS ? 'M' : '.';
    *(p++) = '|';
  }

    /* Edit the Position */
    p += edra(p, pr->ra, opt&17) ;
    p += edsd(p, pr->sd, opt&17) ;
    if (opt&128) {		/* That's all for position-only */
	*p = 0;
        return(buf) ;
    }
    *(p++) = ' ';
    *(p++) = source_letter[pr->abvr[0]];
    *(p++) = ' ';

    /* Edit the sigma's */
    sprintf(p, "%4d %4d|", pr->e_ra, pr->e_sd);
    p += strlen(p) ;

    /* Epoch */
    if (opt&4) {
	sprintf(p, "%d.%d ", pr->epra/10, pr->epra%10);
	p += strlen(p) ;
	sprintf(p, "%d.%d|", pr->epsd/10, pr->epsd%10);
	p += strlen(p) ;
    }

    /* Proper Motions + Sigmas -- Note that pmsd = pmdec */
    sprintf(p, "%+7d", pr->pmra);  p[7] = p[6]; p[6] = '.'; p += 8;
    *(p++) = ' ';
    sprintf(p, "%+7d", pr->pmsd);  p[7] = p[6]; p[6] = '.'; p += 8;
    *(p++) = ' ';
    value = pr->e_pmra; if(value>9999) value=9999;
    sprintf(p, "%4d", value); p[4] = p[3]; p[3] = '.'; p += 5;
    *(p++) = ' ';
    value = pr->e_pmsd; if(value>9999) value=9999;
    sprintf(p, "%4d", value); p[4] = p[3]; p[3] = '.'; p += 5;
    *(p++) = '|';

    /* Magnitudes */
    for (i=0; i<6; i++) {
	p += edmag(p, pr->mag[i]);
	if (i<3) *(p++) = source_letter[pr->abvr[i+1]];
	*(p++) = (i%3) == 2 ? '|' : ' ';
    }

    /* Flags */
    if (opt&8) {
	sprintf(p, "%05x", pr->flags&(~(NOMAD_YB6|NOMAD_2MASS)));
	p += strlen(p);
    }
    /* Special flags: * when problem, R for recommended standard */
    if (pr->flags&(NOMAD_SPIKE|NOMAD_TYCONF|NOMAD_BSCONF|NOMAD_BSART
		  |NOMAD_XRBIT|NOMAD_ITMBIT|NOMAD_IUCBIT|NOMAD_ITYBIT)) 
	 *(p++) = '*';
    else if (pr->flags&NOMAD_USEME) 
	 *(p++) = 'R';
    else *(p++) = ' ';
    *(p++) = '|';

    /* Distance in arcsec */
    if (pr->rho >= 0)  {
	*(p++) = ' ' ; *(p++) = ';' ; *(p++) = ' ' ; 
	if (opt&32) {	/* Edit x,y */
	    p += edmas(p, pr->xy[0]);
	    *(p++) = ' '  ;
	    p += edmas(p, pr->xy[1]);
	}
	else p += edmas(p, pr->rho) ;
    }

    *p = 0 ;
    return(buf) ;
}

static int prec(NOMADrec *rec) { 	/* DEFAULT ACTION ROUTINE */
    printf("%s\n", nomad2a(rec, 3)); 
    return(1); 
}

/*==================================================================
		Operate on a file
 *==================================================================*/

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 int s2dist(int *a1, int *a2)
/*++++++++++++++++
.PURPOSE  Sort routine (called by qsort), 5 integers (min,max(sd) minx,max(ra)
		distance)
.RETURNS  <0 0 >0
-----------------*/
{
  int diff ;
    diff = a1[4] - a2[4] ;
    if (!diff) diff = a1[0] - a2[0] ;
    return(diff) ;
}

int4 nomad_search(int4 ra[2], int4 sd[2], int (*digest_routine)())
/*++++++++++++++++
.PURPOSE  Search USNO stars within range of positions.
.RETURNS  Number of tested records.
.REMARKS  The supplied digest_routine(NOMADrec *record) does whatever,
	return -1 if stop asked. Default is just to print out.
	The search is done from middle to borders, according to
	a crude distance estimate between the center and the point
  ===>	When nomad_options&0x10 ==> Keep Original Order  <===
-----------------*/
{
  static NOMADrec rec ;
  double u[3], u0[3], dsin, dx, du;
  int4 lac[80*2+1], ldc[1800*2+1];	/* Limits on one zone */
  int4 ld[2], la[2], d, a, acen, dcen, high;
  int4 dist2, tested_records ;
  int ja, jd, i, *alist, *p;
  int goon, saved_eof ;

    if (!NOMADcat.root) init((char *)0) ;
    if (!digest_routine) digest_routine = prec ;
    saved_eof = stopat_eof ; 
    stopat_eof = 2 ;			/* Stop after each chunk */
    tested_records = 0 ;

    /* Verify whole range */
    if (ra) {
	if ((ra[0] == 0) && (ra[1] >= 360*3600*1000-1))
	   ra = (int4 *)0;
    }
    if (sd) {
	if ((sd[0] == 0) && (sd[1] >= 180*3600*1000-1)) 
	   sd = (int4 *)0;
    }

    if (!ra) {
	if (!sd) {	/* Whole Sky */
	    if (NOMADcat.fno >= 0) nomad_close();
	    stopat_eof = 0;
	    while (nomad_read(&rec)>0) {
		tested_records++;
		if ((*digest_routine)(&rec) < 0)        /* Ask to stop! */
		    break;
	    }
	    stopat_eof = saved_eof ;
	    return(tested_records) ;
	}
	stopat_eof = 1;	/* No test on RA ==> stop at end-of-file */
	ra = la;
	ra[0] = 0; ra[1] = (360*3600*1000)-1;
    }
    if (!sd) {
        sd = ld ;
	sd[0] = 0;
	sd[1] = 180*3600*1000;
    }

    /* Set the limits in their legal range/order. 
       Remember that the order of RA limits is important !
    */
    ld[0] = sd[0] ;
    ld[1] = sd[1] ;
    if (ld[0] > ld[1]) d=ld[1], ld[1]=ld[0], ld[0]=d ;
    la[0] = ra[0] % (360*3600*1000);
    la[1] = ra[1] % (360*3600*1000);

    /* Direction Cosines of the Center */
    /* WARNING: la[0] close to the 2**31 limit, a[0]+a[1] may be over limit! */
    acen = (la[0]/2) + (la[1]/2) + ((la[0]&1)+(la[1]&1))/2;  
    if(la[0]>la[1]) acen += 180*3600*1000;
    if(acen >= 360*1000*3600) acen -= 360*1000*3600;
    dcen = (ld[0] + ld[1])/2; 
    dsin  = sind(dcen/3.6e6);
    u0[0] = dsin*cosd(acen/3.6e6);
    u0[1] = dsin*sind(acen/3.6e6);
    u0[2] = cosd(dcen/3.6e6);

    /* Limits on the Zones */
    ldc[0] = ld[0];
    jd = 1;
    high = (Dstep*(ldc[0]/Dstep)) + (Dstep-1);
    while (ld[1] > (high+1)) {	/* The pole requires the +1 */
	ldc[jd++] = high;
	ldc[jd++] = high+1;
	high += Dstep;
    }
    ldc[jd++] = ld[1];
    ldc[jd] = -1;		/* End indicator */
    /* jd = 2*Number of zones   */

    /* Distribute RA Range over the 80 chunks */
    lac[0] = la[0];
    ja = 1;
    if (stopat_eof == 2) {
        a = (lac[0] >> 24) ;		/* Chunk index */
        high = (a<<24) | 0xffffff;	/* Upper Limit (mas) */
        if (la[0] > la[1]) {		/* Around 0deg (mas) */
	    while (high < 360*3600*1000) {
	        lac[ja++] = high;
		lac[ja++] = high+1;
	        high += 0x01000000;
	    }
	    lac[ja++] = 360*3600*1000;	/* Upper limit (0h) */
	    lac[ja++] = 0;
	    high = 0xffffff;
        }
        while (la[1] > high) {
	    lac[ja++] = high;
	    lac[ja++] = high+1;
	    high += 0x01000000;
        }
    }
    lac[ja++] = la[1] ;
    lac[ja] = -1;			/* Indicates the End */
    /* ja = 2*Number of RA zones */

    /* Create a table with all the chunks to read.	*/
    alist = (int *)malloc(((ja*jd/4) + 1)*5*sizeof(int));
    for (i=0, d=0; d<jd; d+=2) {
	dcen = (ldc[d]+ldc[d+1]+1)/2;
	u[2] = cosd(dcen/3.6e6);
	dsin = sind(dcen/3.6e6);
	for (a=0; a<ja; a+=2) {
	    alist[i++] = ldc[d]; alist[i++] = ldc[d+1];
	    alist[i++] = lac[a]; alist[i++] = lac[a+1];
	    if (nomad_options&0x10) dist2 = 0;
	    else {
		acen = (lac[a]/2)+(lac[a+1]+1)/2;
		u[0] = dsin*cosd(acen/3.6e6);
		u[1] = dsin*sind(acen/3.6e6);
		dx = u[0]-u0[0];
		du = dx*dx;
		dx = u[1]-u0[1];
		du += dx*dx;
		dx = u[2]-u0[2];
		du += dx*dx;
		dist2 = (du*4.e8) ;
	    }
	    alist[i++] = dist2;
        }
    }
    /* End Indicator */
    alist[i++] = -1; alist[i++] = -1;
    alist[i++] = -1; alist[i++] = -1;
    alist[i++] = -1;

    if (nomad_options&0x10) ;	/* Keep original order */
    else qsort(alist, (ja*jd)/4, 5*sizeof(int), (SORT_FCT)s2dist);

    /* Loop on the zones */
    for (goon=0, p=alist; (goon>=0) && (p[4]>=0); p += 5) {
	/* Ignore the large distances  written when reclen      */
	if (p[4] >= 1999999999) continue;
	if (nomad_seek(p[2], p[0]) < 0) break;
	while((goon>=0) && (nomad_read(&rec) > 0)) {
	    tested_records++;
	    if (rec.ra > p[3]) break;
	    if (rec.sd < p[0]) continue;
	    if (rec.sd > p[1]) continue;
	    goon = (*digest_routine)(&rec);
	    /* goon=1 for good, 0 for rejected, -1 for stop	*/
	}
    }

    free(alist);
    stopat_eof = saved_eof ;
    return(tested_records) ;
}

/*==================================================================
		Main Program
 *==================================================================*/
#ifdef TEST
static char usage[] = "\
Usage: nomad_sub [-c ori_file] compressed_input_zone [-q]\n\
";

static char help[] = "\
		      -c: Compare with original NOMAD file\n\
   compressed_input_zone: one of the files created by nomad_make\n\
                      -q: query (ask for a position)\n\
";

static int4 NOMADid;

#define CMP1(t,n,o)  if(o != n) printf(\
		"++++Rec.#%07d: %s %d != %d (%d) %08x %08x\n", \
		NOMADid, t, n, o, n-o, n, o) 

int nomad_cmp(int4 *rec, NOMADrec *prec)
/*++++++++++++++++
.PURPOSE  Compare with original
.RETURNS  0 (OK) 
.REMARKS  We assume that "zone" and "id" are already filled.
   	Histograms computed.
-----------------*/
{
  static char *Tmag[] = { 
      "B.mag", "V.mag", "R.mag", "J.mag", "H.mag", "K.mag" 
  };
  int i;
    
    /* Swap the 22-long-int record */
    swapping ^= 2; swap4(rec, 22); swapping ^= 2;
    NOMADid = prec->id;

    /* Install RA, SD ... */
#if 0
    prec->ra = rec[0]; 
    prec->sd = rec[1]; 
    prec->e_ra = rec[2];
    prec->e_sd = rec[3];
    prec->pmra = rec[4];
    prec->pmsd = rec[5];
    prec->e_pmra = rec[6];
    prec->e_pmsd = rec[7];
    prec->epra = (rec[8]+50)/100;
    prec->epsd = (rec[9]+50)/100;
    for (i=0; i<6; i++) prec->mag[i] = rec[10+i];
    if (rec[17]) prec->flags |= NOMAD_2MASS;
    if (rec[18]) prec->flags |= NOMAD_YB6;
    prec->UCAC2 = rec[19];
    /* [20]=Tycho-2 */
    prec->flags = rec[21]>>12;
#else
    CMP1("RAmas", prec->ra, rec[0]);
    CMP1("SDmas", prec->sd, rec[1]);
    CMP1("e_RAm", prec->e_ra, rec[2]);
    CMP1("e_SDm", prec->e_sd, rec[3]);
    CMP1("pmRA.", prec->pmra, rec[4]);
    CMP1("pmSD.", prec->pmsd, rec[5]);
    CMP1("epmRA", prec->e_pmra, rec[6]);
    CMP1("epmSD", prec->e_pmsd, rec[7]);
    rec[8] = (rec[8]+50)/100;
    rec[9] = (rec[9]+50)/100;
    CMP1("ep.RA", prec->epra, rec[8]);
    CMP1("ep.SD", prec->epsd, rec[9]);
    for (i=0; i<6; i++) {
	if (rec[10+i] >= 30000) rec[10+i] = -32768;
	CMP1(Tmag[i], prec->mag[i], rec[10+i]);
    }
    CMP1("UCAC2", prec->UCAC2, rec[19]);
    rec[21] >>= 12;
    if (rec[17]) rec[21] |= NOMAD_2MASS;
    if (rec[18]) rec[21] |= NOMAD_YB6;
    CMP1("Flags", prec->flags, rec[21]);
#endif

    return(0);
}

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

int main (argc, argv) int argc; char **argv;
{
  int4 rec22[22];
  FILE *bf; char *bname;
  FILE *of; char *oname;
  char buf[BUFSIZ], *p ;
  int4 id = 0, ra, sd ;
  double q[2];
  int  getpos = 0 ;
  int  z, len, records;
  NOMADrec rec, brec;
  char read_zone = 0;

    nomad_options = 1 ;
    bf = (FILE *)0; bname = (char *)0;
    of = (FILE *)0; oname = (char *)0;
    while (argc > 1) {
	p = *++argv; --argc;
	if (*p == '-') switch(p[1]) {
	  case 'c':	/* Compare with original */
	    p = *++argv; --argc;
	    of = fopen(p, "r");
	    if (!of) {
		fprintf(stderr, "****");
		perror(p);
		exit(1);
	    }
	    continue;
	  case 'q': 	/* Query */
	    getpos = 2;
	    continue;
	  default:	/* Bin file */
	    fprintf(stderr, "%s%s", usage, help);
	    exit(1);
	}
	/* A binary file */
	if (nomad_fopen(p))
	    exit(1);
	continue;
    }

    nomad_stop(1);		/* Stop at EOF */
    records = 0;
    while (1) {
	buf[0] = 0;
	if (of) {
	    if (!fread(rec22, 4, 22, of)) break;
	    nomad_read(&rec);
	    nomad_cmp(rec22, &rec);
	    continue;
	}
	if (bf) fread(&brec, sizeof(brec),1,bf);
    	else if (isatty(0) && (read_zone == 0)) { 
	    if (getpos) 
		 printf("----Give a position (degrees): ") ;
    	    else printf("----Give a zone-NOMAD  number: ") ;
    	    if (! fgets(buf, sizeof(buf), stdin)) break ;
        } else buf[0] = 0 ;
	if (isprint(buf[0])) {	/* I've a zone */
	    for (p=buf; isspace(*p); p++) ;
	    q[0] = atof(p) ;
	    while (isdigit(*p)) p++ ;
	    while (isspace(*p)) p++ ;
	    if (*p == '-') p++;
	    while (isspace(*p)) p++ ;
	    q[1] = atof(p) ;
	    if (getpos) {
	    	ra = q[0]*3.6e6 ; sd = 3.6e6*(90. + q[1]) ;
	    	id = nomad_seek(ra, sd) ;
	    } 
	    else {
		z = q[0] ; id = q[1];
		if (nomad_set(z, id) < 0) continue;
		read_zone = id == 0;
	    }
	}
    	len = nomad_read(&rec);
	if (len>0) records++,
	    printf("%s\n%s\n", nomad_head(7|8), nomad2a(&rec, 7|8)),
	    printf("%s\n%s\n", nomad_head(6|8), nomad2a(&rec, 6|8)),
	    printf("%s\n%s\n", nomad_head(16|8), nomad2a(&rec, 16|8)) ; 
	else { printf("====EOF====\n"); break; }
	if (bf) { 		/* Compare... */
	  int i, *a, *b;
	    if (memcmp(&rec, &brec, sizeof(brec))) {
    		printf("=%02X ", the_status);
    		printf("++ %04d-%d, ra=%010d, sd=%09d ",
	    	   rec.zone, rec.id, rec.ra, rec.sd);
		printf("++ DIFFER!\n");
		a = (int *)&rec; b=(int *)&brec;
		for (i=0; i<sizeof(rec); i+=4, a++, b++) {
		    if (*a != *b) printf("    +%3d: %08X %08X\n", i, *a, *b);
		}
	    }
	}
    	else printf("== %04d-%d, ra=%010d, sd=%09d\n",
	    rec.zone, rec.id, rec.ra, rec.sd);
    }
    printf("====Total of records read: %d\n", records);
    nomad_close() ;
    exit(0);
}
#endif