/*++++++++++++++
.IDENTIFICATION nomad_make.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    NOMAD Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   09-Nov-2002: USNO-B1
.VERSION  1.1   09-Nov-2005: NOMAD
.COMMENTS       Create the vompressed binary files

	The original structure (ordered by 0.1 degrees in Declination)
	is kept.
	Within each file, records are grouped in chuncks

	See the description in file <nomad1def.h>
---------------*/

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

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

#ifndef O_BINARY	/* Required for MS-DOS	*/
#define O_BINARY 0
#endif
#include <stdio.h>

#define ITEMS(a)	(sizeof(a)/sizeof(a[0]))
typedef int (*SORT_FCT)();

static FILE *infile;		/* Input file ------------- */
static FILE *u1file;		/* OutputFile ------------- */
static char *orifile;		/* Input file name	    */
static int  SDzone ;		/* Current Declination zone */
static int  reclen = 0;		/* Current type (17 or 18)  */
static int4 oochunks  ;		/* Seek position of header  */
static int4 ochunks[2*Nchunks];	/* Table of (Offset, ID)    */
static int  ichunks = 0;	/* Position in above table  */
static int4 total_size ;
static int  missing_USNO;
static int4 acc_index[6000] ;	/* Accelerator index 1/Dacc */
static char *tmp_file = "/tmp/NomaD1";
static int  tmp_hand;		/* Temporary File */

static short xtra2[1000];	/* Table of extra values */
static int4  xtra4[300];	/* Table of extra values */
static int   nxtra2, nxtra4;
static int   txtra2, txtra4;
static int4  usnob1_min[3];
static int4  nucac2 = 0;
static int4  ucac2_lim[2] = { 9999999, 0 };
static int4  ntyc2  = 0;
static int   *tycho2_index;	/* recno tyc, RAmas, SDmas */
static int   ntycho2_index; 
static int   tycho2_mas;	/* Max difference Tycho2 (mas) */

static char *store_msg = 		/* Used by store function */
       "Value truncated,";

static int tycho2_diff[4];
static int tycnum[3];		/* Tycho Number */
static int swapping ;	/* Set to 1 or 2 (swap)	*/
static char verbop = 0;	/* Verbose Option	*/

static char intro[] = "\
An index is made of the 3 parts:\n\
   1) One ascii line starting by \"NOMAD-1.0(25)\" or \"NOMAD-1.0(26)\"\n\
   2) A header table of 80 offsets in the case of \"NOMAD-1.0(36)\"\n\
      which contains the NOMAD1-number of first record in chunk\n\
   3) A set of 80 'chunks', each dealing with about 5degrees\n\
Each 'chunk' has a header made of:\n\
   -- A header of 36 bytes (7int + 4shorts)\n\
      + lenHeader + 36+size(Extra)\n\
      + MIN values of ID RA SD + MAX values of ID RA MAX\n\
      + How many 'extras' in fields / mu's / mag's / e_mu\n\
   -- the 'Extra' table of shorts giving the 'extra' values\n\
   -- An index of (ID,RA,bytePos) for fast position (each 1,000 record)\n\
   -- The data for each star, as 17-(18-)byte records + nphot*7\n\
" ;

static char usage[] = "\
Usage: nomad_make [-v] [-r root_name] [-t tmp_file] input_file [tycho2_index]\n\
";

static char help[] = "\
      -v: verbose option\n\
 -r name: defines the name of Root for output (defaut from current directory)\n\
 -t name: defines the name of the temporary file (default NomaD1)\n\
   input: a \"zoneXXXX.cat\" original file\n\
";

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

static void swap1 (array, nshort)
/*++++++++++++++++
.PURPOSE  Swap the bytes in the array of shorts if necessary
.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
.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 int wr2 (short *array, int nshort)
/*++++++++++++++++
.PURPOSE  Write out an array of short int
.RETURNS  Bytes written out
.REMARKS  Writes in Sun (FITS) order
-----------------*/
{
  int n;
    if (swapping == 0) {	/* I'm on a big-endian machine */
	swapping ^= 2;
	swap1(array, nshort);
	n = fwrite(array, 2, nshort, u1file);
	swap1(array, nshort);
	swapping ^= 2;
    }
    else n = fwrite(array, 2, nshort, u1file);
    return(n);
}

