/*++++++++++++++
.IDENTIFICATION pmm80.c
.LANGUAGE       C
.AUTHOR         Francois Ochsenbein [CDS]
.ENVIRONMENT    USNO-B1.0 Catalogue
.KEYWORDS       CDS Catalogue Server
.VERSION  1.0   06-Nov-2002
.VERSION  1.1   05-Dec-2002: Count the stargal outside the [0,11]  interval.
.COMMENTS       List one of the files containing the > 1,000,000,000 stars.
---------------*/

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

#define ITEMS(s)	sizeof(s)/sizeof(s[0])
#ifndef O_BINARY	/* Required for MS-DOS	*/
#define O_BINARY 0
#endif
#include <stdio.h>

static char verbop = 0;	/* Verbose Option	*/
static int swapping ;	/* Set to 1 or 2 (swap)	*/

static char usage[] = "\
Usage: pmm80 [-o binfile] input_file\n\
";

static char help[] = "\
      -o: write a binary copy\n\
   input: a \"bXXXX.cat\" original file\n\
";

/* Histogrammes */
static int hepo[1000];	/* Epoch Values	*/
static int hpos[1000];	/* Sigmas Position */
static int hpm[10000];	/* Proper Motions */
static int hspm[1000];	/* Sigmas */

static int xmag[32];	/* Existence */
static int hgal[20];	/* star/gal */
static int hgal15[6];	/* star/gal>11 */
static char *hgal15tit[6] = { 
   "blu#1", "red#1", "blu#2", "red#2", "ir(N)"
};
static int hclass[10000];	/* Star/gal class function of Survey */

static int hmag[10000];	/* Magnitudes */
static int hfld[1000];	/* Field  */
static int hxie[10000];	/* xi-eta */

/*==================================================================
		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 ;
    }
}

/*==================================================================
		Histogram Summary
 *==================================================================*/
static void show_histo(char *title, int *freq, int size)
/*++++++++++++++++
.PURPOSE  Histogram Summary
.RETURNS  ---
.REMARKS  
-----------------*/
{
  int *p, *e, *bounds[2], *bound1[2];
  long val1[2], tot = 0;
  int bins = 0;
    e = freq + size;
    /* Compute the Total */
    for(p=freq; p<e; p++) tot += *p;

    /* Filled bins */
    for(p=freq; p<e; p++) { if(*p) bins++; }

    fprintf(stderr, "%s %5dbins%3d%%filled.", title, bins, 100*bins/size);

    /* Find the Ranges */
    for (p=freq; (p<e) && (*p == 0); p++) ;    bounds[0] = p;
    for (--e; (e>=freq) && (*e == 0); e--) ;   bounds[1] = e;
    fprintf(stderr, " Range %5d'%-5d", bounds[0]-freq, bounds[1]-bounds[0]);

    /* Find the 1%-99% Range */
    bound1[0] = bounds[0]; val1[0] = tot/100;		/* 1% */
    bound1[1] = bounds[1]; val1[1] = tot - val1[0];	/* 99% */
    for (p=bounds[0], tot=0, bins=0; p<e; p++) {
	tot += *p;
	if (tot<val1[0]) { if (*p) bins++; bound1[0] = p; }
	if (tot>val1[1]) break;
	bound1[1] = p;
    }
    while (p<e) { if (*p) bins++; p++; }
    fprintf(stderr, " [-2%%:%4dbins %5d'%-5d]", 
	bins, bound1[0]-freq, bound1[1]-bound1[0]);

    fprintf(stderr, "\n");
}


/*==================================================================
		Convert the original record
 *==================================================================*/

/* Convert the 80-byte into the Structure */
static int tyc_set (long *rec, USNOBtyc *prec)
/*++++++++++++++++
.PURPOSE  Interpret the TYCHO-2
.RETURNS  0 (OK)
.REMARKS  The magnitudes are already set.
-----------------*/
{
  long val;

    prec->flags = USNOB_TYC;

    /* Zero the errors */
    prec->e_ra = prec->e_sd = 0;
    prec->epoch = 20000;
    prec->fit_ra = prec->fit_sd = 0;
    prec->e_pmra = prec->e_pmsd = 0;

    /* Tycho Number */
    val = rec[3];
    if      (val == 21) prec->flags |= USNOB_TS1;
    else if (val == 22) prec->flags |= USNOB_TS2;
    prec->fileno = rec[3];
    prec->recno  = rec[4];

    val = rec[10];
    prec->TYC1 = val;

    val = rec[11];
    prec->TYC3 = val%10; val /= 10;
    prec->TYC2 = val;

    return(0);
}

