/*++++++++++++++
.IDENTIFICATION pmmake.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-A2.0 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   01-Mar-1997
.VERSION  2.0   21-Oct-1988: Version USNO-A2.0
.COMMENTS       Create the hieracrhy of compressed binary files
	containing the 525,000,000 stars.

	The original files are split into pieces of 2.5*2.5 degrees
	named  [ns]DDMM/[nsm]HHMM.pmm
	
	Data for 1 star are compressed in 7 bytes, containing
	field index (3bits), Gflag (1b), Q flag (1b),
	oRA10mas (11bits), oSD10mas (22b), mB (9bits), mR (9bits).
	
	The header consists in:
	1) an ascii line 
	   USNO-A2.0(7) ..  or USNO-A2.0(8)  	terminated by \n
	   made of a multiple of 4 bytes
   	2) ntab, otab, zone, id_min, RA_min, SPD_min, id_max, RA_max, SPD_max, 
	    i4    i2    i2     i4     i4       i4      i4       i4     i4
   	    8*field     8*epoch
   	  = 80 bytes; 312*I2 special magnitudes
	    => ntab is the number of elements in the offsets table
	       (normally 3*512+1)
	    => otab is the starting seek byte of index
	    => zone is the number of file on CD-ROM (0000, 0075,... 1725)
	    => id0  is number of first star 
        3) An index table of ntab offsets providing the starting seek byte 
           (from the beginning of the file) of a subregion
   	4) The records of 7* or 8*(number of stars).
	   The original values can be retrieved as, if
	       i = subregion index (from 0 to ntab-2)
	       o = byte offset from beginning of file (seek position)
	   RA = RA_min + (i<<11) + oRA10mas  
	     or RA_min + (i<<18) + oRA10mas
	   SD = SD_min + oSD10mas
	   ID = USNO + (o-seek0)/7
	5) Magnitudes are stored with an offset of 50 (i.e. 20th magnitude
	   star has a stored vaLue of 150). Values outside the interval
	   [0,199] (i.e. magnitudes 5 to 24.9) are stored as an index
	   in the "special magnitude" table.
---------------*/

#include <pmm_def.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 char *orifile;		/* Input file name	    */
static HEADER headers[3] ;	/* Headers of current files */
static HEADER8 *header8 = (HEADER8 *)headers ;
static FILE *fpmm[3] ;		/* File handles of output   */
static char fnames[3][20] ;	/* Name of pmm files	    */
static int  SDzone ;		/* Current Declination zone */
static int  reclen = 7;		/* Current type (7 or 8)    */
static int  mplates ; 		/* 8 or 16 */
static long oheader;		/* Seek position of header  */
static char *ep_red, *ep_blue ;
static long total_size ;
static long ostart ;		/* Starting byte in input   */

static char verbop = 0;	/* Verbose Option	*/
static long histo[1000] ;
static long misto[1000] ;
static FILE *fhisto ;

static char intro[] = "\
An index is made of the 3 parts:\n\
   1) One ascii line starting by \"USNO-A2.0(7)\" or \"USNO-A2.0(8)\"\n\
   2) A header containing the limits RA, SPD in units of 10mas,\n\
      and the list of the 8 (or 16) plates (and plate epochs) used there\n\
   3) A table of 312 short integers with \"special\" (>=200) magnitudes\n\
   4) A table of seek positions for the 1319(11)+1 \"blocks\" \n\
      each containing a zone of 2**11(2**18) * 10mas\n\
   5) The data for each star, as 7-(8-)byte records\n\
      Gflag (1b) Qflag(1b) Plate(3(4)b) dRA (11(18)b)\n\
      dSD(22b) Bmag(9b) Rmag(9b)\n\
" ;

static char usage[] = "\
Usage: pmmake [-v] [-r root_name] [-om mag_histo_file] input_zone\n\
";

static char help[] = "\
      -v: verbose option\n\
 -r name: defines the name of Root (defaut from current directory)\n\
 -om fil: write out an histogram of imput magnitudes\n\
   input: a \"zoneXXXX.cat\" original file\n\
";

static long index_id[145] ;

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

