/*++++++++++++++
.IDENTIFICATION pmmsub.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-A2.0 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   01-Mar-1997
.VERSION  1.1   01-Apr-1997: Error in pmm_search ....
.VERSION  2.0   20-Oct-1998: USNO-A2.0
.VERSION  2.1   12-Nov-2000: Use Byte Swaps (for PCs)
.VERSION  2.2   03-Mar-2001: Reorder the blocks for more efficiency
.VERSION  2.3   06-Mar-2001: ***ERROR detected by P. Fernique***
			(Syndrome: when center close to a 7.5degrees limit in RA
			a part of the stars is missing)
.COMMENTS       Edition of PMM data.

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

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

#include <stdio.h>
#include <stdlib.h>	/* Malloc   */
#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 PMM queries.......	*/
int pmm_options ;		/* 0x10 for OrigOrder	*/

static struct {
    char *root;			/* Root name of PMM	*/
    char *name ;		/* Name of opened file	*/
    short SDz, RAz ;		/* Zone opened		*/
    int fno;			/* Opened file number 	*/
    long obuf;			/* Byte position loaded	*/
    int mbuf, nbuf, ibuf ;	/* Max, used, offset	*/
    char *abuf;			/* Loaded data		*/
    char *line ;		/* First line of file	*/
    HEADER *h ;
    HEADER8 *h8 ;		/* Only when reclen==8	*/
    long *oblock ; 		/* From HEADER		*/
    short *mags ;		/* From HEADER		*/
    long index[145]; 		/* Index of USNO in SDz	*/
} PMMroot ;

static int stopat_eof ;		/* Used by pmm_read	*/
static int swapping ;		/* Set to 1 or 2 (swap)	*/

/*==================================================================
		Simple transformations
 *==================================================================*/

char *pmm_dir(int zone)
/*++++++++++++++++
.PURPOSE  Convert the zone number (0 .. 71) into directory
.RETURNS  Its equivalence
.REMARKS  Internal buffer returned
-----------------*/
{
  static char name[16];
    zone *= 25 ;
    if (zone < 900)
	sprintf(name, "S%02d%02d", 82-(zone/10), zone%10 ? 0 : 30) ;
    else
	sprintf(name, "N%02d%02d", (zone-900)/10, zone%10 ? 30 : 0) ;
    return(name) ;
}

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

static void swap1 (array, 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
-----------------*/
  short *array ;  /* IN: The array to convert   */
  int nshort ;    /* IN: The number of integers */
{
  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 ;
    }
}