int usnob1_ori(long *rec, USNOBrec *prec)
/*++++++++++++++++
.PURPOSE  Interpret the original USNO catalog 
.RETURNS  0 (OK) 
.REMARKS  We assume that "zone" and "id" are already filled.
   	Histograms computed.
-----------------*/
{
  long word;
  int i, m, x, val;
    
    /* Install RA, SD ... */
    prec->flags = 0;
    prec->ndet  = 0;
    prec->ra = rec[0]*10; 
    prec->sd = rec[1]*10; 

    /* Proper Motions = iPSSSSAAAA */
    word = rec[2];
    prec->pmra = (val=word%10000)*2 - 10000; word /= 10000; hpm[val] += 1;
    prec->pmsd = (val=word%10000)*2 - 10000; word /= 10000; hpm[val] += 1;
    prec->muprob = word%10; word /= 10;
    prec->spare = 0;
    if (word) prec->flags |= USNOB_PM;

    /* Magnitudes are common with Tycho records */
    m = 1; x = 0;
    for (i=0; i<5; i++) {
	if (rec[5+i]) x |= m; 
	prec->mag[i] = val = rec[5+i]%10000;
	if (val) hmag[val] += 1;
	m <<= 1;
    }
    xmag[x] += 1;

    /* Zero the offsets in position */
    prec->theta = prec->xy[0] = prec->xy[1] = 0;
    prec->rho = -1;

    /* Zero the photometries */
    memset(&prec->phot, 0, sizeof(prec->phot));

    /* Test for a Tycho-2 entry */
    word = rec[3];
    if (word < 100000000) 
	return(tyc_set(rec, (USNOBtyc *)prec));
    prec->e_pmra = (val=word%1000); hspm[val] += 1; word /= 1000;
    prec->e_pmsd = (val=word%1000); hspm[val] += 1; word /= 1000;
    prec->fit_ra = (val=word%10);   word /= 10;
    prec->fit_sd = (val=word%10);   word /= 10;
    prec->ndet   = (val=word%10);   word /= 10;
    if (word) prec->flags |= USNOB_SPK;

    /* Sigmas */
    word = rec[4];
    prec->e_ra = (val=word%1000); hpos[val] += 1; word /= 1000;
    prec->e_sd = (val=word%1000); hpos[val] += 1; word /= 1000;
    prec->epoch = (val=word%1000) + 19500; hepo[val] += 1; word /= 1000;
    if (word) prec->flags |= USNOB_YS4;

    /* Residuals, etc */
    /* We Add in hgal15 which photometries have value>11 */
    for (i=0; i<5; i++) {
	if (rec[5+i] == 0) {	/* EMPTY Magnitude !!!	*/
	    if (rec[10+i]) fprintf(stderr, "++++Wmag=%d, Wphot=%d\n",
		rec[5+i], rec[10+i]) ;
	    continue;
	}
	word = rec[5+i]/10000;	/* Magnitudes already used */
	x = word%10000;		/* Field + Survey */
	prec->phot[i].field   = val=word%1000; word /= 1000; hfld[val] += 1;
	prec->phot[i].survey  = val=word%10; word /= 10;
	prec->phot[i].stargal = val=word; hgal[val] += 1; 
	if ((val>11) && (val<19)) { hgal15[i] += 1; hclass[x] += 1; }
	word = rec[10+i];
	prec->phot[i].xi  = (val=word%10000)-5000; word /= 10000; hxie[val]+=1;
	prec->phot[i].eta = (val=word%10000)-5000; word /= 10000; hxie[val]+=1;
	prec->phot[i].calib = word;
    }

#if 0
    /* Lookback indexes */
    for (i=0; i<5; i++) 
	prec->index[i] = rec[15+i];
#endif

    return(0);
}


/*==================================================================
		Main Program
 *==================================================================*/