static char *load_plates(char *name) 
/*++++++++++++++++
.PURPOSE  Load the list of plates (file "plates.yr" with plate number epoch)
.RETURNS  ---
-----------------*/
{
  long bytes;
  char *plates_yr ;
  int file;

    file = open(name, O_BINARY);
    if (file <= 0) {  perror(name);  return ((char *)0); }
    bytes = lseek(file, 0L, 2); lseek(file, 0L, 0);
    if (plates_yr = malloc(bytes+1)) {
    	if (read(file, plates_yr, bytes) != bytes) {
    	    perror(name);
    	    free(plates_yr);
	    plates_yr = (char *)0 ;
    	    return(plates_yr);
    	}
	plates_yr[bytes] = 0 ;
    }
    else perror(name) ;
    close(file);
    return(plates_yr);
}

static int code_mag(int mag, short *amag)
/*++++++++++++++++
.PURPOSE  Reserve a cell for a magnitude outside standard 0..200 range
.RETURNS  [0,199]=standard ; >=200=value saved in header.
-----------------*/
{
  short *s ; int m ;
    m = mag-50 ;
    if ((m >= 0) && (m < 200)) return(m) ;
    if (*amag == 0) amag[0] = -50, amag[1] = 949 ;
    for (s=amag; (*s != m) && *s ; s++ ) ;
    if (!*s) *s = m ;
    m = 200 + (s-amag) ;
    return(m) ;
}

/*==================================================================
		Convert the Input Record(s)
 *==================================================================*/
static int trec(HEADER *h, long in[3], unsigned char *out)
/*++++++++++++++++
.PURPOSE  Convert the input (RA+DE+MAG) into 7- or 8-byte record
.RETURNS  -1=Error / 0 .. INDSIZ-1: block number
-----------------*/
{
  int mb, mr, pl, f1 ;	/* f1 = flags ACT / Q */
  long mg, ra, sd ;
  short *plate, *mags ;
  long *epoch ;
  double yb, yr, atof() ;

    /* Positions */
    if ((in[0] < h->ra0) || (in[0] > h->ra1)) {
    	fprintf(stderr, "****Incompatible RA: %ld(%.6f) [%ld %ld]\n", 
    	    in[0], in[0]/3.6e5, h->ra0, h->ra1) ;
        return(-1) ;
    }
    if ((in[1] < h->sd0) || (in[1] > h->sd1)) {
    	fprintf(stderr, "****Incompatible SD: %d(%+.6f) [%ld %ld]\n",
    	    in[1], in[1]/3.6e5-90. -90., h->sd0, h->sd1) ;
        return(-1) ;
    }
    pl = f1 = 0 ;
    ra = in[0] - h->ra0 ;
    sd = in[1] - h->sd0 ;
    mg = in[2] ;
    if (mg < 0) f1 |= PMM_ACT, mg = -mg ;
    if (mg >= 1000000000) f1 |= PMM_m, mg -= 1000000000 ;

    /* Plate number */
    pl = mg/1000000 ; mg -= (pl*1000000) ;	/* Plate number */
    if ((pl <= 606) && (in[1] <= 75*3600*100))  /* Southern plate */
        pl = -pl ;	
    if ((pl == 0) && (!f1)) fprintf(stderr,	/* No ACT/Plate */
	"++++FFF is zero for RA=%09ld SD=%08ld MG=%11ld\n",
    	in[0], in[1], in[2] ) ;
    plate = h->plate ;
    if (pl != 0) for (++plate; *plate && (*plate != pl); plate++) ;
    if (pl && (*plate==0))  { 			/* New Plate 	*/
      char aplate[6], *p; 	
    	*plate = pl ; 
	sprintf(aplate, "%4d", pl);
	for (p = ep_blue; p; p = strchr(p, '\n')) {
	    if (*p == '\n') p++ ;
	    if (strncmp(p, aplate, 4) == 0) break ;
	}
	if (p)  yb = atof(p+4) ;	/* Blue Epoch */
	else    yb = 0, fprintf(stderr,
	    "++++Unknown blue plate %d for RA=%09ld SD=%08ld MG=%11ld\n",
	     pl, in[0], in[1], in[2]) ;
	for (p = ep_red; p; p = strchr(p, '\n')) {
	    if (*p == '\n') p++ ;
	    if (strncmp(p, aplate, 4) == 0) break ;
	}
	if (p)  yr = atof(p+4) ;	/* Blue Epoch */
	else    yr = 0, fprintf(stderr,
	    "++++Unknown red  plate %d for RA=%09ld SD=%08ld MG=%11ld\n",
	     pl, in[0], in[1], in[2]) ;
	if (yb && yr) yr = (yb+yr)/2. ;
	else if (yr == 0) yr = yb ;
	if (yr) {
	    epoch = (long *)(&(h->plate[mplates])) ;
	    epoch[plate - h->plate] = 1000.*yr ;	/* Take 3 decimals */
	}
    }
    pl = plate - h->plate ;			/* Plate index 	*/
    if (pl >= mplates) {
	if (verbop) fprintf(stderr, " (abort)\n") ;
    	fprintf(stderr, 
	    "++++Too many (%d) plates for RA=%09ld SD=%08ld MG=%11ld\n",
    	    pl, in[0], in[1], in[2]) ;
    	for (plate = h->plate; --pl>=0; plate++) fprintf(stderr, "%8d", *plate);
    	fprintf(stderr, "\n") ;
    	pl = plate - h->plate ;
	return(-1) ;
    }

    /* Blue magnitude */
    mb = mg/1000 ; mg -= (mb*1000) ;
    if ((mb == 0) && (*plate == 0)) mb = 999 ;
    histo[mb] += 1 ;

    /* Red magnitude */
    mr = mg ;
    histo[mr] += 1 ;

    /* Compact the result */
    if (reclen == 8) {		/* 18 bits for RA */
	mags = ((HEADER8 *)h)->mags ;
	out[0] = f1 | (pl<<2) | ((ra&0x030000)>>16) ;
	out[1] = (ra&0x00ff00) >> 8 ;
	out[2] = ra & 0xff ;
	out++ ;
	ra >>= 18 ;
    }
    else {			/* 11 bits for RA */
	mags = h->mags ;
	out[0] = f1 | (pl<<3) | ((ra&0x700)>>8) ;
	out[1] = ra & 0xff ;
	ra >>= 11 ;
    }

    /* Code the Magnitudes */
    mb = code_mag(mb, mags);
    if (mb > 511) {
	if (verbop) fprintf(stderr, "\n") ;
	fprintf(stderr, 
    	"++++Bad mB flux (%d) for RA=%09ld SD=%08ld MG=%11ld\n",
    	mb, in[0], in[1], in[2]) ;
	mb = 201 ;	/* --> 999 */
    	if (verbop) fprintf(stderr, "....filling '%s': %8s", fnames[0], ""), 
		fflush(stderr) ;
    }
    mr = code_mag(mr, mags);
    if (mr > 511) {
	if (verbop) fprintf(stderr, "\n") ;
	fprintf(stderr, 
    	"++++Bad mR flux (%d) for RA=%09ld SD=%08ld MG=%11ld\n",
    	mr, in[0], in[1], in[2]) ;
	mr = 201;	/* --> 999 */
    	if (verbop) fprintf(stderr, "....filling '%s': %8s", fnames[0], ""), 
		fflush(stderr) ;
    }

    sd <<= 2 ;
    mb <<= 1 ;
    out[2] = (sd&0xff0000)>>16 ;
    out[3] = (sd&0x00ff00)>>8 ;
    out[4] = (sd&0xff) | ((mb&0x300)>>8) ;
    out[5] = (mb&0xff) | ((mr&0x100)>>8) ;
    out[6] = mr & 0xff ;

    return(ra) ;
}

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