static int wr4 (int4 *array, int nint)
/*++++++++++++++++
.PURPOSE  Write out an array of int
.RETURNS  Bytes written out
.REMARKS  Writes in Sun (FITS) order
-----------------*/
{
  int n ;
    if (swapping == 0) {	/* I'm on a big-endian machine */
	swapping ^= 2;
	swap(array, nint);
	n = fwrite(array, 4, nint, u1file);
	swap(array, nint);
	swapping ^= 2;
    }
    else n = fwrite(array, 4, nint, u1file);
    return(n);
}

static int verify(char *text, int value, int nbits)
/*++++++++++++++++
.PURPOSE  Verify the value has the correct number of bits
.RETURNS  0 (OK) / -1 (bad)
.REMARKS  When text id NULL ==> NO ERROR MESSAGE
-----------------*/
{
  static int maxis[32];
  int i,m;
    if (!maxis[1]) for (i=0, m=1; m; i++, m<<=1) maxis[i] = m-1;
    if (value & (~maxis[nbits])) {
	if (text) {
	    fprintf(stderr, "++++%s too large for %d bits: %d = 0x%x",
	      text, nbits, value, value);
	    if (tycnum) fprintf(stderr, ", TYC %d-%d-%d", 
	      tycnum[0], tycnum[1], tycnum[2]);
	    fprintf(stderr, "\n");
	}
	return(-1);
    }
    else return(0);
}

static int store(unsigned char *a, int b, int value, int len)
/*++++++++++++++++
.PURPOSE  Store the 'value' made of 'len' bits
	  to byte position 'a' offseted 'b' bits.
.RETURNS  The number of bits free in last byte (0 to 7) / <0 for Truncated
.REMARKS  Independent of the architecture. If number too large, returns -1
	*** The unused part of the rightmost byte is set to zero.
	*** Therefore only use from left to right !
-----------------*/
{
  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 } ;
  unsigned char *ac;
  int fb, /* free right bits */
      nb, /*remaining bits to write */
      status;

    /* Test the number is not too large  */
    if (value & (~mask[len])) {
	if (store_msg) fprintf(stderr, 
	  "++++%s too large for %d bits: %d = 0x%x",
	  store_msg, len, value, value);
	if (value < 0) value = 0;
	else value = mask[len];
	if (store_msg) fprintf(stderr, " => %d\n", value);
	status = -1;
    }
    else status = 0;

    ac = a + (b>>3);	/* Byte position  */
    nb = len;		/* Bits to move   */

    if (b&7) { 	  /* Fill incomplete byte */
        fb = 8 - (b&7);	/* Bits available */
        if (fb >= nb) { fb -= nb; *ac |= value<<fb; nb=0; }
        else          { nb -= fb; *ac |= value>>nb; fb=0; }
	ac++;
    }

    /* Fill the complete bytes */
    while (nb >= 8) { nb -= 8; *(ac++) = value>>nb; }

    if (nb) {	/* Last byte: rightmost buts are set to zero */
	fb = 8 - nb;
	*ac = value<<fb;
    }
    else fb = 0;

    return (status|fb);	/* Negative number if error */
}


static int extra_value(int value, int max, int max2, int edit)
/*++++++++++++++++
.PURPOSE  Reserve a cell for a value out-of-bounds, as short int or int4
.RETURNS  [0,max] = standard; >max => saved in table; >max2 => saved in tab4
.REMARKS  Message edited if edit!=0
-----------------*/
{
  int i;
    if ((value >=0) && (value <= max)) return(value);
    if (edit) fprintf(stderr, "....use extra for %s >%d: %d\n", store_msg,
	    max, value);
    if ((value >= -32767) && (value <= 32767)) {
        for (i=0; (i<nxtra2) && (xtra2[i] != value); i++) ;
	if (i >= nxtra2) {
	    if (i >= ITEMS(xtra2)) {
	        fprintf(stderr, "****xtra2 array (%d) too small: %d\n", 
			    ITEMS(xtra2));
		exit(1);
	    }
	    xtra2[nxtra2++] = value;
	}
        value = max+1 + i;
	if (value >= max2) {
	    fprintf(stderr, "****Too many (%d) short extra values (%s)\n", 
		    i, store_msg);
	    exit(1);
	}
    }
    else {
        for (i=0; (i<nxtra4) && (xtra4[i] != value); i++) ;
	if (i >= nxtra4) {
	    if (i >= ITEMS(xtra4)) {
	        fprintf(stderr, "****xtra4 array (%d) too small: %d\n", 
			    ITEMS(xtra4));
		exit(1);
	    }
	    xtra4[nxtra4++] = value;
	}
       	value = max2+1 + i;
    }
    return(value) ;
}