static void swap (array, 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
-----------------*/
  int *array ;  /* IN: The array to convert   */
  int nint ;    /* IN: The number of integers */
{
  char *p, *e ;
  int n ;
  
    p = (char *)array ; e = p + 4*nint ;
    if (swapping == 1) swap1((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 void swap_h (h)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of integers if necessary
.RETURNS  0/1/2 (type of swap)
.REMARKS  Useful for big-endian machines
-----------------*/
  HEADER *h;
{
  HEADER8 *h8 = (HEADER8 *)0;
  int n = 8 ;
    swap1(&(h->reclen), 2) ;
    swap (&(h->ntab), 8) ;
    if (h->reclen == 8) {
        h8 = (HEADER8 *) h ;
	swap1(h8->plate, ITEMS(h8->plate)) ;
	swap (h8->epoch, ITEMS(h8->epoch)) ;
	swap1(h8->mags , ITEMS(h8->mags )) ;
	swap(h8->oblock, INDSIZ8) ;
    }
    else {
	swap1(h->plate, ITEMS(h->plate)) ;
	swap (h->epoch, ITEMS(h->epoch)) ;
	swap1(h->mags , ITEMS(h->mags )) ;
	swap(h->oblock, INDSIZ) ;
    }
}

static long ra8(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  long ra ;
    ra = ((pc[0]&0x3)<<16) | (pc[1]<<8) | pc[2] ;
    return(ra) ;
}

static long ra7(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  long ra ;
    ra = ((pc[0]&7) <<8) | pc[1] ;
    return(ra) ;
}

static int binloc(long *start, int items, long value)
/*++++++++++++++++
.PURPOSE  Find index in table such that    start[i] <= value < start[i+1]
.RETURNS  index / -1 (< *start) / items (>= start[items])
-----------------*/
{
  long *low, *high, *med ;
    if (value < *start) return(-1) ;
    low=start, high=start+items-1; 
    if (value >= *high) return(items) ;
    while ((high-low) > 1) {
    	med = low + (high-low)/2 ;
    	if (value < *med) high = med;
    	else low = med ;
    }
    return(low - start) ;
}

static int binlocf(char *start, int bytes, int reclen, long (*f)(), long 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) ;
}

static int decode_mag(int mag, short *mags)
/*++++++++++++++++
.PURPOSE  Decode the magnitude
.RETURNS  Its original value
-----------------*/
{
  int m ;
    if (mag >= 200) m = mags[mag-200] ;
    else m = mag ;
    return(m+50) ;
}

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

    if (PMMroot.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 && (pmm_options&1)) fprintf(stderr, 
        "#...swapping type=%d\n", swapping) ;

    /* Choose the root name */
    if (root)  PMMroot.root = root ;
    if (!PMMroot.root) PMMroot.root = getenv("PMMroot") ;
    if (!PMMroot.root) PMMroot.root = "/PMM" ;
    n = 20 + strlen(PMMroot.root);	/* Size of file name */
    if (n < 128) n = 128 ;
    PMMroot.name = malloc(n) ;
    PMMroot.SDz  = PMMroot.RAz = -1 ;
    PMMroot.mbuf = 7*2048 ;
    PMMroot.abuf = malloc(PMMroot.mbuf) ;
    PMMroot.line = malloc(n + sizeof(header_text)) ;
    PMMroot.h    = (HEADER *)malloc(sizeof(HEADER)) ;
}

static int pmm_fopen(char *name)
/*++++++++++++++++
.PURPOSE  Open a name-specified file
.RETURNS  0 = OK / -1 = error
-----------------*/
{
  char buffer[128] ;
  int file, len, reclen, i ;

    if (pmm_options&1) printf("#....pmm_fopen(%s)\n", name) ;
    if (!PMMroot.root) init((char *)0) ;
    if (PMMroot.fno > 0) {
	if (close(PMMroot.fno)<0) perror(PMMroot.name) ;
	PMMroot.fno = -1 ;
    }
    file = open(name, O_BINARY);
    if (file <= 0) {  perror(name);  return (-1) ; }
    len = read(file, buffer, sizeof(buffer)) ;
    if (len <= 0) {  perror(name);  return (-1); }

    /* Verify magic */
    reclen = 0 ;
    if (strncmp(buffer, header_text, strlen(header_text)) == 0) 
	reclen = 7 ;
    else if (strncmp(buffer, header8text, strlen(header8text)) == 0)
	reclen = 8 ;
    if (reclen == 0) {
    	fprintf(stderr,  "****File %s: bad magic, should start by\n\"%s\"\n", 
    	    name, header_text) ;
    	return (-1);
    }
    for (i=3; (i<len) && (buffer[i] != '\n'); i += 4) ;
    if (i >= len) {
    	fprintf(stderr,  "****File %s: bad magic, should start by\n\"%s\"\n", 
    	    name, header_text) ;
    	return (-1);
    }
    len = i+1 ;		/* Length of first line */
    lseek(file, len, 0) ;

    /* Save now everything in PMMroot */
    if (PMMroot.name != name) strcpy(PMMroot.name, name) ;
    PMMroot.fno  = file ;
    PMMroot.nbuf = PMMroot.ibuf = 0 ;
    PMMroot.obuf = 0 ;
    memcpy(PMMroot.line, buffer, len); 
    PMMroot.line[len-1] = 0 ;		/* Terminate 1st line	*/
    i = read(file, PMMroot.h, sizeof(HEADER)) ;
    if (i != sizeof(HEADER)) perror(name) ;

    /* Swap the Header */
    if (swapping) swap_h(PMMroot.h) ;

    /* Set the record characteristics */
    PMMroot.mags   = PMMroot.h->mags,
    PMMroot.oblock = PMMroot.h->oblock ;
    if (PMMroot.h->reclen == 8) {
    	PMMroot.h8 = (HEADER8 *)PMMroot.h ;
	PMMroot.mags   = PMMroot.h8->mags, 
	PMMroot.oblock = PMMroot.h8->oblock ;
    }
    else PMMroot.h8 = (HEADER8 *)0 ;

    return(0) ;
}

static long aseek(HEADER *h, long value)
/*++++++++++++++++
.PURPOSE  Position the file to specified RA value (in 10mas)
.RETURNS  The UNSO-ID / -1
-----------------*/
{
  long *oblock, *start, o, opos, da, Da, N, n ;
  long (*fra)() ;
  char *p, *e ;
  int k, nbuf ;

    if ((value < h->ra0) || (value >  h->ra1)) return(-1) ;
    da = value - h->ra0 ;

    /* Test first in the offset index table */
    if (h->reclen==8) 
	fra = ra8, Da = 1<<18, k = (da>>18),
	start = ((HEADER8 *)h)->oblock ;
    else
	fra = ra7, Da = 1<<11, k = (da>>11),
	start = h->oblock ;
    da -= k*Da ;
    oblock = start + k ;

    /* When the size is very large, try to choose a position */
    N = (oblock[1] - oblock[0])/h->reclen ;
    n = (da*N/Da) ;
    if (n > 2048) 	/* Assume constant repartition */
    	n -= n/8 ;
    else n = 0 ;
    opos = oblock[0] + n*h->reclen ;

  Get_1_Block:
    if ((opos >= PMMroot.obuf) && (opos < PMMroot.obuf+PMMroot.nbuf)) {
    	PMMroot.ibuf = opos -  PMMroot.obuf ;
    }
    else {	/* Need to read a block ... */
    	if (lseek(PMMroot.fno, opos, 0) < 0) { 
	    perror(PMMroot.name); 
	    return(-1) ; 
	}
    	PMMroot.obuf = opos ;
    	PMMroot.nbuf = read(PMMroot.fno, PMMroot.abuf, PMMroot.mbuf) ;
    	if (PMMroot.nbuf < 0) { perror(PMMroot.line); return(-1) ; }
    	PMMroot.ibuf = 0 ;
    }
    p = PMMroot.abuf + PMMroot.ibuf ;

    if (((*fra)((unsigned char *)p) > da) && (opos > oblock[0])) {
	opos = oblock[0] ;
	PMMroot.nbuf = 0 ;
	goto Get_1_Block ;
    }

    /* Check now that the required RA is effectively in the read block */
    o = PMMroot.obuf + PMMroot.nbuf ;
    e = PMMroot.abuf + PMMroot.nbuf - h->reclen ;	/* Last record */
    while ((o < oblock[1]) && ((*fra)((unsigned char *)e) < da)) {
	opos = PMMroot.obuf = o ;
	PMMroot.nbuf = read(PMMroot.fno, PMMroot.abuf, PMMroot.mbuf) ;
	if (PMMroot.nbuf < 0) { perror(PMMroot.line); return(-1) ; }
	PMMroot.ibuf = 0 ;
	o = opos + PMMroot.nbuf ;
	e = PMMroot.abuf + PMMroot.nbuf - h->reclen ;
    }

    /* Find position in current buffer */
    p = PMMroot.abuf + PMMroot.ibuf ;
    n = PMMroot.nbuf - PMMroot.ibuf ;	/* Number of bytes in current buffer */
    if ((PMMroot.obuf + PMMroot.nbuf) >= oblock[1]) 	/* Buf. spans blocks */
	n = oblock[1] - PMMroot.obuf - PMMroot.ibuf ;
    n = binlocf(p, n,  h->reclen, fra, da) ;
    if (n < 0) n = 0 ;
    /* if ((*fra)(PMMroot.abuf+n) < da) n += h->reclen ; */
    PMMroot.ibuf += n ;
    o = PMMroot.obuf + PMMroot.ibuf ;

    return(h->id0 + (o - *start)/h->reclen) ;
}

static int zseek(int RAz, int SDz)
/*++++++++++++++++
.PURPOSE  Seek to specified zone
.RETURNS  0 / --1
-----------------*/
{
  int file ;
  char *p ;

    if (!PMMroot.root) init((char *)0) ;
    SDz = 3*(SDz/3) ; 

    p = (char *)0 ;
    if ((PMMroot.fno < 0) || (PMMroot.SDz != SDz)) {
	p = PMMroot.name ;
	strcpy(p, PMMroot.root) ; p += strlen(p) ;
	*(p++) = '/' ; 
	strcpy(p, pmm_dir(SDz)) ; p += strlen(p) ;
	*(p++) = '/' ; 

	/* Load index */
	PMMroot.SDz = SDz ; PMMroot.RAz = -1 ;
    	memset(PMMroot.index, 0, sizeof(PMMroot.index)) ;
	strcpy(p, "index.145") ;
	file = open(PMMroot.name, 0) ;
	if (file > 0) {
	    read(file, PMMroot.index, sizeof(PMMroot.index)) ,
	    close(file) ;
	    if (swapping) swap(PMMroot.index, 145) ;
	}
    }

    /* Go to RA zone */
    if (RAz < 0) return(0) ;

    RAz = 3*(RAz/3) ;
    if (PMMroot.RAz != RAz) {
	if (!p) {
	    p = PMMroot.name ;
	    strcpy(p, PMMroot.root) ; p += strlen(p) ;
	    *(p++) = '/' ; 
	    strcpy(p, pmm_dir(SDz)) ; p += strlen(p) ;
	    *(p++) = '/' ; 
	}
	sprintf(p, "%02d%d0.pmm", RAz/6, RAz%6) ;
	if (pmm_fopen(PMMroot.name) < 0) return(-1) ;
	PMMroot.RAz = RAz ;
    }
    return(0) ;
}

/*==================================================================
		Convert the Input Record(s)
 *==================================================================*/
static int ed_rec(PMMrec *pr)
/*++++++++++++++++
.PURPOSE  Convert the current compressed PMM record into its standard structure
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
-----------------*/
{
  long opos, value ;
  static int oo = 0 ;
  unsigned char *pc ;
  int i ;

    /* Find out which of the RA subparts */
    opos = PMMroot.obuf + PMMroot.ibuf ;
    if ((oo < (PMMroot.h->ntab)-1) && (opos >= PMMroot.oblock[oo])) {
    	while (opos >= PMMroot.oblock[oo+1]) oo++ ;
    }
    else oo = PMMroot.h->ntab;
    if (oo >= PMMroot.h->ntab) 	/* Binary Search */
    	oo = binloc(PMMroot.oblock, PMMroot.h->ntab, opos) ;

    /* Convert the compacted record */
    pc = (unsigned char *)(PMMroot.abuf + PMMroot.ibuf) ;
    pr->zone = PMMroot.h->zone ;  
    pr->id   = PMMroot.h->id0 + (opos-PMMroot.oblock[0])/PMMroot.h->reclen ;
    pr->flags[0] = pc[0]&PMM_GSC ? 'A' : ' ' ; 
    pr->flags[1] = pc[0]&PMM_m   ? '*' : ' ' ; 
    if (PMMroot.h8) {	/* 8 bytes/rec */
	i = (pc[0]>>2) & 0x0f ;
    	pr->epoch = PMMroot.h8->epoch[i] ; 
    	pr->plate = PMMroot.h8->plate[i] ; 
	value = (oo<<18) | ra8(pc) ;
	pc++ ;
    }
    else {	/* 7 bytes/rec */
    	i = (pc[0]>>3) & 0x07 ;
    	pr->epoch = PMMroot.h->epoch[i] ; 
    	pr->plate = PMMroot.h->plate[i] ; 
	value = (oo<<11) | ra7(pc) ;
    }
    pr->ra = PMMroot.h->ra0 + value ;
    value  = (pc[2]<<16) | (pc[3]<<8) | pc[4] ;
    pr->sd = PMMroot.h->sd0 + (value>>2) ;

    value  = ((pc[4]&0x3)<<8) | pc[5] ;
    pr->mB = decode_mag(value>>1, PMMroot.mags) ;
    value  = ((pc[5]&0x1)<<8) | pc[6] ;
    pr->mR = decode_mag(value, PMMroot.mags) ;
    pr->ci = pr->mB - pr->mR ;

    pr->rho = -1 ;

    return(0) ;
}

static int set_id(long id)
/*++++++++++++++++
.PURPOSE  Set the current record to id 
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
-----------------*/
{
  long opos, value ;
  unsigned char *pc ;
  int i ;

    if ((id < PMMroot.h->id0) || (id > PMMroot.h->id1)) {
    	fprintf(stderr, 
    	  "++++Bad USNO number (%d.%08ld), should be in range [%ld,%ld]\n", 
    	  PMMroot.h->zone, id, PMMroot.h->id0, PMMroot.h->id1) ;
    	fprintf(stderr, "    (%s)\n", PMMroot.line);
  	return(-1) ;
    }
    opos = PMMroot.oblock[0] + (id - PMMroot.h->id0)*PMMroot.h->reclen ;

    /* Load the part of file */
    if ((opos < PMMroot.obuf) || (opos >= (PMMroot.obuf + PMMroot.nbuf))) {
    	PMMroot.obuf = opos ;
    	if (lseek(PMMroot.fno, opos, 0) < 0) perror(PMMroot.name) ;
    	PMMroot.nbuf = read(PMMroot.fno, PMMroot.abuf, PMMroot.mbuf) ;
    	if (PMMroot.nbuf < 0) { perror(PMMroot.line); return(-1) ; }
	PMMroot.ibuf = 0 ;
    }
    else PMMroot.ibuf = (opos - PMMroot.obuf) ;
    
    return(0) ;
}

int pmm_set(int zone, long id)
/*++++++++++++++++
.PURPOSE  Position in specified zone + id
.RETURNS  0 / -1
-----------------*/
{
  int i, z ;

    if (!PMMroot.root) init((char *)0) ;

    /* Go to specified zone */
    z = zone/25;
    if (zseek(-1, z) < 0) return(-1) ;

    if ((PMMroot.RAz >= 0) && (id >= PMMroot.h->id0) && (id <= PMMroot.h->id1))
	return(set_id(id)) ;
    if (!PMMroot.index[0]) {
	fprintf(stderr, "****Can't pmm_set(%d, %ld): no index %s\n",
	    zone, id, PMMroot.name) ;
	return(-1) ;
    }
    i = binloc(PMMroot.index, 145, id) ;
    if ((i < 0) || (i >= 144)) {
	fprintf(stderr, "****Can't pmm_set(%d, %ld): not in range [0,%ld[\n", 
	zone, id, PMMroot.index[144]) ;
	return(-1) ;
    }
    if (zseek(i, z) < 0) return(-1) ;
    return(set_id(id)) ;
}

int pmm_read(PMMrec *rec)
/*++++++++++++++++
.PURPOSE  Read next record in PMM
.RETURNS  0(End) / -1 (Err) / 1(OK)
.REMARKS  When end if a RA file is found, action depends on stopat_eof
-----------------*/
{
  int i ; 
  long id, opos ;
    if (!PMMroot.abuf) {
	fprintf(stderr, "****Can't pmm_read before setting\n") ;
	return(-1) ;
    }
    if (PMMroot.ibuf >= PMMroot.nbuf) {		/* Load next record */
	PMMroot.obuf += PMMroot.ibuf ;
    	PMMroot.nbuf = read(PMMroot.fno, PMMroot.abuf, PMMroot.mbuf) ;
    	if (PMMroot.nbuf < 0) { perror(PMMroot.line); return(-1) ; }
	PMMroot.ibuf = 0 ;
	if (PMMroot.nbuf == 0) {		/* End of file !    */
	    if (stopat_eof) return(0) ;
	    id = PMMroot.h->id1 + 1 ;
	    if (id >= PMMroot.index[144]) return(0) ;
	    if (pmm_set(PMMroot.h->zone, id) < 0) return(-1) ;
	}
    }
    ed_rec(rec) ;
    PMMroot.ibuf += PMMroot.h->reclen ;
    return(1) ;
}

int pmm_get(int zone, long id, PMMrec *rec)
/*++++++++++++++++
.PURPOSE  Read specified record in PMM
.RETURNS  0 / -1
-----------------*/
{
    if (pmm_set(zone, id) < 0) return(-1) ;
    return(pmm_read(rec)) ;
}

static int edra(char *buf, long ra, int opt /* 0=deg, 1=sexa */)
/*++++++++++++++++
.PURPOSE  Edit a RA (unit=10mas) into HH:MM:SS.SSS
.RETURNS  Number of bytes (12 or 10)
-----------------*/
{
  int i ; 
  long value ;
    if (opt) {		/* Sexagesimal */
    	value = ra ;
    	value = value - (value/3) ;	/* 2/3 for number in ms */
    	for (i=11; 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(12) ;
    }
    else {
    	value = 3*ra ;	/* 10^6/360000 = 3 * (1 - 2/27) */
    	value -= (2*value)/27 ;
    	for (i=9; i>3; i--) buf[i] = '0' + value%10, value /= 10 ;
    	buf[3] = '.' ;
    	buf[2] = '0' + value%10, value /= 10 ;
    	buf[1] = '0' + value%10, value /= 10 ;
    	buf[0] = '0' + value ;
    	return(10) ;
    }
}

static int edsd(char *buf, long sd, int opt /* 0=deg, 1=sexa */)
/*++++++++++++++++
.PURPOSE  Edit a SPD (unit=10mas) into +DD:MM:SS.SS
.RETURNS  Number of bytes (12 or 10)
-----------------*/
{
  int i ; 
  long value ;
  value = sd - 90*3600*100 ;
  if (value > 0) buf[0] = '+' ;
  else buf[0] = '-', value = -value ;
    if (opt) {		/* Sexagesimal */
    	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(12) ;
    }
    else {
    	value = 3*value ; 		/* 10^6/360000 = 3 * (1 - 2/27) */
    	value -= ((2*value)/27) ; 
    	for (i=9; i>3; i--) buf[i] = '0' + value%10, value /= 10 ;
    	buf[3] = '.' ;
    	buf[2] = '0' + value%10, value /= 10 ;
    	buf[1] = '0' + value ;
    	return(10) ;
    }
}

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

static int edplate(char *buf, int plate) 
/*++++++++++++++++
.PURPOSE  Edit a plate
.RETURNS  Number of bytes (4)
-----------------*/
{
  int value ;
  value = plate ;
    if (value >= 0) buf[0] = ' ' ;
    else buf[0] = '-', value = -value ;
    buf[2] =  buf[1] = ' ' ; buf[3] = '0' ;
    if (value == 0) return(4) ;
    buf[3] = '0' + value%10, value /= 10 ;
    buf[2] = '0' + value%10, value /= 10 ;
    buf[1] = '0' + value ;
    return(4) ;
}

static int edepoch(char *buf, long epoch) 
/*++++++++++++++++
.PURPOSE  Edit the epoch of plate (originally in 0.001yr)
.RETURNS  Number of bytes (8)
-----------------*/
{
  long value ; int i ;
    value = epoch ;
    for (i=7; i>4; i--) buf[i] = '0' + value%10, value /= 10 ;
    buf[i--] = '.' ;
    while (value>0) buf[i--] = '0' + value%10, value /= 10 ;
    while (i>=0)    buf[i--] = ' ' ;
    return(7-i) ;
}

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

char *pmm_head(int opt /*0=deg, 1=sexa, 2=ID, 4=Ep, 8=x,y*/)
/*++++++++++++++++
.PURPOSE  Get the Header for one edition
.RETURNS  The Header, starting by #
-----------------*/
{
  static char buf[80] ;
  char *p, *e ; int i, b ;
  long value ;

    p = buf ; b = 0 ;
    if (opt&2) {	/* Edit the USNO-ID */
    	sprintf(p, "#%-12s ", "USNO-A2.0") ;
    	p += strlen(p) ;
    }
    else b = 1, *(p++) = '#' ;

    /* Edit the Position */
    i = edra(p, 0, opt&1) + edsd(p, 0, opt&1) - b ;
    b = (i-14)/2 ;
    while(--b >= 0) *(p++) = ' ', i-- ;
    strcpy(p, "RA  (ICRS) Dec"), i -= 14, p += 14 ;
    while(--i >= 0) *(p++) = ' ' ;

    /* Flags */
    *(p++) = 'A' ;
    *(p++) = '*' ;

    /* Magnitudes */
    i = edmag(p,0) - 4;
    while (--i >= 0) *(p++) = ' ' ;
    strcpy(p, "Bmag"), p += 4 ;
    i = edmag(p,0) - 4;
    while (--i >= 0) *(p++) = ' ' ;
    strcpy(p, "Rmag"), p += 4 ;

    /* Field */
    if (opt&4) strcpy(p, "  Epoch  ") ;
    else       strcpy(p, " Pl.") ;
    p += strlen(p) ;

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

    *p = 0 ;

    return(buf) ;
}

static int ed10mas(char *buf, long value)
/*++++++++++++++++
.PURPOSE  Edit a value in units of 10mas
.RETURNS  Number of bytes written
-----------------*/
{
  static char aval[16] ;
  char s, *p, *b ;  long x ;
    if (value < 0) s = '-', x = -value ;
    else s = 0, x = value  ;
    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) ;
}

char *pmm2a(PMMrec *pr, int opt /*0=deg, 1=sexa, 2=ID, 4=Date, 8=(x,y)*/)
/*++++++++++++++++
.PURPOSE  Edit (in a static buffer) one record.
.RETURNS  The edited record
-----------------*/
{
  static char buf[100] ;
  char *p ; int i ;
  long value ;
    p = buf ;

    if (opt&2) {	/* Edit the USNO-ID */
    	sprintf(p, "%04d-%08ld ", pr->zone, pr->id) ;
    	p += strlen(p) ;
    }
    /* Edit the Position */
    p += edra(p, pr->ra, opt&1) ;
    p += edsd(p, pr->sd, opt&1) ;

    /* Flags */
    *(p++) = pr->flags[0] ;
    *(p++) = pr->flags[1] ;

    /* Magnitudes */
    p += edmag(p, pr->mB) ;
    p += edmag(p, pr->mR) ;
 
    /* Field */
    if (opt&4) *(p++) = ' ' , p += edepoch(p, pr->epoch) ;
    else       p += edplate(p, pr->plate) ;

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

    *p = 0 ;

    return(buf) ;
}

static int prec(PMMrec *rec) { printf("%s\n", pmm2a(rec, 3)); return(1); }

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

long pmm_seek(long ra, long sd)
/*++++++++++++++++
.PURPOSE  Seek to specified position
.RETURNS  The first USNO-ID number
-----------------*/
{
  long id ;

    if (!PMMroot.root) init((char *)0) ;
    if (zseek(ra/XIZE, sd/XIZE) < 0) return(-1) ;

    /* Seek now until the specified RA found */
    id = aseek(PMMroot.h, ra) ;

    return(id) ;
}

int pmm_close()
/*++++++++++++++++
.PURPOSE  Close the currently opened part
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    if (close(PMMroot.fno) < 0) perror(PMMroot.name) ;
    PMMroot.SDz = PMMroot.RAz = -1 ;
    PMMroot.fno = -1 ;
    PMMroot.obuf = 0 ;
    PMMroot.nbuf = PMMroot.ibuf = 0 ;
    PMMroot.line[0] = 0 ;
    memset(PMMroot.h, 0, sizeof(HEADER)) ;
    memset(PMMroot.index, 0, sizeof(PMMroot.index)) ;
    return(0) ;
}

static int s2int(int *a1, int *a2)
/*++++++++++++++++
.PURPOSE  Sort routine (called by qsort), 2 integers (index distance)
.RETURNS  <0 0 >0
-----------------*/
{
  int diff ;
    diff = a1[1] - a2[1] ;
    if (!diff) diff = a1[0] - a2[0] ;
    return(diff) ;
}

long pmm_search(long ra[2], long sd[2], int (*digest_routine)())
/*++++++++++++++++
.PURPOSE  Search USNO stars within range of positions.
.RETURNS  Number of tested records.
.REMARKS  The supplied digest_routine(PMMrec *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 pmm_options&0x10 ==> Keep Original Order  <===
-----------------*/
{
  static int sin1000[] = {	/* Approx. sine in 10**-3 of 3.75+d*7.5deg */
     65, 195, 321, 442, 555, 659, 751, 831, 896, 946, 980, 997,
    997, 980, 946, 896, 831, 751, 659, 555, 442, 321, 195,  65
  } ;
  static PMMrec rec ;
  int zones[24*48*2], nzones, iz;	/* Table index d*48+a, distance   */
  long ld[2], la[2], d, amin, dmin, dmax, a, amax, dd[2], dsd ;
  long id, tested_records ;
  int goon, reclen, od, bytes;
  int ja, jd, ja0, jd0, dja, djd ;
  int saved_eof ;
  unsigned char *p, *e ;

    if (!PMMroot.root) init((char *)0) ;
    if (!digest_routine) digest_routine = prec ;
    saved_eof = stopat_eof ; stopat_eof = 1 ;
    nzones = 0 ;

    if (!ra) {
	if (!sd)	/* Whole Sky */
	    for (jd=0; jd<24; jd++) for (ja=0; ja<48; ja++, nzones += 2) 
	    zones[nzones] = nzones<<1;
	ra = la ;
        ra[0] = 0 ;
	ra[1] = (48*Dstep)-1 ;
    }
    if (!sd) {
        sd = ld ;
	sd[0] = 0;
	sd[1] = 24*Dstep ;
    }
    ld[0] = sd[0] ;
    ld[1] = sd[1] ;
    la[0] = ra[0] % (48*Dstep) ;
    la[1] = ra[1] % (48*Dstep) ;
    tested_records = 0 ;

    if (ld[0] > ld[1]) d=ld[1], ld[1]=ld[0], ld[0]=d ;
    dmin = Dstep*(ld[0]/Dstep);
    dmax = Dstep*(ld[1]/Dstep) ;
    amin = Dstep*(la[0]/Dstep);
    jd0 = (sd[0]+sd[1])/2/Dstep;
    ja0 = (la[0]+la[1])/2/Dstep;
    if (la[0] > la[1]) ja0 = (ja0+24)%48 ;

    /* Loop to define the zones, with their distance */
    if (nzones) 		/* Whole Sky */
        goto Search ;
    for (d=dmin ; d <= dmax ; d += Dstep) {
	if (d >= (Dstep*24)) continue ;	/* >= 90deg */
	for (a=amin, goon=1; goon; a = (amax+1)%(48*Dstep)) {
	    amax = a+Dstep-1 ;
	    if ((la[1] >= a) && (la[1] <= amax)) goon = 0;
	    jd = d/Dstep ; ja = a/Dstep ;
	    djd = (jd - jd0)*1000; 
	    dja = ja - ja0; if (dja >= 24) dja -= 24 ;
	    dja = dja*sin1000[jd]; 
	    zones[nzones++] = 48*jd + ja ;
	    zones[nzones++] = dja*dja + djd*djd ;
	}	/* EOL on  RAzone */
    }		/* EOL on  SDzone */

    /* Sort according to distances */
    if (pmm_options&0x10) ;
    else qsort(zones, nzones/2, 2*sizeof(int), (SORT_FCT)s2int);

    /* Loop on search in the sorted cells */
  Search:
    for (iz=0, goon=1; goon && iz<nzones; iz += 2) {
	jd = zones[iz]/48 ;
	ja = zones[iz]%48 ;
	d = jd * Dstep ;
	a = ja * Dstep; amax = a+Dstep-1 ;
	if ((la[0] > a) && (la[0] <= amax)) a = la[0]; 
	id = pmm_seek(a, d) ;
	if (id < 0) continue ;
	if ((la[1] >= PMMroot.h->ra0) && (la[1] <= PMMroot.h->ra1))
	    amax = la[1] ;
	else amax = PMMroot.h->ra1 ;
	/* Define limits relative to current file */
	reclen = PMMroot.h->reclen ;
	od = reclen - 5 ;
	dd[0] = sd[0] - PMMroot.h->sd0 ;
	dd[1] = sd[1] - PMMroot.h->sd0 ;
	if (dd[0] < 0) dd[0] = 0 ;
	if (dd[1] > PMMroot.h->sd1) dd[1] = PMMroot.h->sd1 ;

	/* Test the loaded block for declination */
	bytes = PMMroot.nbuf ;
	while(bytes > 0) {
	    for (p = (unsigned char *)PMMroot.abuf + PMMroot.ibuf + od, 
		 e = (unsigned char *)PMMroot.abuf + PMMroot.nbuf ; p < e ;
		 p += reclen, tested_records++) {
		dsd = (p[0]<<16)|(p[1]<<8)|p[2]; dsd >>= 2 ;
		if (dsd<dd[0]) continue ;
		if (dsd>dd[1]) continue ;
		break ;
	    }
	    PMMroot.ibuf = (p-(unsigned char *)PMMroot.abuf)-od;
	    bytes = pmm_read(&rec) ;
	    if (bytes <= 0) continue ;
	    tested_records++;
	    if (rec.ra > amax) break ;
	    if (rec.sd < ld[0]) continue ; 
	    if (rec.sd > ld[1]) continue ;
	    if ((*digest_routine)(&rec) < 0) { 	/* Ask to stop ! */
	        goon = 0 ; 
	        break ;
	    }
	}	/* EOL on a block */
    }

    stopat_eof = saved_eof ;
    return(tested_records) ;
}

/*==================================================================
		Main Program
 *==================================================================*/
#ifdef TEST
static char usage[] = "\
Usage: pmmsub [.|-q]\n\
";

static char help[] = "\
		      -q: search by position\n\
   compressed_input_zone: one of the files created by pmmake\n\
";
main (argc, argv) int argc; char **argv;
{
  char buf[BUFSIZ], *p ;
  long id = 0, ra, sd ;
  double q[2], radius ;
  int  getpos = 0 ;
  int  z;
  PMMrec rec ;

    pmm_options = 0 ;
    if (argc != 2) {
	fprintf(stderr, "%s%s", usage, help) ;
	exit(1) ;
    }
    ++argv;
    if (strcmp(*argv, "-q") == 0) getpos = 2 ;
    else ; 

    while (1) {
    	if (isatty(0)) { 
	    if (getpos) 
		printf("----Give a position (degrees): ") ;
    	    else printf("----Give a zone USNO number: ") ;
    	    if (! gets(buf)) break ;
        } else buf[0] = 0 ;
	if (buf[0]) {
	    for (p=buf; isspace(*p); p++) ;
	    q[0] = atof(p) ;
	    while (isgraph(*p)) p++ ;
	    q[1] = atof(p) ;
	    if (getpos) {
	    	ra = q[0]*3.6e5 ; sd = 3.6e5*(90. + q[1]) ;
	    	id = pmm_seek(ra, sd) ;
	    } 
	    else {
		z = q[0] ; id = q[1];
		pmm_set(z, id) ;
	    }
	}
    	if (pmm_read(&rec)>0) printf("%s\n", pmm2a(&rec, 3)) ;
	else printf("====EOF====\n");
    }
    pmm_close() ;
    exit(0);
}
#endif