static int open3(int RAzone, int SDzon, int rectyp /* 7 or 8 */)
/*++++++++++++++++
.PURPOSE  Create the file for the specified RA zone
.RETURNS  0 = OK, -1 = Error
.REMARKS  RAzone in range [0,144[ (360/2.5), SDzone in range [0,72[
	headers and fpmm statis are set.
-----------------*/
{
  static char s[3] = {'s', 'm', 'n'} ;
  static char blanks[] = "       \n" ;
  long RAmin, SDmin ;
  int i, k ;
  char *p ;
  long bytes ;
  
    reclen = rectyp ;
    SDzone = SDzon ;
    if (reclen == 7) mplates = 8;
    else	    mplates = 16;

    RAmin = RAzone * XIZE ;
    SDmin = SDzone * XIZE ;
    if ((RAzone < 0) || (RAzone >= 144) || (SDzone < 0) || (SDzone >= 72)) {
    	fprintf(stderr, "****open3: bad zone RAmin=%ld, SDmin=%ld\n",
    	    RAmin, SDmin) ;
    	return(-1) ;
    }

    memset(histo, 0, sizeof(histo)) ;
    k = 0 ;
    	sprintf(fnames[k], "%02d%d0.pmm", RAzone/6, RAzone%6) ;
    	fpmm[k] = fopen(fnames[k], "wb") ;
    	if (! fpmm[k]) {
    	    perror(fnames[k]) ;
    	    return(-1) ;
    	}
    	if (verbop) fprintf(stderr, "....filling '%s': ", fnames[k]) ;

    	/* Write the header text, must be a multiple of 8 bytes */
    	bytes = sizeof(header_text)-1 ;
    	fwrite(reclen == 8 ? header8text : header_text, 1, bytes, fpmm[k]) ;
    	i = strlen(orifile) ; fwrite(orifile, 1, i, fpmm[k]) ; bytes += i ;
    	p = blanks + (bytes&7) ;
    	i = strlen(p) ;
    	fwrite(p, 1, i, fpmm[k]) ; bytes += i ;
	oheader = bytes;			/* Keep the header position */

	/* Write the header */
    	memset(&headers[k], 0, sizeof(HEADER)) ;
    	headers[k].reclen = reclen ;
    	headers[k].zone = 25*SDzone ;

	ostart  = ftell(infile) ;
	headers[k].id0 = 1 + ostart/12 ;
    	headers[k].ra0 = headers[k].ra1 = RAmin ;
    	headers[k].sd0 = headers[k].sd1 = SDmin + k*XIZE ;
    	headers[k].ra1 += 3*XIZE-1;
    	headers[k].sd1 += 3*XIZE;

    	if (reclen==8) {
	    mplates = 16 ;
	    headers[k].ntab = INDSIZ8 ,
	    headers[k].otab = bytes + sizeof(HEADER8) - INDSIZ8*sizeof(long) ;
	    fwrite(&headers[k], sizeof(HEADER8), 1, fpmm[k]) ;
	    bytes += sizeof(HEADER8) ;
	}
	else {
	    mplates = 8 ;
	    headers[k].ntab = INDSIZ ,
	    headers[k].otab = bytes + sizeof(HEADER) - INDSIZ*sizeof(long) ;
	    fwrite(&headers[k], sizeof(HEADER), 1, fpmm[k]) ;
	    bytes += sizeof(HEADER) ;
	}

	index_id[RAzone] = headers[k].id0 ;

    	if (verbop) fprintf(stderr, "%8ld", bytes), fflush(stderr),
	    fflush(fpmm[k]);

    return(0) ;
}