/*==================================================================
		Compute distance between 2 points
 *==================================================================*/
static int dist(int ra1, int de1, int ra2, int de2)
/*++++++++++++++++
.PURPOSE  Cmpute distance between 2 points (mas)
.RETURNS  Difference
-----------------*/
{ 
  double u1[3], u2[3], ra, de;
  int mas;
    ra = ra1; ra *= M_PI/(3.6e6*180.);
    de = de1; de *= M_PI/(3.6e6*180.); de -= M_PI/2.;
    u1[0] = cos(de); u1[2] = sin(de);
    u1[1] = u1[0] * sin(ra); u1[0] *= cos(ra);
    ra = ra2; ra *= M_PI/(3.6e6*180.);
    de = de2; de *= M_PI/(3.6e6*180.); de -= M_PI/2.;
    u2[0] = cos(de); u2[2] = sin(de);
    u2[1] = u2[0] * sin(ra); u2[0] *= cos(ra);
    /* Difference = 4 * sin^2(r/2) */
    u1[0] -= u2[0];
    u1[1] -= u2[1];
    u1[2] -= u2[2];
    de = sqrt(u1[0]*u1[0] + u1[1]*u1[1] + u1[2]*u1[2]);
    de = 2.*asin(de/2.);
    de *= (3.6e6*180.)/M_PI;
    mas = de+0.5;
    return(mas);
}

/*==================================================================
		Load the USNO range
 *==================================================================*/
static void load_usno(char *fname, int zone)
/*++++++++++++++++
.PURPOSE  Load the "usnob1.counts" file
.RETURNS  ---
.REMARKS  Set the usnob1_min array (0=current, 1=nex, 2=prev)
-----------------*/
{ 
  char buf[128], *p;
  FILE *f;
  int z;
    memset(usnob1_min, 0, sizeof(usnob1_min));
    f = fopen(fname, "r");
    if (!f) {
	fprintf(stderr, "****load_usno: ");
	perror(fname);
	exit(1);
    }
    while (fgets(buf, sizeof(buf), f)) {
	if (buf[0] != 'z') continue;
	if (atoi(buf+2) == zone) break;
        if (p = strchr(buf, '+')) 
            usnob1_min[2] = atoi(p+1);	/* Preceding zone */
    }
    if (p = strchr(buf, '+')) 
        usnob1_min[0] = atoi(p+1);	/* Current zone   */
    if (fgets(buf, sizeof(buf), f)) {
	if (p = strchr(buf, '+')) 
	    usnob1_min[1] = atoi(p+1);	/* Next zone */
    }
    else {	/* PROBLEM for the 1799 zone !!! 
		   In NOMAD file, starts from 1 ... */
	fprintf(stderr, "++++Very last zone, restart from 1\n");
	usnob1_min[1] = usnob1_min[0];
	usnob1_min[0] = 0;
    }

    fprintf(stderr, "....load_usno(%s): usnob1_min(%d,%d,%d)=(%d,%d,%d)\n",
	    fname, zone, zone+1, zone-1, 
	    usnob1_min[0], usnob1_min[1], usnob1_min[2]);
    fclose(f);
}

/*==================================================================
		Load one Tycho Index
 *==================================================================*/
static int diff_int(int *i1, int *i2) 
/*++++++++++++++++
.PURPOSE  Sort function
.RETURNS  Difference
-----------------*/
{ 
    return(*i1 - *i2);
}

