/*++++++++++++++
.IDENTIFICATION usnob_make.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-A2.0 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   09-Nov-2002
.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 <usnob1def.h>
---------------*/

#include <usnob1def.h>	/* Structure definitions */
#include <usnob1.h>	/* Structure definitions */

#include <stdio.h>
#include <stdlib.h>	/* Malloc   */
#include <string.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]))

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 long oochunks  ;		/* Seek position of header  */
static long ochunks[2*Nchunks];	/* Table of (Offset, ID)    */
static int  ichunks = 0;	/* Position in above table  */
static long total_size ;
static long acc_index[6000] ;	/* Accelerator index 1/Dacc */
static char *tmp_file = "/tmp/UsnoB1";
static int tmp_hand;		/* Temporary File */

static short tf[300];		/* List of Fields */
static short tpm[2000];		/* Extra pm's */
static short tmag[2000];	/* Extra mag's */
static short txi[2000];		/* Extra xi's */
static int   htf[300];		/* Counts of Extras */
static int   htpm[2000];	/* Extra Counts */
static int   htmag[2000];	/* Extra Counts */
static int   htxi[2000];	/* Extra xi's */
static int   hlenph[8];		/* How many photometries of each sort */

static char *store_msg = 		/* Used by store function */
       "Value truncated,";
static short *xtratabs[] = { tf, tpm, tmag, txi } ;

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 \"USNO-B1.0(17)\" or \"USNO-B1.0(18)\"\n\
   2) A header table of 128 offsets in the case of \"USNO-B1.0(17)\"\n\
   3) A set of 128 '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 / xi-eta \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: usnob_make [-v] [-r root_name] [-t tmp_file] input_file\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 UsnoB1)\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
.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 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 of %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
-----------------*/
{
  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 of %d bits: %d = 0x%x",
	  store_msg, b, value, value);
	value &= mask[len];
	status = -1;
    }
    else status = 0;

    ac = a + (b>>3);	/* Byte position  */
    nb = b;		/* 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) {	/* Fill last incomplete byte */
	fb = 8 - nb;
	*ac |= value<<fb;
    }
    else fb = 0;

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


static int extra_value(int value, int max, short *table, int *counts)
/*++++++++++++++++
.PURPOSE  Reserve a cell for a value out-of-bounds
.RETURNS  [0,max] = standard; >max => saved in header
-----------------*/
{
  short *s ;
    if ((value >=0) && (value <= max)) return(value);
    for (s=table; (*s != value) && *s ; s++ ) ;
    if (!*s) *s = value ;
    if (counts) { int i; i=s-table; counts[i] += 1; }
    value = max+1 + (s-table) ;
    return(value) ;
}

static int extra_count(short *table)
/*++++++++++++++++
.PURPOSE  Count the number of values
.RETURNS  The number
-----------------*/
{
  short *s ;
    for (s=table; *s ; s++ ) ;
    return(s-table) ;
}

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

