/*++++++++++++++
.IDENTIFICATION ucacsub.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    UCAC1 Catalogue, compressed version
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   21-Apr-2001
.VERSION  1.1   13-Jan-2004: Problem of 0h ???
.COMMENTS       Edition of UCAC data.
---------------*/

#include <ucac_def.h>	/* Structure definitions */
#include <ucac.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 EDITED_LEN 	150	/* Maximal length of edited record */
#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 UCAC queries.......	
   Flags used:
   0x01 = Verbose option
   0x10 = Read files by increasing distance from the center.
   
*/
int ucac_options ;	/* Shared with main prog */

static struct {
    char *root;			/* Root name of UCAC	*/
    char *name ;		/* Buffer for full name	*/
    int fno;			/* Opened file number 	*/
    short zone ;		/* Opened zone		*/
    short chunkno;		/* Current chunk 	*/
    int mbuf, nbuf ;		/* Max, used, offset	*/
    unsigned char *abuf;	/* Allocated data	*/
    char line[80] ;		/* First line of file	*/
    long oblock[80] ; 		/* List of Chuncks	*/
    CHUNK dat ;			/* Full details chunk	*/
} UCat ;

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

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

static void swap1 (short *array, int nshort)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of shorts if necessary
.RETURNS  ---
.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 ;
    }
}