static void load_tycho2(char *fname)
/*++++++++++++++++
.PURPOSE  Load a Tycho index (generated by get_tycho2)
.RETURNS  ---
.REMARKS  Creates the tycho2 = (recno, 14+14+3, RA, SD)
-----------------*/
{
  int val, lines, alines, i;
  char buf[128], *p;
  FILE *f;

    f = *fname == '-' ? stdin : fopen(fname, "r");
    if (!f) {
	fprintf(stderr, "****load_tycho2: ");
	perror(fname);
	exit(1);
    }

    /* Read line-by-line */
    lines = alines = 0;
    while (fgets(buf, sizeof(buf), f)) {
	if (lines <= alines) {
	    alines += 250;
	    if (!tycho2_index) tycho2_index = malloc(16*alines);
	    else tycho2_index = realloc(tycho2_index, 16*alines);
	}
	i = lines*4;
	p = buf;
	tycho2_index[i++] = atoi(p);	/* recno */
	while (isgraph(*p)) p++; while (isspace(*p)) p++;

	tycnum[0] = atoi(p);
	while (isdigit(*p)) p++; if (*p == '-') p++;
	tycnum[1] = atoi(p);
	while (isdigit(*p)) p++; if (*p == '-') p++;
	tycnum[2] = atoi(p);
	while (isdigit(*p)) p++; while (isspace(*p)) p++;
	/* Verify components of TYC2 */
	verify("TYC1", tycnum[0], 14);
	verify("TYC2", tycnum[1], 14);
	if ((tycnum[2]<1) || (tycnum[2]>5)) {
	    fprintf(stderr, "****Invalid Tycho-2 Number (load_tycho2)\n\t%s",
		    buf);
	    exit(2);
	}
	val = (tycnum[0]<<17) | (tycnum[1]<<3) | tycnum[2];
	tycho2_index[i++] = val;	/* Tycho-2 Name */

	tycho2_index[i++] = atoi(p);	/* RAmas */
	while (isdigit(*p)) p++; while (isspace(*p)) p++;
	tycho2_index[i++] = atoi(p);	/* SDmas */
	lines++;
    }
    if (f != stdin) fclose(f);

    fprintf(stderr, "....load_tycho2: %d TYC2 stars\n", lines);
    ntycho2_index = lines;

    /* Sort it ! */
    qsort(tycho2_index, lines, 16, (SORT_FCT)diff_int);
    
#if 0
    fprintf(stderr, "====List of Tycho2 stars:\n");
    for (i=0; i<ntycho2_index*4; i+=4) 
	fprintf(stderr, "%11d %08X %10d%09d\n",
	    tycho2_index[i],   tycho2_index[i+1],
	    tycho2_index[i+2], tycho2_index[i+3]);
#endif

    /* To get the record: 
     * bsearch(&recno, tycho2_index, lines, 16, diff_int) 
     */
}

/*==================================================================
		Convert the Input Record(s)
 *==================================================================*/

