/* tyc2xe2.c -- process Hipparcos/Tycho-2 catalog into xe2 files */
/*
Copyright 2000, Steve VanDevender <stevev@hexadecimal.uoregon.edu>
Permission is granted to redistribute source code for this program if
and only if all of the following conditions are met:

1.  This copyright notice must not be removed or modified.
2.  Any changes you make to this source code must be documented.
3.  Binary versions must be distributed with complete source code.
4.  No money will be charged for distribution of binaries or source code.
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "readcat.h"
#include "packxe2.h"

/* data from Hipparcos catalog not available in Tycho-2 */
struct {
	double vmag;
	char sc[2];
} hipdata[121000];

unsigned char tyc2xe2[2539913][XE2SIZE]; /* urp */

/* helper function for qsort, compares xe2 record declination field */
static int sortxe2(unsigned char *e1, unsigned char *e2)
{
	if (e1[8] > e2[8]) return 1;
	if (e1[8] < e2[8]) return -1;
	if (e1[9] > e2[9]) return 1;
	if (e1[9] < e2[9]) return -1;
	if (e1[10] > e2[10]) return 1;
	if (e1[10] < e2[10]) return -1;
	return 0;
}

int main(int argc, char *argv[])
{
	FILE *hip, *tyc2;
	char line[1024];
	struct xe2data in;
	unsigned long i, n;

	hip = fopen("hip_main.dat", "r");
	if (hip == 0) {
		fprintf(stderr, "cannot open hip_main.dat\n");
		return 1;
	}
	printf("Reading Hipparcos catalog data from hip_main.dat ...\n");

	/* first, read in Hipparcos catalog data for things not
	   available in Tycho-2 */
	while (fgets(line, sizeof(line), hip) != 0) {
		char sc[13];

		get_integer(line, 9, 14, &i);
		if (i == BAD_INT || i > sizeof(hipdata)/sizeof(*hipdata)) continue;
		/* visual magnitude of star */
		get_fp(line, 42, 46, &hipdata[i].vmag);
		/* spectral class (must be of form [OBAFGKMWRNS][0-9 ]) */
		get_string(line, 436, 447, sc);
		if (strchr("OBAFGKMWRNS", sc[0]) != 0) {
			hipdata[i].sc[0] = sc[0];
		}
		else {
			hipdata[i].sc[0] = ' ';
		}
		if (strchr("0123456789 ", sc[1]) != 0) {
			hipdata[i].sc[1] = sc[1];
		}
		else {
			hipdata[i].sc[1] = ' ';
		}
	}
	fclose(hip);

	tyc2 = fopen("tyc2.dat", "r");
	if (tyc2 == 0) {
		fprintf(stderr, "can't open tyc2.dat\n");
		return 1;
	}
	printf("Reading Tycho-2 catalog data from tyc2.dat ...\n");

	n = 0;
	while (fgets(line, sizeof(line), tyc2) != 0) {
		char mult;

		in.typecode = STAR;
		get_char(line, 201, &mult);
		if (mult != ' ') {
			in.typecode = DOUBLE;
		}
		/* check for Hipparcos number in catalog entry */
		get_integer(line, 143, 148, &in.num1);
		/* if there is a Hipparcos number with a corresponding record 
		   in the previously-read data, give it a HIP namecode */
		if (in.num1 != BAD_INT && hipdata[in.num1].sc[0] != '\0' && hipdata[in.num1].vmag != BAD_DOUBLE) {
			in.namecode = HIP;
		}
		else {
			/* give it Tycho-2 namecode and catalog numbers */
			in.namecode = TYCHO2;
			get_integer(line, 1, 4, &in.num1);
			get_integer(line, 6, 10, &in.num2);
			get_integer(line, 12, 12, &in.num3);
		}
		get_fp(line, 16, 27, &in.ra);
		if (in.ra == BAD_DOUBLE) continue;
		in.ra /= 15.0; /* convert degrees to hours of RA */
		get_fp(line, 29, 40, &in.dec);
		if (in.dec == BAD_DOUBLE) continue;
		get_fp(line, 42, 48, &in.pmra);
		if (in.pmra == BAD_DOUBLE) in.pmra = 0.0;
		get_fp(line, 50, 56, &in.pmdec);
		if (in.pmdec == BAD_DOUBLE) in.pmdec = 0.0;
		if (in.namecode == HIP) {
			/* if this is a valid Hipparcos star, use the
			   Hipparcos spectral class and magnitude */
			in.sc[0] = hipdata[in.num1].sc[0];
			in.sc[1] = hipdata[in.num1].sc[1];
			in.mag = hipdata[in.num1].vmag;
		}
		else {
			double btmag, vtmag;

			/* Tycho-2 doesn't include spectral class */
			in.sc[0] = in.sc[1] = ' ';
			/* compute visual magnitude from BT and VT mags */
			get_fp(line, 111, 116, &btmag);
			get_fp(line, 124, 129, &vtmag);
			if (btmag != BAD_DOUBLE && vtmag != BAD_DOUBLE) {
				in.mag = vtmag - 0.09 * (btmag - vtmag);
			}
			else if (btmag != BAD_DOUBLE) {
				in.mag = btmag;
			}
			else if (vtmag != BAD_DOUBLE) {
				in.mag = vtmag;
			}
			else continue;
		}
		/* pack this into an xe2 structure */
		if (packxe2(&in, tyc2xe2[n])) {
			n++;
			if (n % 1000 == 0) printf("%ld\r", n);
		}
	}
	fclose(tyc2);
	printf("%ld Tycho-2 records.\n", n);

	printf("Sorting data ...\n");
	/* sort all records by declination value */
	qsort(tyc2xe2[0], n, XE2SIZE, sortxe2);

	printf("Writing output files hip.xe2 and tycho2.xe2.\n");
	hip = fopen("hip.xe2", "wb");
	tyc2 = fopen("tycho2.xe2", "wb");
	fprintf(hip, "XE2.\n");
	fprintf(tyc2, "XE2.\n");
	for (i = 0; i < n; i++) {
		if ((tyc2xe2[i][0] & 7) == HIP) {
			fwrite(tyc2xe2[i], XE2SIZE, 1, hip);
		}
		fwrite(tyc2xe2[i], XE2SIZE, 1, tyc2);
	}
	fclose(hip);
	fclose(tyc2);
	return 0;
}