main (argc, argv) int argc; char **argv;
{
  long in[20], counts[10];
  USNOBrec rec;
  FILE *f, *of ;
  char *p, *oname ;
  int i;

    of = (FILE *)0;
    oname = (char *)0;
    while (argc > 2) {
	p = *++argv; --argc;
	if (*p != '-') break; 
	if (p[1] != 'o') break;
	p = *++argv; argc--;
	of = fopen(p, "w");
	if (!of) perror(p);
	oname = p;
    }
    if (argc != 2) {
	fprintf(stderr, "%s%s", usage, help) ;
	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;		/* ASSUME THIS TYPE OF SWAPPING */
    fprintf(stderr, "#...swapping type=%d\n", swapping) ;

    memset(&rec, 0, sizeof(rec)) ;
    while (--argc > 0) {
	fprintf(stderr, "----Dealing with file: %s", *++argv) ;
        if (!(f = fopen(*argv, "rb"))) {
	    perror(*argv) ;
	    exit(1) ;
	}
	fflush(stderr) ;

	/* Get the Zone Number (from file xxxx.cat) */
	p = *argv; p += strlen(p) -8 ;
	rec.zone = atoi(p);

	/* Zero all histogram values */
	memset(counts, 0, sizeof(counts));
	memset(hepo, 0, sizeof(hepo));  /* Epoch Values */
	memset(hpos, 0, sizeof(hpos));  /* Sigmas Position */
	memset(hpm, 0, sizeof(hpm));    /* Proper Motions */
	memset(hspm, 0, sizeof(hspm));  /* Sigmas */
	memset(xmag, 0, sizeof(xmag));  /* Existence */
	memset(hgal, 0, sizeof(hgal));  /* star/gal */
	memset(hgal15, 0, sizeof(hgal15));  /* star/gal */
	memset(hmag, 0, sizeof(hmag));  /* Magnitudes */
	memset(hfld, 0, sizeof(hfld));  /* Field  */
	memset(hxie, 0, sizeof(hxie));  /* xi-eta */


        for (rec.id=0; fread(in, sizeof(long), 20, f) > 0 ; ) {
	    rec.id++ ;
	    swap(in, 20);
	    usnob1_ori(in, &rec);
	    if(of) fwrite(&rec, sizeof(rec), 1, of);
	    counts[rec.ndet] += 1;
	}
	fprintf(stderr, " => %d objects.\n", rec.id);
	close(f) ;
	/* Existence of Magnitudes */
	fprintf(stderr, "   ndet:");
	for (i=0; i<=5; i++) fprintf(stderr, " %d=%-7d", i, counts[i]);
	fprintf(stderr, "\n");

	show_histo(" ErrPos:", hpos, ITEMS(hpos));
	show_histo("  Epoch:", hepo, ITEMS(hepo));
	show_histo("     pm:", hpm, ITEMS(hpm));
	show_histo("    spm:", hspm, ITEMS(hspm));
	show_histo("  Field:", hmag, ITEMS(hfld));
	show_histo("    mag:", hmag, ITEMS(hmag));
	show_histo(" xi-eta:", hxie, ITEMS(hxie));
	show_histo("stargal:", hgal, ITEMS(hgal));  /* star/gal */
	for (i=0; i<ITEMS(hgal); i++) {
	    if (hgal[i]) fprintf(stderr, "%5d=%-10d", i, hgal[i]);
	}
	counts[0] = 0; for (i=0; i<ITEMS(hgal15); i++) counts[0] += hgal15[i];
	if (counts[0]) {
	    fprintf(stderr, "\n%-16s", "---stargal>11:") ;
	    for (i=0; i<ITEMS(hgal15); i++) {
	        if (hgal15[i]) fprintf(stderr, "%s=%-10d", 
	  	  hgal15tit[i], hgal15[i]);
	    }
	    fprintf(stderr, "\n%-16s", "---stargal>11:") ;
	    for (i=0; i<ITEMS(hclass); i++) {
	        if (hclass[i]) fprintf(stderr, "F%04d=%-10d", i, hclass[i]);
	    }
	}
	fprintf(stderr, "\n");
	if(oname) fprintf(stderr, "====Binary Copy to: %s\n", oname);
	fprintf(stderr, "\n");
    }

#if 0
    fprintf(stderr, "==== Total Histo Bmag, Rmag:\n") ;
    for (i=0; i<1000; i++)  {
	if ((b[i] == 0) && (r[i] == 0)) continue ;
	fprintf(stderr, "%8d %10d %10d\n", i, bbb[i], rrr[i]) ;
    }
#endif
}