static int trec(int4 rec[22], unsigned char *out, int flag)
/*++++++++++++++++
.PURPOSE  Convert the input (RA+DE+MAG) into 36- or 8-byte record
.RETURNS  Number of bytes in out / -1=Error
.REMARKS  flag=0 ==> reclen=26 (NOMAD1reclen+1)
-----------------*/
{
  static int4 id = 0;
  unsigned char *o, *o23;
  int i, b, mas;
  int4 word;

    o = out;
    o23 = out+23;	/* Where to write the 2MASS/YB6 flags */
    b = 0;		/* Bit offset */
    *(o++) = 0;		/* Reserved for flags */

    /* The USNO-B1 number, should be OK with 22 bits */
    if (rec[16] == 0) {
	*(o++) = 0;
	*(o++) = 0;
	*(o++) = 0;
	missing_USNO ++;
    }
    else {
	if (rec[16]>usnob1_min[1]) i=1;		/* Upper zone */
	else if (rec[16]>usnob1_min[0]) i=0;	/* Same  zone */
	else i = 2;				/* Lower zone */
        word = rec[16]-usnob1_min[i];
	store_msg = "USNO-B1";
	store(o, b, i, 2); b += 2;
	store(o, b, word, 22); b += 22;
	o += b/8;  b = 0;
    }

    word = rec[0];	/* RA */
    *(o++) = word>>24;
    if (flag) --o; 	/* Only 3 bytes for RA */
    else ++o23;
    *(o++) = word>>16;
    *(o++) = word>> 8;
    *(o++) = word;

    store_msg = "  e_RA";
    word = rec[2];	/* sigma(RA), 10b */
    store(o, b, word, 10); b += 10;

    /* Epoch: round to 0.1yr since 1901, 10b */
    store_msg = "  EpRA";
    word = ((rec[8]+50)/100) - Oep;
    store(o, b, word, 10); b += 10;

    word = rec[1] - Dstep*SDzone;	/* SPD */
    if ((word<0) || (word>Dstep)) {
	fprintf(stderr, "****SPD is outside range: %d\n", word);
        exit(1);
    }
    store_msg = "   SPD"; 
    store(o, b, word, 20); b += 20;

    if (b&7) {
	fprintf(stderr, "****Should be a byte boundary: b\n", b);
	exit(1);
    }
    o += b/8; b = 0;

    store_msg = "  e_DE";
    word = rec[3];	/* sigma(DE), 10b */
    store(o, b, word, 10); b += 10;

    /* Epoch: round to 0.1yr since 1901, 9b */
    store_msg = "  EpDE";
    word = ((rec[9]+50)/100) - Oep;
    store(o, b, word, 10); b += 10;

    /* pm(RA) + 7500, 14b   */
    store_msg = "  pmRA";
    word = extra_value(rec[4]-Opm, EXTRA_14_2, EXTRA_14_4, verbop);
    store(o, b, word, 14); b += 14;

    /* pm(DE) + 7500, 14b   */
    store_msg = "  pmDE"; 
    word = extra_value(rec[5]-Opm, EXTRA_14_2, EXTRA_14_4, verbop);
    store(o, b, word, 14); b += 14;

    /* sigma(pmRA), 13b */
    store_msg = "e_pmRA";
    word = extra_value(rec[6], EXTRA_13_2, EXTRA_13_4, 1);
    store(o, b, word, 13); b += 13;

    /* sigma(pmDE), 13b */
    store_msg = "e_pmDE";
    word = extra_value(rec[7], EXTRA_13_2, EXTRA_13_4, 1);
    store(o, b, word, 13); b += 13;

    /* Write out the flags -- keep the net-order */
    word = rec[21];
    store_msg = " flags";
    if (word&0xc000c000) fprintf(stderr, 
	"****Status has reserved bits set: %x08X\n", word);
    store(o, b, word, 30); b += 30;

    if (b&7) {
	fprintf(stderr, "****SHOULD be a byte boundary: b\n", b);
	exit(1);
    }
    o += b/8; b = 0;

    /* The 2MASS number -- replace by a simple bit */
    if (rec[17]) *o23 |= 0x80;

    /* The YB6 number --  replace by a simple bit */
    if (rec[18]) *o23 |= 0x40;

    /* The UCAC2 number -- 24 bits, offset *ucac2 */
    if (rec[19]) {
	nucac2++;
	word = rec[19];
	out[0] |= 0x80;
	if (word < ucac2_lim[0]) ucac2_lim[0] = word;
	if (word > ucac2_lim[1]) ucac2_lim[1] = word;
	*(o++) = word>>24;
	*(o++) = word>>16;
	*(o++) = word>>8;
	*(o++) = word;
    }

    /* The Tycho-2 Number -- save as 14b+14b+3b -- can be Hipparcos ! */
    if (rec[20]) {
	ntyc2++;
	word = mas = 0;
	if (tycho2_index) {	/* Convert to Tycho-ID */
	  int *atyc, da, dd;
	    atyc = bsearch(rec+20, 
		   tycho2_index, ntycho2_index, 16, (SORT_FCT)diff_int);
	    if (atyc) {
		word = atyc[1];	
		da = rec[0] - atyc[2]; 
		dd = rec[1] - atyc[3];
		if ((da>=100) || (da<=-100) || (dd>=100) || (dd<=-100))
		    mas = dist(rec[0], rec[1], atyc[2], atyc[3]);
		else mas = 0;
		if (mas > tycho2_mas) tycho2_mas = mas;
		if (mas>=300) fprintf(stderr, 
		    "++++Difference NOMAD/Tycho-2 (mas): %5d (%+4d %+4d)\n", 
		    mas, da, dd);
		if (mas >= 3000) mas = 3;
		else if (mas >= 1000) mas = 2;
		else if (mas >= 300) mas = 1;
		else mas = 0;
		tycho2_diff[mas] += 1;
	    }
	    else { double ra, dec;
		ra  = rec[0]; ra /= 3.6e6;
		dec = rec[1]; dec = (dec/3.6e6) - 90.;
		if (rec[21]&NOMAD_HIPAST) ;	/* Hipparcos astrometry */
		else fprintf(stderr, 
			"++++Tycho-2 missing: %d (%09.5f%+09.5f)\n", 
			rec[20], ra, dec);
		word = rec[20]<<3;		/* Room for TYC3 */
	    }
	}
	else word = rec[20]<<2;
	out[0] |= 0x40;

	store_msg = "Tyc-diff"; store(o, b,  mas/2, 1);  b += 1;
	store_msg = "Tyc-num."; store(o, b, word, 31); b += 31;
	o += b/8; b = 0;
    }

    /* The 6 pass-bands in photometry */
    for (i=10, b=0x20; i<16; i++, b>>=1) {
	if (rec[i] >= 30000) continue;	/* Inexistent magnitude */
	word = rec[i]-Omag;
	if (word&0xffff0000) {
	    fprintf(stderr, "****Mag#%d out-of-bounds: %d\n", word);
	    exit(1);
	}
	*(o++) = word>>8;
	*(o++) = word;
	out[0] |= b;
    }

    /* The length in bytes */
    return(o-out) ;
}

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