static long close3(FILE *ofile, HEADER *h) 
/*++++++++++++++++
.PURPOSE  Arrange the file (move from 8-byte records to 7-byte records)
.RETURNS  Number of Objects in the file
.REMARKS  Output the offset position
-----------------*/
{
  long bytes, stars, *ob, *ob0, *ob1 ;
  int file, i, len, reclen ;

    fflush(ofile) ;
    file = fileno(ofile) ;
    bytes = lseek(file, 0L, 2) ;
    reclen = h->reclen ;
    
    /* Update statistics */
    for (i=0; i<1000; i++) misto[i] += histo[i] ;
    if (fhisto) {
	for (i=0; i<1000; i++)  if (histo[i])
	    fprintf(fhisto, "%3d %3d %8d\n", h->ra0/XIZE, i, histo[i]);
	fprintf(fhisto, "\n") ;
    }

    /* Update the Header Text */
    lseek(file, oheader, 0) ;
    if (verbop) fprintf(stderr, "\b\b\b\b\b\b\b\b%8db (header rewriting", 
    	bytes), fflush(stderr) ;

    if (reclen == 8) { HEADER8 *h8 ;
	h8 = (HEADER8 *) h;
	len = sizeof(HEADER8) ;
	ob0 = h8->oblock ;
	ob1 = ob0 + INDSIZ8-1 ;
    }
    else {
	len = sizeof(HEADER) ;
	ob0 = h->oblock ;
	ob1 = ob0 + INDSIZ-1 ;
    }

    stars = (bytes - *ob0) / reclen;	/* Number of stars 	*/
    h->id1 = h->id0 + (stars - 1) ;
    *ob1 = bytes ; 		/* Upmost value (= file size)	*/
    /* Fill out the offsets: some regions may be empty ...	*/
    for ( --ob1; ob1 > ob0; ob1--) 
	if (*ob1 == 0) ob1[0] = ob1[1] ;

    i = write(file, h, len) ;
    if (i != len) { perror("\n****Rewriting Header") ; exit(1) ; }
    else if (verbop) fprintf(stderr, "\b\b\bten)") ;
    close(file) ;
    total_size += bytes ;
    
    if (verbop) fprintf(stderr, " %7ld obj. (%d)\n", stars,
	h->reclen) ;
    
    return(stars) ;
}

/*==================================================================
		Main Program
 *==================================================================*/