static int photostore(long *rec, unsigned char *out, int bytes)
/*++++++++++++++++
.PURPOSE  Code one Photometry on 2, 6 or 7 bytes
.RETURNS  0 (OK) / -1 (failed)
.REMARKS  We assume the mag is not NULL...
-----------------*/
{
  long word;
  int mag, stargal, sgcs, survey, field, xi, eta;

    word = rec[0];
    mag     = word%10000; word/=10000;	/* The value of the magnitude  */
    field   = word%1000; word /= 1000;
    survey  = word%10; word /= 10; 
    stargal = word;
    if (stargal >= 20) fprintf(stderr, "****stargal=%d !!\n", stargal);

    if (bytes == 2) {			/* Just store the Magnitude */
	out[0] = mag>>8;
	out[1] = mag;
	return(0);
    }

    word = rec[5]; 
    xi   = word%10000-5000; word /= 10000;
    eta  = word%10000-5000; word /= 10000;
    sgcs = stargal*100 + word*10 + survey;
    verify("S/G+C,S", sgcs, 11);

    /* Move to the appropriate ranges */
    mag = extra_value(mag-Omag, Nmag, tmag, htmag);
    xi  = extra_value( xi-Oxi , Nxi,  txi,  htxi);
    eta = extra_value(eta-Oxi , Nxi,  txi,  htxi);
    field = extra_value(field, 0, tf, htf);

    out[0]  = sgcs>>3;
    out[1]  = sgcs<<5;

    if (bytes == 6) {			/* Verify Coding on 6 bytes */
	if (verify((char *)0, mag, 11)) return(-1);
	if (verify((char *)0,  xi, 11)) return(-1);
	if (verify((char *)0, eta, 11)) return(-1);
	if (verify((char *)0, field,4)) return(-1);

	out[1] |= field<<1;
	out[1] |= mag>>10;
	out[2]  = mag>>2;
	out[3]  = mag<<6;
	out[3] |= xi>>5;
	out[4]  = xi<<3;
	out[4] |= eta>>8;
	out[5]  = eta;

	return(0);
    }

    /* Coding on 7 bytes */
    if (
        verify("--mag", mag, 13) ||
        verify("---xi",  xi, 13) ||
        verify("--eta", eta, 13) ||
        verify("--fld", field,5)
    ) exit(1);
    out[1] |= field>>1;
    out[2]  = field<<7;
    out[2] |= mag>>6;
    out[3]  = mag<<2;
    out[3] |=  xi>>11;
    out[4]  =  xi>>3;
    out[5]  =  xi<<5;
    out[5] |= eta>>8;
    out[6]  = eta;

    return(0) ;
}