static int4 store_chunk(int RAzone, int id0, int id1, int tmp_hand)
/*++++++++++++++++
.PURPOSE  The file tmp_hand contains a chunk -- write it.
.RETURNS  Number of Objects in the file
.REMARKS  Output the offset position
-----------------*/
{
  int4 bytes, minmax[3], filepos[3];
  char buffer[8192];
  short extra_counts[2];
  int4 length, nobj, i, xalign, iacc;

    nobj = id1-id0+1;
    if (nobj< 100) printf(
	"++++Only %d objects in chunk %04d#%03d\n",
	nobj, SDzone, RAzone);
    printf("----Zone %04d, RA#=%03d: %8ldobj.", 
	SDzone, RAzone, nobj);

    /* Write the text */
    sprintf(buffer, chunk_text, RAzone, id0, id1);
    fwrite(buffer, strlen(buffer), 1, u1file);
    if (strlen(buffer)&3) {
	fprintf(stderr, 
	    "****store_chunk: length(%s) not multiple of 4?\n",
	    buffer);
	exit(1);
    }
    fflush(u1file);
    filepos[0] = ftell(u1file);

    /* Save the position + first ID of this CHUNK. */
    ochunks[ichunks++] = filepos[0];
    ochunks[ichunks++] = id0;

    /* Count the number of 'extras' */
    extra_counts[0] = nxtra2;
    extra_counts[1] = nxtra4;
    length = nxtra2*2 + nxtra4*4;
    printf(" Extras: short=%-3d int4=%d", nxtra2, nxtra4);
    printf(" (%d)\n", length);

    /* Write the header length */
    length += 32;
    xalign = length&2 ;
    if (xalign) length += 2; 
    wr4(&length, 1);

    /* Write the 3 minima, maxima, number of extras */
    minmax[0] = id0; 
    minmax[1] = RAzone<<24; 
    minmax[2] = SDzone*360000;
    wr4(minmax, 3);
    minmax[0] = id1; 
    if (reclen == NOMAD1reclen) 
	 minmax[1] |= 0xffffff;
    else minmax[1] = 360*3600*1000 -1;
    minmax[2] += 360*1000;	/* 0.1deg */
    wr4(minmax, 3);
    wr2(extra_counts, 2);

    /* Write all 'extras' */
    if (nxtra4>0) wr4(xtra4, nxtra4);
    if (nxtra2>0) wr2(xtra2, nxtra2);

    /* Align the extras   */
    if (xalign) { xalign = 0; fwrite(&xalign, 1, 2, u1file); }

    /* Complete the Accelerator */
    bytes = lseek(tmp_hand, 0L, 1);	/* Number of Data Bytes */
    lseek(tmp_hand, 0L, 0) ;

    /* How many elements in the accelerator? Count on [1] = ID */
    for (iacc=1; acc_index[iacc]; iacc += 3) ;
    /* No need of an extra index when just a few remain */
    if ((iacc>4) && ((id1-acc_index[iacc-3]) < (Dacc/3)))  
	iacc -= 3;
    acc_index[--iacc] = bytes;		/* Position of End */
    acc_index[++iacc] = id1+1;		/* Last ID...	   */
    acc_index[++iacc] = 0x7fffffff;	/* Ending RA!!	   */
    length += 4*(++iacc);		/* length now includes Accelerator */
    for (i=0; i<iacc; i += 3) acc_index[i] += length;

    /* Write the Accelerator */
    wr4(acc_index, iacc);

    /* The chunk must be a multiple of 4 bytes */
    xalign = bytes&3; if (xalign) xalign = 4-xalign;
    filepos[1] = ftell(u1file);		/* For Debugging */
    /* Copy now the temporary file */
    while (bytes > 0) {
	i = sizeof(buffer);
	if (i > bytes) i = bytes;
	i = read(tmp_hand, buffer, i);
	if (i<=0) fprintf(stderr,
	    "****Probleme re-reading temporary file ??\n") ;
	else bytes -= i;
	fwrite(buffer, 1, i, u1file);
    }
    /* Align... */
    memset(buffer, 0, 4);
    if (xalign) fwrite(buffer, 1, xalign,  u1file);
    filepos[2] = ftell(u1file);		/* For Debugging */

    /* Reinitialize */
    lseek(tmp_hand, 0L, 0) ;
    memset(xtra2, 0, sizeof(xtra2));
    memset(xtra4, 0, sizeof(xtra4));
    txtra2 += nxtra2;
    txtra4 += nxtra4;
    nxtra2 = nxtra4 = 0;
    memset(acc_index, 0, sizeof(acc_index));

    return(nobj);
}