main (argc, argv) int argc; char **argv;
{
  int i, k, zone, RAzone;
  char *b, *p, *root;
  char name[20] ;
  long rec12[3], total, nobj, *oblock ;
  unsigned char recb[8] ;

    /* Load the Epochs */
    ep_red  = load_plates("red.yr") ;
    ep_blue = load_plates("blue.yr") ;

    /* Check the Arguments */
    total = 0 ;
    root = "." ;
    while (-- argc > 0) {
	p = *++argv;
	if (*p == '-') switch (*++p) {
	  case 'h':
	    printf("%s%s", usage, help) ;
	    exit(0);
	  case 'r':
	    root = *++argv; -- argc;
    	    if (chdir(root)) { perror(root); exit(1); }
	    continue;
	  case 'o':
	    p = *++argv; -- argc;
	    fhisto = fopen(p, "w") ;
	    if (!fhisto) perror(p);
	    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) ;
    }

    /* Check Input File, at poles we've a different structure */
    p = *argv;			/* File Name 	*/
    p += strlen(p) - 8;
    zone = atoi(p) ;		/* Zone in 0.1deg	*/
    SDzone = zone/25 ;		/* Zone in range 0..71 	*/
    orifile = p - 4 ;		/* Name original file	*/
    				/* Type as 7 or 8 bytes	*/
    infile = fopen(*argv, "rb");
    if (!infile) { perror(*argv); exit(1); }

    /* Define Output Names */
    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) ;
    
    if (chdir(name)) { perror(name); exit(1); }
    if (fhisto) 
	fprintf(fhisto, "Zone: %s\nRAz mag      N\n\n", name);

    /* Open the 3 files dealing with this 7.5 degrees zone */
    RAzone = 0 ;
    reclen = 7, oblock = headers[0].oblock ;
    open3(RAzone, SDzone, 7) ;

    /* Loop reading the input file, sorted in RA */
    while((i = fread(rec12, sizeof(long), 3, infile)) > 0) {
    	zone = rec12[0] / (3*XIZE) ;
    	if (zone != RAzone) {	/* A new zone appeared there */
	    nobj = close3(fpmm[0], headers) ;
	    printf("----File '%s'(%d): %8ld obj.\n", fnames[0], reclen, nobj) ;
	    fflush(stdout) ;
	    total += nobj ;
	    fseek(infile, -12L, 1) ;
	    RAzone = zone ;
    	    reclen = 7, oblock = headers[0].oblock ;
	    open3(3*RAzone, SDzone, 7) ;
	    continue;
	}
    	k = trec(headers, rec12, recb) ;
	if (k < 0) {		/* Doesn't wirk with 7 bytes/record */
	    if (reclen == 8) {
		fprintf(stderr, "****Terminal Error****") ;
		exit(1) ;
	    }
	    fclose(fpmm[0]) ;
	    fseek(infile, ostart, 0) ;
	    reclen = 8, oblock = header8->oblock ;
	    open3(3*RAzone, SDzone, 8) ;
	    continue ;
	}
	if (! oblock[k]) {
	    oblock[k] = ftell(fpmm[0]) ;
	    if (verbop) fprintf(stderr, "\b\b\b\b\b\b\b\b%8ld", 
	    	oblock[k]), fflush(stderr) ;
    	}
    	fwrite(recb, reclen, 1, fpmm[0]) ;
    }
    if (i < 0) perror(*argv) ;

    /* Terminate last file */
    	fclose(infile) ;
	nobj = close3(fpmm[0], headers) ;
	printf("----File '%s'(%d): %8ld obj.\n", fnames[0], reclen, nobj) ;
	fflush(stdout) ;
    	total += nobj ;
	index_id[144] = headers[0].id1 + 1 ;
	for (i=0; i<145; i++) {
	    if (index_id[i] != 0) continue ;
	    for (k=i; index_id[k] == 0; k++) ;
	    index_id[i] = index_id[k] ;
	}
    
    /* Write out the index.145 file */
    i = open("index.145", 1|O_CREAT|O_BINARY, 0644) ;
    if (i < 0) perror("index.145");
    else write(i, index_id, sizeof(index_id)), close(i) ;

    if (fhisto) 
	for (i=0; i<1000; i++) if (misto[i])
	    fprintf(fhisto, "=== %3d %8d\n", i, misto[i]);

    fprintf(stderr, "====Total: %d files, %ld obj. (%.3f Mbytes)\n", 
	24, total, total_size/1024./1000.) ;
    exit(0);
}