static int trec(long rec[20], unsigned char *out, int flag)
/*++++++++++++++++
.PURPOSE  Convert the input (RA+DE+MAG) into 7- or 8-byte record
.RETURNS  Number of bytes in out / -1=Error
.REMARKS  flag=0 ==> reclen=18
-----------------*/
{
  static long int id = 0;
  unsigned char *head, *o, *uflags, *arec;
  long word;
  int i, m, val, lenph;

    o = out;
    head = o++;		/* Where to write the Flags */
    *head = 0;
    tycnum[0] = 0;
    id++;
    if (rec[3] < 100000000) {	/* A TYCHO Record */
	*head = USNOB_TYC;
	lenph = 2;
    }
    else lenph = 6;

    /* Positions: Right Ascension */
    arec = o;		/* Where to write number of detections */
    word = rec[0];	/* RA in 10mas */
    *(o++) = word>>24;
    if (flag) --o;	/* Only 17 bytes */
    *(o++) = word>>16;
    *(o++) = word>> 8;
    *(o++) = word;

    /* South Polar Distance */
    word = rec[1] - (SDzone*36000);
    if ((word<0) | (word>36000)) {
	fprintf(stderr, "****Incompatible SD: %d => offset=%d\n", word, rec[1]);
	return(-1);
    }
    *(o++) = word>> 8;
    *(o++) = word;

    /* Mean Errors on Position, or Tycho Number */
    if (*head) {
	*arec = *arec&0xf;		/* No detection value */
	tycnum[0] = rec[10];
	tycnum[1] = rec[11]/10;
	tycnum[2] = rec[11]%10;
	verify("TYC1", tycnum[0], 15);
	verify("TYC2", tycnum[1], 15);
	*(o++) = tycnum[0]>> 8;
	*(o++) = tycnum[0];
	*(o++) = tycnum[1]>>8;
	*(o++) = tycnum[1];
	*(o++) = tycnum[2];
    }
    else {				/* keeevvvuuu */
	word = rec[4];
	val = word%1000; word /= 1000;	/* e_RA, 10bits */
	*(o++) = val>>2;
	*o     = val<<6;
	val = word%1000; word /= 1000;	/* e_SD, 10bits */
	*(o++)|= val>>4; 		/* (6b) */
	*o     = val<<4;
	uflags = o;			/* --Position to store USNOB_flags */
	val = word%1000; word /= 1000;  /* Epoc */
	verify("Epoch", val, 9);
	*(o++)|= val>>8; 		/* Keep a single bit */	
	*(o++) = val;
	if (word) *uflags |= USNOB_YS4;

	word = rec[3];			/* jMRQyyyxxx */
	word /= 1000000;		/* Skip the Proper Motions */
	val = word%10; word /= 10;	/* fit_RA, 4bits */
	*o  = val<<4;
	val = word%10; word /= 10;	/* fit_SD, 4bits */
	*(o++) |= val;
	val = word%10; word /= 10; 	/* ndet   */
	*arec = (*arec&0xf) | val<<4;
	if (word) *uflags |= USNOB_SPK;
    }

    /* Proper Motions  iPSSSSAAAA */
    word = rec[2];
    val = word%10000;  word /= 10000;	/* pmRA */
    val -= (5000+Opm);			/* Code until ±1arc/yr */
    val = extra_value(val, Npm, tpm, htpm);
    verify("pmra", val, 12);
    *(o++) = val>>4;
    *o     = val<<4;
    val = word%10000;  word /= 10000;	/* pmSD */
    val -= (5000+Opm);			/* Code until ±1arc/yr */
    val = extra_value(val, Npm, tpm, htpm);
    verify("pmsd", val, 12);
    *(o++)|= val>>8; 
    *(o++) = val; 
    if (*head) {			/* TYCHO: Take File + Seq */
	val = rec[3];			/* File Number */
	verify("fileno", val, 5);
	*o = val<<2;
	val = rec[4];			/* Sequence Number */
	if (val < 0) {			/* Strange ... */
	    val = -val;
	    *o |= 0x80;
	}
	verify("recno", val, 18);
	*(o++) |= (val>>16)&3; 
	*(o++)  = val>> 8;
	*(o++)  = val    ;
    }
    else {				/* MUPROB. + e_pm */
	val = word%10; word /= 10;
	*o  = val<<4;
	if (word)			/* USNOB_PM */
	    *uflags |= USNOB_PM;

	/* Errors on Proper Motions */
	word = rec[3];
	val  = word%1000; word /= 1000;	/* e_pmra */
	*(o++) |= val>>6; 
	*o   = val<<2;
	val  = word%1000; word /= 1000;	/* e_pmsd */
	*(o++) |= val>>8; 
	*(o++)  = val; 
    }

    /* Code now the Magnitudes */
    arec = o;				/* Beginning of Photometries */
  Enlarge_Phot:
    for (i=0, m=16; i<5; i++, m >>= 1) {
	if (rec[5+i] == 0)	/* EMPTY Magnitude !!!  */
	    continue ;
	*head |= m;		/* Set the flag indicating a magnitude */
	if (photostore(rec+5+i, o, lenph) == 0) {
	    o += lenph;
	    continue;
	}
	/* ERROR, the storage is too small */
	o = arec;
	if (lenph == 6) { 
	    lenph++; 
	    *head |= USNOB_ph7 ;
	    goto Enlarge_Phot; 
	}
	fprintf(stderr, "****Problem encoding photometry at id=%d\n",
	  id);
	exit(1);
    }
    /* Update the histograms on photometry length */
    for (m=16; m; m >>= 1) {
        if (*head&m) hlenph[lenph] += 1;
    }

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

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

static long 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
-----------------*/
{
  long bytes, minmax[3], filepos[3];
  char buffer[8192];
  short fpmx[4];
  int length, nobj, i, xalign, iacc;

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

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

    /* Count the number of 'extras' */
    printf(" Extras(f,pm,mag,xi): ");
    for (i=0, length=0; i<4; i++) {
	fpmx[i] = extra_count(xtratabs[i]) ;
	if (i) putchar(',');
	printf("%d", fpmx[i]);
	length += 2*fpmx[i];
    }
    printf(" (%d)\n", length);

    /* When there are many additional magnitudes (>=100), 
       list them in the "/tmp/histomag" file
    */
    if (verbop && (fpmx[2] >= 100)) {
	for (i=0; i<fpmx[2]; i++) {
	    printf("....Zone %04d, RA#=%03d magCount=%4d mag=%4d\n",
		SDzone, RAzone, htmag[i], tmag[i]);
	}
    }

    /* Write the header length */
    length += 36;
    xalign = length&2 ;
    if (xalign) length += 2; 
    fwrite(&length, 4, 1, u1file);

    /* Write the 3 minima, maxima, number of extras */
    minmax[0] = id0; minmax[1] = RAzone<<20; minmax[2] = SDzone*36000;
    fwrite(minmax, 4, 3, u1file);
    minmax[0] = id1; 
    if (reclen == 17) minmax[1] += (1<<20)-1;
    else minmax[1] = (1<<30);
    if (minmax[1] > 360*3600*100) minmax[1] = 360*3600*100;
    minmax[2] += 36000;
    fwrite(minmax, 4, 3, u1file);
    fwrite(fpmx, 2, 4, u1file);

    /* Write all 'extras' */
    for (i=0; i<4; i++) 
	fwrite(xtratabs[i], 2, fpmx[i],  u1file);
    /* Align the extras   */
    if (xalign) { xalign = 0; fwrite(&xalign, 2, 1, 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 */
    fwrite(acc_index, 4, iacc, u1file);

    /* 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(tf, 0, sizeof(tf));
    memset(tpm, 0, sizeof(tpm));
    memset(txi, 0, sizeof(txi));
    memset(tmag, 0, sizeof(tmag));
    memset(htf, 0, sizeof(htf));
    memset(htpm, 0, sizeof(htpm));
    memset(htxi, 0, sizeof(htxi));
    memset(htmag, 0, sizeof(htmag));
    /* memset(hlenph, 0, sizeof(hlenph)); */
    memset(acc_index, 0, sizeof(acc_index));

    return(nobj);
}

/*==================================================================
		Main Program
 *==================================================================*/
main (argc, argv) int argc; char **argv;
{
  static char blanks[] = "       \n" ;
  int i, k, raz, razp, RAmask, id0, id, iacc;
  char *p, *root;
  char name[40] ;
  long rec80[20], bytes, inbytes, total, nobj;
  unsigned char recb[18+5*7] ;
  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) {
    	fprintf(stderr, "%s", usage) ;
    	exit(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 USNO-B 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	*/

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

    /* Define Output Names */
    if (root && chdir(root)) { perror(root); exit(1); }
    strcpy(name, orifile); name[4] = 'U';
    p = strchr(name, '.'); strcpy(p, ".bin");
    /* Go to the subdirectory */
    name[3] = 0;
    if (mkdir(name, -1)) perror(name);
    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 */
    fprintf(u1file, header_text, reclen, name, Ofld, Nfld, Opm, Npm, 
	Omag, Nmag, Oxi, Nxi);
    fflush(u1file); bytes = ftell(u1file);
    /* Add blanks to comply with alignment requirements */
    p = blanks + (bytes&7) ;
    fwrite(p, 1, strlen(p), 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);
    fwrite(ochunks, 1, sizeof(ochunks), u1file);
    fflush(u1file);

    /* Change of Chunk depends on reclen. */
    RAmask = 0x7ff << 20 ;
    if (reclen == 18) RAmask = 0;

    /* Loop reading the input file, sorted in RA */
    id0 = 1; id = 0;
    while((i = fread(rec80, sizeof(long), 20, infile)) > 0) {
	swap(rec80, 20);
    	raz = rec80[0] & RAmask ;
    	if (raz != razp) {	/* A new zone appeared there */
	    nobj = store_chunk(razp>>20, id0, id, tmp_hand);
	    total += nobj ;
	    razp = raz ; id0 = id+1;
	    iacc = 0;
	}
	id++;
    	k = trec(rec80, 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++] = rec80[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>>20, 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); 

    /* Rewrite the List of Chunks (was just zeroes) */
    if (oochunks) {
        fseek(u1file, oochunks  , 0) ;
	fwrite(ochunks, 1, sizeof(ochunks), u1file);
	if (verbop) fprintf(stderr, "    (header rewritten)\n") ;
    }
    
    /* Edit some statistics */
    mean = bytes; mean /= total;
    printf("====> %10ld bytes, %8ldobj. %.1fb/o", bytes, total, mean);
    /* Compare photometries 6/7 bytes */
    mean = hlenph[2]; mean = 100.*mean/(hlenph[6]+hlenph[7]+hlenph[2]);
    printf(", ph2=%.1f%%", mean);
    /* Compute the compression factor */
    mean = hlenph[7]; mean = 100.*mean/(hlenph[6]+hlenph[7]);
    printf(", ph7=%.1f%%", mean);
    /* Compute the compression factor */
    mean = inbytes-bytes;
    mean = 100.*mean/inbytes;
    printf(" (red. %.0f%%)", mean);
    printf("\n");
    exit(0);
}