/*==================================================================
		Main Program
 *==================================================================*/
main (argc, argv) int argc; char **argv;
{
  static char header_line[Lheader];
  int i, k, raz, razp, RAmask, id0, id, iacc;
  char *p, *root;
  char name[40] ;
  int4 rec22[22], bytes, inbytes, total, nobj;
  unsigned char recb[60] ;
  double mean;

    /* Check the Arguments */
    total = 0 ;
    root = "." ;
    raz = razp = iacc = 0;
    while (-- argc > 0) {
	p = *++argv;
	if (*p == '-') switch (*++p) {
	  case 'h':
	    printf("%s%s\n%s", usage, help, intro) ;
	    exit(0);
	  case 'r':
	    root = *++argv; -- argc;
	    continue;
	  case 't':
	    tmp_file = *++argv; -- argc;
	    continue;
	  case 'v':		/* Verbose		*/
	    verbop = 1; 
	    continue;
	  default:
	  bad_option:
	    fprintf(stderr, "\n****Bad option: -%s\n", p);
	    exit(1);
	}
	break ;
    }

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

    if (argc == 2) {		/* LOAD Tycho-2 Index */
	load_tycho2(argv[1]);
    }

    /* Open the temporary file */
    tmp_hand = open(tmp_file, O_TRUNC|O_CREAT|O_RDWR, -1);
    if (tmp_hand < 0) {
	perror(tmp_file);
	exit(1);
    }

    /* Verify the Swapping ! */
    { static int value ; char *v ;
        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) ;
        }
    }
    swapping ^= 2;		/* Original NOMAD is big-endian */
    fprintf(stderr, "#...swapping type=%d\n", swapping) ;

    p = *argv;			/* File Name 	*/
    p += strlen(p) - 8;
    SDzone = atoi(p) ;		/* Zone in 0.1deg	*/
    orifile = p - 5 ;		/* Name original file	*/

    load_usno("usnob1.counts", SDzone);

    infile = fopen(*argv, "rb");
    if (!infile) { perror(*argv); exit(1); }
    fseek(infile, 0, 2);	/* Size in bytes */
    inbytes = ftell(infile);
    fseek(infile, 0, 0);
    /* < 8 Mbytes is a single chunk */
    reclen = NOMAD1reclen;
    if (inbytes < 8*1024*1024) reclen++;

    /* Define Output Names */
    if (root && chdir(root)) { perror(root); exit(1); }
    strcpy(name, orifile); name[4] = 'N';
    p = strchr(name, '.'); strcpy(p, ".bin");
    /* Go to the subdirectory */
    name[3] = 0;
    if (mkdir(name, -1)) {
	if (errno == EEXIST) ;
	else { fprintf(stderr, "****"); perror(name); exit(1); }
    }
    name[3] = '/';

    /* Open the output file dealing with this 0.1 degrees zone */
    u1file = fopen(name, "wb") ;
    if (! u1file) {
    	perror(name) ;
    	exit(-1) ;
    }

    /* Write the header text, must be a multiple of 8 bytes */
    sprintf(header_line, header_text, reclen, name, 0, 0, 
		    Opm, Omag, Oep, 0, 0, 0, 0, 0, 0);
    for (i=strlen(header_line); i<sizeof(header_line); i++) 
	header_line[i] = ' ';
    header_line[sizeof(header_line)-1] = '\n';
    fwrite(header_line, 1, sizeof(header_line), u1file);

    /* At this file position the List of Chunks will be written.
       Keep this position, will be rwritten when all chunks written
    */
    fflush(u1file); oochunks = ftell(u1file);
    wr4(ochunks, ITEMS(ochunks));
    fflush(u1file);

    /* Change of Chunk depends on reclen. */
    RAmask = 0x7ff << 24 ;
    if (reclen > NOMAD1reclen) RAmask = 0;

    /* Loop reading the input file, sorted in RA */
    id0 = 1; id = 0;
    while((i = fread(rec22, sizeof(int4), 22, infile)) > 0) {
	swap(rec22, 22);
    	raz = rec22[0] & RAmask ;
    	if (raz != razp) {	/* A new zone appeared there */
	    nobj = store_chunk(razp>>24, id0, id, tmp_hand);
	    total += nobj ;
	    razp = raz ; id0 = id+1;
	    iacc = 0;
	}
	id++;
    	k = trec(rec22, recb, RAmask) ;
	if (k < 0) {		/* Doesn't work !! */
	    exit(1);
	}
	if (((id-id0)%Dacc) == 0) {
	    acc_index[iacc++] = lseek(tmp_hand, 0L, 1);
	    acc_index[iacc++] = id;
	    acc_index[iacc++] = rec22[0];	/* RA of this entry */
	    if (iacc >= ITEMS(acc_index)-3) {
		fprintf(stderr, 
		  "****acc_index too small (must be >%d)\n", iacc);
		exit(1);
	    }
	}
    	write(tmp_hand, recb, k); 
    }
    if (i < 0) perror(*argv) ;

    /* Terminate last Chunk */
    nobj = store_chunk(razp>>24, id0, id, tmp_hand);
    total += nobj ;

    /* Indicates the End of Last Chunk (end-of-file) */
    ochunks[ichunks++] = ftell(u1file);
    ochunks[ichunks++] = id+1;

    /* Update the Header Text */
    fflush(u1file);
    bytes = ftell(u1file); 
    fseek(u1file, 0, 0) ;
    sprintf(header_line, header_text, reclen, name, 1, id, 
		    Opm, Omag, Oep, txtra2, txtra4, 
		    nucac2, ucac2_lim[0], ucac2_lim[1], ntyc2);
    if (header_line[sizeof(header_line)-1] != '\n') {
	printf("****Header line too long:\n%s\n", header_line);
	exit(1);
    }
    header_line[strlen(header_line)] = ' ';
    fwrite(header_line, 1, sizeof(header_line), u1file);

    /* Rewrite the List of Chunks (was just zeroes) */
    if (oochunks) {
        fseek(u1file, oochunks  , 0) ;
	wr4(ochunks, ITEMS(ochunks));
	if (verbop) fprintf(stderr, "    (header rewritten)\n") ;
    }
    
    /* Edit some statistics */
    printf(
	"====Tycho-2 connections [0.3,1,3]\": %d %d %d %d - max=%dmas\n",
	tycho2_diff[0], tycho2_diff[1], tycho2_diff[2], tycho2_diff[3], 
	tycho2_mas);
    fprintf(stderr, 
	"====Tycho-2 connections [0.3,1,3]\": %d %d %d %d - max=%dmas\n",
	tycho2_diff[0], tycho2_diff[1], tycho2_diff[2], tycho2_diff[3], 
	tycho2_mas);
    fprintf(stderr, 
	"    IDs: UCAC2=%d, Tyc2=%d; noUSNO=%d; xtra=%d,%d\n",
	nucac2, ntyc2, missing_USNO, txtra2, txtra4);
    printf("    IDs: UCAC2=%d, Tyc2=%d; noUSNO=%d; xtra=%d,%d\n",
	nucac2, ntyc2, missing_USNO, txtra2, txtra4);
    mean = bytes; mean /= total;
    fprintf(stderr, "====> %8ld records in %10ld bytes, %.1fb/o", 
	    total, bytes, mean);
    printf("====> %8ld records in %10ld bytes, %.1fb/o", 
	    total, bytes, mean);

    /* Compute the compression factor */
    mean = inbytes-bytes;
    mean = 100.*mean/inbytes;
    printf(" (reduction=%.0f%%)\n", mean);
    fprintf(stderr, " (reduction=%.0f%%)\n", mean);
    exit(0);
}