static void swap (int *array, int nint)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of integers if necessary
.RETURNS  ---
.REMARKS  Useful for big-endian machines
-----------------*/
{
  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 long ra4(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  long ra ;
    ra = (pc[0]<<24) | (pc[1]<<16) | (pc[2]<<8) | pc[3] ;
    return(ra) ;
}

static long ra3(unsigned char *pc)
/*++++++++++++++++
.PURPOSE  Compute Right Ascension from buffer
.RETURNS  Its value
-----------------*/
{
  long ra ;
    ra = (pc[0]<<16) | (pc[1]<<8) | pc[2] ;
    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 (offseted by dat.omag)
-----------------*/
{
  int m ;
    if (mag >= 950) m = mags[mag-950] ;
    else m = mag ;
    return(m) ;
}

static int decode_pm(int pm, short *apm2, int *apm4)
/*++++++++++++++++
.PURPOSE  Decode the proper motion (13 bits)
.RETURNS  Its original value 
-----------------*/
{
  int m ;
    if (pm&(0x1000)) {	/* Out of bounds */
        m = pm&0xfff ;
	if (m&0x800) 	/* 4-byte value	 */
	     m = apm4[m&0x7ff] ;
	else m = apm2[m&0x7ff] ;
    }
    else m = pm ;
    m -= 2048 ;
    return(m) ;
}

static int decode_no(int no, char *overno)
/*++++++++++++++++
.PURPOSE  Decode the number of observations
.RETURNS  Its original value 
-----------------*/
{
  int m ;
    if (no&0xf0) 	/* Outside range */
         m = overno[no&0xf]; 
    else m = no ;
    m += 2 ;	/* At least 2 obs. / catalogues */
    return(m) ;
}

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

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

    /* Choose the root name */
    if (root)  UCat.root = root ;
    if (!UCat.root) UCat.root = getenv("UCACroot") ;
    if (!UCat.root) UCat.root = "/UCAC" ;
    n = 20 + strlen(UCat.root);	/* Size of file name */
    if (n < 128) n = 128 ;
    UCat.name = malloc(n) ;
    UCat.zone = -1 ;
    UCat.chunkno = -1 ;
    memset(&UCat.dat, 0, sizeof(CHUNK)) ;
}

int ucac_fopen(char *filename)
/*++++++++++++++++
.PURPOSE  Open a file reprsenting one of the UCAC1 files
.RETURNS  0 = OK / -1 = error
-----------------*/
{
  int file, len, reclen, i ;
  char *p ;

    if (ucac_options&1) printf("#...ucac_fopen(%s)\n", filename) ;
    if (!UCat.root) init((char *)0) ;
    if (UCat.fno > 0) {
	if (close(UCat.fno)<0) fprintf(stderr, "****"), 
	    perror(UCat.name) ;
	UCat.fno = -1 ;
	UCat.zone = -1 ;
	UCat.chunkno = -1 ;
    }
    file = open(filename, O_BINARY);
    if (file <= 0) {  
        fprintf(stderr, "****"), perror(filename);  
	return (-1) ; 
    }
    len = read(file, UCat.line, 
          sizeof(UCat.line)+sizeof(UCat.oblock)) ;
    if (len <= 0) {  
        fprintf(stderr, "****"), perror(filename);  
	return (-1); 
    }

    /* Verify magic */
    reclen = 0 ;
    if (len != sizeof(UCat.line)+sizeof(UCat.oblock) ||
        strncmp(UCat.line, header_text, 8) || 
        strncmp(UCat.line+9, header_text+9, strlen(header_text+9))) {
        fprintf(stderr,  
	    "****File %s: got %d/%d bytes, magic should start by\n\"%s\"\n",
	    filename, len, sizeof(UCat.line)+sizeof(UCat.oblock), header_text);
	return(-1) ;
    }

    /* Terminate the Line	*/
    if (p = strchr(UCat.line, '\r')) *p = 0 ;

    /* Fill with file parameters */
    UCat.fno = file ;
    if (p = strchr(UCat.line, 'Z'))
         UCat.zone = atoi(p+1) ;
    else UCat.zone = -1 ;
    UCat.chunkno = -1 ;

    /* Fill with known values */
    UCat.dat.zone   = UCat.zone ;
    UCat.dat.reclen = atoi(UCat.line+7) ;
    UCat.dat.ep0 = 1998000 ;
    UCat.dat.omag = 750 ;
    UCat.nbuf = UCat.dat.len = 0 ;

    /* Swap the Offsets */
    if (swapping) swap((int *)UCat.oblock, ITEMS(UCat.oblock)) ;

    return(0) ;
}

int ucac_zopen(int zone)
/*++++++++++++++++
.PURPOSE  Open a speficied zone
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    if (ucac_options&1) printf("#...ucac_zopen(z%03d)\n", zone) ;
    if (!UCat.root) init((char *)0) ;
    if (UCat.zone == zone) return(0) ;
    sprintf(UCat.name, "%s/z%03d.bin", UCat.root, zone) ;
    return(ucac_fopen(UCat.name));
}

int ucac_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 ~25000 stars, i.e. < 500Kbytes
-----------------*/
{
  long o, len ;
  int  *p; short *s; 

    /* The file must be started and loaded */
    if (ucac_options&1) printf("#...read_chunk(%d) in zone: %d\n", 
        chunkno, UCat.zone) ;
    if (UCat.zone <= 0) return(-1) ;

    /* The table of Offsets gives the location in the File */
    o = UCat.oblock[chunkno] ;
    len = UCat.oblock[chunkno+1] - o ;
    if (len < 0) return(-1) ;
    if (len == 0) {	
        if (UCat.dat.reclen == 18) { 	/* Single Chunk */
	    o = UCat.oblock[chunkno = 0]; 
	    len = UCat.oblock[chunkno+1] - o ;
	}
	else return(-1) ;
    }

    /* Verify chunck not already loaded -- not necessary to reload ! */
    if (UCat.chunkno == chunkno)
        return(0) ;

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

    /* Need to read the whole chunck -- enough memory ? */
    if (UCat.mbuf < len) {
        if (UCat.abuf) free(UCat.abuf) ;
	UCat.mbuf = 1+(len|0xffff) ;	/* Multiple of 64Kbytes */
	UCat.abuf = malloc(UCat.mbuf) ;
    }
    UCat.nbuf = read(UCat.fno, UCat.abuf, len) ;
    if (UCat.nbuf != len) {
	fprintf(stderr, "****Got %d/%d bytes, ", UCat.nbuf, len) ;
        if (UCat.nbuf < 0) perror(UCat.name); 
	else fprintf(stderr, "file: %s\n", UCat.name) ;
	return(-1) ; 
    }

    /* Get the integer*1 -- no swapping problem there ! */
    UCat.dat.ra0  = UCat.abuf[10]<<24 ;
    UCat.dat.sd0  = (UCat.dat.zone-1)*Dstep; 	/* 0.5deg */
    UCat.dat.npm4 = UCat.abuf[11];
    UCat.dat.npm2 = UCat.abuf[12];
    UCat.dat.nmag = UCat.abuf[13];
    memcpy(UCat.dat.overno, &UCat.abuf[14], 16) ;

    /* Interpret the dat bytes */
    if (swapping) {
	p = (int *)UCat.abuf ;
        swap(p, 2) ;		/* length of chunck + IDstart */
	s = (short *)(p+2) ; swap1(s, 1) ;	/* Head length*/
	p = (int *)(UCat.abuf+32) ;
	swap(p, UCat.dat.npm4) ; p += UCat.dat.npm4 ;
	s = (short *)p ;
	swap1(s, UCat.dat.npm2); s += UCat.dat.npm2 ;
	swap1(s, UCat.dat.nmag); 
    }

    /* Terminate the Header Assignations */
    p = (int *)UCat.abuf ;
    UCat.dat.len = *(p++) ;
    UCat.dat.id0 = *(p++) ;
    s = (short *)p ;
    UCat.dat.headlen = *s ;
    UCat.dat.id1 = 
        (UCat.dat.len - UCat.dat.headlen)/UCat.dat.reclen
	+ UCat.dat.id0 - 1 ;
    UCat.dat.ra1 = UCat.dat.ra0 | MASK(24) ;
    if (UCat.dat.reclen == 18) 	/* Single Chunk */
        UCat.dat.ra1 = 360*3600*1000 -1;
    UCat.dat.sd1 = UCat.dat.sd0 + Dstep-1 ;

    /* Out-of-range arrays */
    p = (int *)(UCat.abuf+32) ;
    UCat.dat.amu4 = p ; p += UCat.dat.npm4 ;
    s = (short *)p ;
    UCat.dat.amu2 = s ; s += UCat.dat.npm2 ;
    UCat.dat.amag = s ;

    /* Finalisation: len refers to length of data only */
    UCat.dat.pos   = 0 ;
    UCat.dat.len   = (1+UCat.dat.id1-UCat.dat.id0) * UCat.dat.reclen;
                     /* UCat.nbuf - UCat.dat.headlen + extra_bytes !! */
    UCat.dat.chunk = UCat.abuf + UCat.dat.headlen ;
    UCat.chunkno = chunkno ;

    return(0) ;
}

static long aseek(long value)
/*++++++++++++++++
.PURPOSE  Position the file to specified RA value (in mas)
.RETURNS  The UCAC-ID / -1
-----------------*/
{
  long (*fra)(), a ;
  int da, n ;

    if (read_chunk(value>>24)) return(-1) ;

    if ((value < UCat.dat.ra0) || (value >  UCat.dat.ra1)) 
        return(-1) ;
    da = value - UCat.dat.ra0 ;

    /* Find position in current buffer : binlocf finds
       RA(n) <= value < RA(n+1)) 
    */
    fra = (UCat.dat.reclen == 18 ? ra4 : ra3) ;
    n = binlocf((char *)UCat.dat.chunk, 
    		UCat.dat.len,
                UCat.dat.reclen, fra, da) ;
    if (n < 0) n = 0 ;
    /* We have to have   RA(UCAC-ID) >= value */
    a = (*fra)(UCat.dat.chunk + n) ;
    if (a < da) n += UCat.dat.reclen ;
    if (n >= UCat.dat.len)
	return(-1) ;

    UCat.dat.pos = n ;

    return(UCat.dat.id0 + n/UCat.dat.reclen) ;
}

/*==================================================================
		Convert the Input Record(s)
 *==================================================================*/
static int ed_rec(UCACrec *pr)
/*++++++++++++++++
.PURPOSE  Convert the current compressed UCAC record into its standard structure
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
-----------------*/
{
  unsigned char *pc ;

    /* Convert the compacted record */
    pc = (unsigned char *)(UCat.dat.chunk + UCat.dat.pos) ;
    pr->zone  = UCat.dat.zone ;  
    pr->id    = UCat.dat.id0 + UCat.dat.pos/UCat.dat.reclen ;
    pr->ra    = (pc[0]<<16) | (pc[1]<<8) | pc[2] ;
    if (UCat.dat.reclen == 18) 
	  pr->ra = (pr->ra<<8) | pc[3], pc++ ;
    pr->ra   += UCat.dat.ra0 ;
    pr->sd    = UCat.dat.sd0 + ((pc[3]<<13) | (pc[4]<<5) | (pc[5]>>3));
    pr->epoch = UCat.dat.ep0 + (((pc[5]&MASK(3))<<8) | pc[6]);
    pr->sra   = pc[7] ;
    pr->ssd   = pc[8] ;
    pr->mag   = UCat.dat.omag
              + decode_mag((pc[9]<<2)|(pc[10]>>6), UCat.dat.amag) ;
    pr->no    = decode_no((pc[10]>>1)&MASK(5), UCat.dat.overno) ;
    pr->nc    = 2 + (((pc[10]&1)<<4) | (pc[11]>>4)) ;
    pr->pmra  = decode_pm(((pc[11]&MASK(4))<<9) | (pc[12]<<1) | (pc[13]>>7), 
    		UCat.dat.amu2, UCat.dat.amu4) ;
    pr->spmra = ((pc[13]&MASK(7))<<2) | (pc[14]>>6) ;
    pr->pmsd  = decode_pm(((pc[14]&MASK(6))<<7) | (pc[15]>>1),
    		UCat.dat.amu2, UCat.dat.amu4) ;
    pr->spmsd = ((pc[15]&MASK(1))<<8) | pc[16] ;

    pr->pmtot = -1 ;
    pr->rho   = -1 ;

    return(0) ;
}

int ucac_set(long id)
/*++++++++++++++++
.PURPOSE  Set the current UCAC1number to id 
.RETURNS  -1=Error / 0=OK / 1 = Error (too high)
-----------------*/
{
  long did, o ;
  int i, zone ;

    /* To which zone does the UCAC number correspond to ? */
    i = binloc(UCAC1, ITEMS(UCAC1), id) ;
    if ((i<0) || (i>=ITEMS(UCAC1))) {
        fprintf(stderr, "++++Bad UCAC1 number %08l, should be <%ld\n",
	  id, UCAC1[ITEMS(UCAC1)-1]) ;
	return(-1) ;
    }
    zone = i+1 ;
    if (UCat.zone != zone) {
        if (ucac_zopen(zone)) return(-1) ;
    }
    did = id - UCAC1[i] ;	/* The ID starting from zero in zone */

    /* The currently loaded chunk could be the good one ! */
    if ((UCat.chunkno >= 0) && (id >= UCat.dat.id0) && (id <= UCat.dat.id1)) {
        UCat.dat.pos = (id - UCat.dat.id0)*UCat.dat.reclen ;
	return(0) ;
    }

    /* Move now in the zone to load the appropriate chunck */
    if (UCat.dat.reclen == 18) {
        if (read_chunk(0)) return(-1) ;
    }
    else {	/* A chunck header is >= 32 bytes, i.e. about 2-3 records */
        o = UCat.oblock[0] + did*UCat.dat.reclen ;
	i = binloc(UCat.oblock, ITEMS(UCat.oblock), o) ;
	/* Close to the liit: go to next chucnk */
        if ((UCat.oblock[i+1]-o) < (2*i*UCat.dat.reclen)) i++ ;
        if (read_chunk(i)) return(-1) ;
	if (id > UCat.dat.id1) {
	    if (read_chunk(++i)) return(-1) ;
	}
    }

    if ((id < UCat.dat.id0) || (id > UCat.dat.id1)) {
    	fprintf(stderr, 
    	  "****Bad UCAC1 number %08ld, should be in range [%ld,%ld]\n", 
    	  id, UCat.dat.id0, UCat.dat.id1) ;
    	fprintf(stderr, "    (%s)\n", UCat.line);
  	return(-1) ;
    }

    UCat.dat.pos = (id - UCat.dat.id0)*UCat.dat.reclen ;
    return(0) ;
}

int ucac_read(UCACrec *rec)
/*++++++++++++++++
.PURPOSE  Read next record in UCAC
.RETURNS  0(End) / -1 (Err) / 1(OK)
.REMARKS  When end of chunk is found, action depends on stopat_eof
-----------------*/
{
  int i ; 

    /* First Record in File */
    if (UCat.zone <= 0) {	/* Beginning = Zone 1 */
	if (ucac_zopen(1)) return(-1) ;	
	if (read_chunk(0)) return(-1) ;
	UCat.dat.pos = 0 ;
    }

    /* Buffer exhausted */
    if (UCat.dat.pos >= UCat.dat.len) {	/* Load next record */
	if (stopat_eof>1) return(0) ;	/* Stop at EndOfChunk 	*/
	i = UCat.chunkno+1 ;
	if (i<0) i = 0 ;
	while ((i<78) && (UCat.oblock[i] == UCat.oblock[i+1])) i++ ;
	if (i==78) {
	    /* Zone exhausted! */
	    if (stopat_eof) return(0) ;
	    i = UCat.zone+1 ;
	    if (i >= ITEMS(UCAC1)-1) 	/* Last zone was read */
	        return(0) ;
	    if (ucac_zopen(i)) return(-1) ;
	    i = 0 ;			/* Chunk # to read .. */
	}
	if (read_chunk(i)) return(-1) ;
    }

    ed_rec(rec) ;
    UCat.dat.pos += UCat.dat.reclen ;
    return(1) ;
}

int ucac_get(long id, UCACrec *rec)
/*++++++++++++++++
.PURPOSE  Read specified record in UCAC
.RETURNS  0 / -1
-----------------*/
{
    if (ucac_set(id) < 0) return(-1) ;
    return(ucac_read(rec)) ;
}

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

static int edsd(char *buf, long sd, int opt /* 0=deg, 1=sexa */)
/*++++++++++++++++
.PURPOSE  Edit a SPD (unit=mas) into +DD:MM:SS.SSS
.RETURNS  Number of bytes (13 or 11)
-----------------*/
{
  int i ; 
  long 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*3600*1000 ;
    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)/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
.RETURNS  Number of bytes (5)
-----------------*/
{
  int value ;
  value = mag ;
    buf[0] = buf[1] = buf[2] = buf[4] = ' ' ;
    buf[2] = '.' ;
    if (value < 0) {
	value = -value ;
	buf[0] = '-' ;
    }
    buf[4] = '0' + value%10, value /= 10 ;
    buf[3] = '0' + value%10, value /= 10 ;
    buf[1] = '0' + value%10, value /= 10 ;
    if (value) buf[0] = '0' + value ;
    return(5) ;
}

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 *ucac_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[EDITED_LEN] ;
  char *p, *e ; int i, b ;
  long value ;

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

    /* Edit the Position */
    i = edra(p, 0, opt&17) + edsd(p, 0, opt&17) - b ;
    b = (i-14)/2 ;
    while(--b >= 0) *(p++) = ' ', i-- ;
    if (opt&16) strcpy(p, "RAmas(ICRS)SPD");
    else        strcpy(p, "RA  (ICRS) Dec");
    i -= 14, p += 14 ;
    while(--i >= 0) *(p++) = ' ' ;
    strcpy(p, " +/- +/- ") ; p += strlen(p) ;

    /* Magnitudes */
    i = edmag(p,0) - 5;
    while (--i >= 0) *(p++) = ' ' ;
    strcpy(p, "UCmag"), p += 5 ;

    /* Field */
    strcpy(p, " No  Epoch  ") ;
    p += strlen(p) ;

    /* Proper Motions */
    strcpy(p, "  pmRA(mas/yr)pmDE") ;
    p += strlen(p) ;
    strcpy(p, "  +/-  +/- Nc"); 
    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 edmas(char *buf, long value)
/*++++++++++++++++
.PURPOSE  Edit a value in units of mas to arcsec (with 2 decimals only)
.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  ;
    x = (x+5)/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) ;
}

char *ucac66(UCACrec *pr)
/*++++++++++++++++
.PURPOSE  Edit (in a static buffer) one record in original format
.RETURNS  The edited record
-----------------*/
{
  static char buf[100] ;
    sprintf(buf, "%10d%10d%5d%4d%4d%3d%5d%7d%7d%4d%4d%3d",
      pr->ra, pr->sd, pr->mag, pr->sra, pr->ssd, pr->no, pr->epoch-1997000, 
      pr->pmra, pr->pmsd, pr->spmra, pr->spmsd, pr->nc) ;
    return(buf);
}

char *ucac2a(UCACrec *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[EDITED_LEN] ;
  char *p ; int i ;
  long value ;

    p = buf ;

    if (opt&2) {	/* Edit the UCAC1 */
    	sprintf(p, "%08d ", pr->id) ;
    	p += strlen(p) ;
    }

    /* Edit the Position */
    p += edra(p, pr->ra, opt&17) ;
    p += edsd(p, pr->sd, opt&17) ;

    /* Edit the Sigmas on Position */
    sprintf(p, " %3d %3d", pr->sra, pr->ssd) ;
    p += strlen(p) ;

    /* Magnitudes */
    *(p++) = ' ' ;
    p += edmag(p, pr->mag) ;
 
    /* Number of Observations */
    sprintf(p, " %2d ", pr->no) ;
    p += strlen(p) ;

    /* Epoch */
    p += edepoch(p, pr->epoch) ;

    /* Proper Motions in mas/yr */
    sprintf(p, " %+7d", pr->pmra) ;
    p += strlen(p)-1 ;
    p[1] = p[0]; p[0] = '.'; p += 2 ;
    sprintf(p, " %+7d", pr->pmsd) ;
    p += strlen(p)-1 ;
    p[1] = p[0]; p[0] = '.'; p += 2 ;

    sprintf(p, " %3d", pr->spmra) ;
    p += strlen(p)-1 ;
    p[1] = p[0]; p[0] = '.'; p += 2 ;
    sprintf(p, " %3d", pr->spmsd) ;
    p += strlen(p)-1 ;
    p[1] = p[0]; p[0] = '.'; p += 2 ;

    /* Number of catalogues */
    sprintf(p, " %2d", pr->nc) ;
    p += strlen(p) ;

    /* Distance in arcsec */
    if (pr->rho >= 0)  {
	*(p++) = ' ' ; *(p++) = ';' ; *(p++) = ' ' ; 
	if (opt&8) {	/* 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(UCACrec *rec) { printf("%s\n", ucac2a(rec, 3)); return(1); }

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

long ucac_seek(long ra, long sd)
/*++++++++++++++++
.PURPOSE  Seek to specified position
.RETURNS  The first UCAC1 number / 0
-----------------*/
{
  long id ; int i ;

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

    i = sd/Dstep ;
    if ((i<0) || (i >= ITEMS(UCAC1)-1))
        return(0) ;			/* Outside Boundaries */
    if (ucac_zopen(i+1)) return(-1) ;

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

    return(id) ;
}

int ucac_close()
/*++++++++++++++++
.PURPOSE  Close the currently opened part
.RETURNS  0 = OK / -1 = error
-----------------*/
{
    if (close(UCat.fno) < 0) perror(UCat.name) ;
    UCat.zone = UCat.chunkno = -1 ;
    UCat.fno = -1 ;
    UCat.nbuf = UCat.dat.pos = 0 ;
    UCat.line[0] = 0 ;
    memset(&UCat.dat, 0, sizeof(CHUNK)) ;
    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 = a2[1] - a1[1] ;
    if (!diff) diff = a2[0] - a1[0] ;
    return(diff) ;
}

long ucac_search(long ra[2], long sd[2], int (*digest_routine)())
/*++++++++++++++++
.PURPOSE  Search UCAC1 stars within range of positions.
.RETURNS  Number of tested records.
.REMARKS  The supplied digest_routine(UCACrec *record) does whatever,
	return -1 if stop asked. Default is just to print out.
-----------------*/
{
  static UCACrec rec ;
  long lac[80*2+1], tmp[80*2+1], high, *pamax ;
  long la[2], ld[2], d, dmin, dmax, d0, a ;
  long id, tested_records ;
  int goon, reclen, bytes;
  int ja, jd, jd0, dja, djd, djd1 ;
  int saved_eof ;
  unsigned char *p, *e ;

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

    if (!ra) {
	if (!sd) {	/* Whole Sky */
	    stopat_eof = 0 ;
	    ucac_set(1) ;
	    while(ucac_read(&rec)>0) {
		tested_records++;
		if ((digest_routine)(&rec) < 0) break ;
	    }
	    stopat_eof = saved_eof ;
	    return(tested_records) ;
	}
	stopat_eof = 1 ;
	ra = la ;
        ra[0] = 0 ;
	ra[1] = 360*3600*1000 -1;
    }
    if (!sd) {
        sd = ld ;
	sd[0] = 0;
	sd[1] = 1800*1000*ITEMS(UCAC1) ;
    }
    ld[0] = sd[0] ;
    ld[1] = sd[1] ;
    la[0] = ra[0] % (360*3600*1000);
    la[1] = ra[1] % (360*3600*1000);

    /* Distribute the RA range in the chunks -- 
       high = highest RA value in current chunk
    */
    lac[0] = la[0] ;
    high = (lac[0]&0xff000000) | 0xffffff ;
    ja = 1 ;
    if (la[0] > la[1]) {	/* Around 0h   */
	/* lac[ja++] = high;  -- Modified V1.2 */
	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 zones     */
	
    /* Reshuffle the RA chunks if query around center  */
    if (ucac_options&0x10) {
	ja &= ~1 ;		/* Middle lower bound  */
	dja = 0 ;
	pamax = tmp ;
	for (goon=2; goon; ja += dja) {
	    dja = -dja ;
	    if (dja == 0) dja = lac[ja] & 0x800000 ? 2 : -2 ;
	    else if (dja < 0) dja -= 2 ;
	    else dja += 2 ;
	    if (ja < 0) { goon--; continue; }
	    if (lac[ja] < 0) { goon--; continue; }
	    goon = 2 ;
	    *(pamax++) = lac[ja] ;
	    *(pamax++) = lac[ja+1];
	}
	*pamax = -1 ;
	memcpy(lac, tmp, (pamax - tmp)*sizeof(long)) ;
    }

    /* Compute the middle (in intervals with the step) */
    d0 = (ld[0]+ld[1])/2 ;
    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) ;
    jd0 =  d0/Dstep;

    /* When looking around a center (0x10), search first in strip
       containing the center, then move North/South.
       djd = -1 when the closest strip is below.
    */
    if (ucac_options&0x10) {
	djd = 0 ;
        if ((d0 - jd0*Dstep) <= (Dstep/2)) djd1 = -Dstep ;
	else djd1 = Dstep ;
	d = (jd0*Dstep) ;	/* Starting zone */
    }
    else {
	djd = Dstep ;
	d = dmin ;		/* Starting zone */
    }

    /* Loop on search by Declination Strips 
       goon is 2, which meqns that after 2 consecutive mismatching
       strips, the loop will stop.
    */
    for (goon=2; goon; d += djd) {
	if (ucac_options&0x10) {	/* djd are 0 +1 -2 +3 ... */
	    if (djd == 0) djd = djd1 ;
	    else  {
	        djd = -djd;
		if (djd < 0) djd -= Dstep ;
		else djd += Dstep ;
	    }
	}
	if (d < dmin) { goon--; continue; }
	if (d > dmax) { goon--; continue; }
	goon = 2 ;

	/* Searching around a Position: optimize the chunks to 
	   minimize the comparisons.
	*/
	for (pamax = lac; *pamax >= 0; pamax++) {
	    a = *(pamax++); 
	    if (ucac_options&1) fprintf(stderr, 
		"....Search RA [%d,%d]\n", a, *pamax);
	    id = ucac_seek(a, d) ;
	    if (id <= 0) continue ;
	    while ((ucac_read(&rec)>0) && goon) {
	        tested_records++;
	        if (rec.ra > *pamax) break ;
	        if (rec.sd < ld[0]) continue ; 
	        if (rec.sd > ld[1]) continue ;
	        if ((*digest_routine)(&rec) < 0)  	/* Ask to stop ! */
	            goon = 0 ; 
	    }
	}	/* EOL on a block */
    }

    stopat_eof = saved_eof ;
    return(tested_records) ;
}

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

static char help[] = "\
		      -q: search by position\n\
   compressed_input_zone: one of the files created by ucacake\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;
  UCACrec rec ;

    ucac_options = 1 ;
    if(!getenv("UCACroot")) putenv("UCACroot=../bin");
    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 UCAC1 (zxxx) or UCAC1 number: ") ;
    	    if (! gets(buf)) break ;
        } else buf[0] = 0 ;
	if (buf[0]) {
	    for (p=buf; isspace(*p); p++) ;
	    z = 0 ; id = 0;
	    if (*p == 'z') ucac_zopen(z = atoi(p+1)) ;
	    else {
	        q[0] = atof(p) ;
	        while (isgraph(*p)) p++ ;
		while (isspace(*p)) p++ ;
	        q[1] = atof(p) ;
	    }
	    if (getpos) {
	    	ra = q[0]*3.6e5 ; sd = 3.6e5*(90. + q[1]) ;
	    	id = ucac_seek(ra, sd) ;
	    } 
	    else {
		id = q[0];
		if (z == 0) {
		    if (ucac_set(id)) continue ;
		}
	    }
	}
    	/*if (ucac_read(&rec)>0) printf("%s\n", ucac2a(&rec, 3)) ;*/
    	if (ucac_read(&rec)>0) printf("%s\n", ucac66(&rec)) ;
	else printf("====EOF====\n");
    }
    ucac_close() ;
    exit(0);
}
#endif