/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/

#include <arpa/inet.h>
#include <netinet/in.h>
#include <assert.h>
#include <stdio.h>
#include <stdint.h>

#include "usnob.h"
#include "an-endian.h"

int unsob_get_survey_epoch(int survey, int obsnum) {
	switch (survey) {
	case USNOB_SURVEY_POSS_I_O:
	case USNOB_SURVEY_POSS_I_E:
		return 1;

	case USNOB_SURVEY_POSS_II_J:
	case USNOB_SURVEY_POSS_II_F:
	case USNOB_SURVEY_POSS_II_N:
	case USNOB_SURVEY_SERC_J:
		//case USNOB_SURVEY_SERC_EJ:
	case USNOB_SURVEY_AAO_R:
	case USNOB_SURVEY_SERC_I:
	case USNOB_SURVEY_SERC_I_OR_POSS_II_N:
		return 2;

	case USNOB_SURVEY_ESO_R:
		//case USNOB_SURVEY_SERC_ER:
		{
			if (obsnum == 1)
				return 1;
			if (obsnum == 3)
				return 2;
			return -1;
		}
	default:
		return -1;
	}
}

unsigned char usnob_get_survey_band(int survey) {
	// from Tables 1 and 3 (esp footnote h) of the USNO-B paper.
	switch (survey) {
	case USNOB_SURVEY_POSS_I_O:
		return 'O'; // blue (350-500 mm)
	case USNOB_SURVEY_POSS_I_E:
		return 'E'; // red  (620-670 mm)
	case USNOB_SURVEY_POSS_II_J:
	case USNOB_SURVEY_SERC_J:
		//case USNOB_SURVEY_SERC_EJ:
		return 'J'; // blue (385-540 mm)
	case USNOB_SURVEY_POSS_II_F:
	case USNOB_SURVEY_ESO_R:
		//case USNOB_SURVEY_SERC_ER:
	case USNOB_SURVEY_AAO_R:
		return 'F'; // red  (590, 610 or 630-690 mm)
	case USNOB_SURVEY_POSS_II_N:
	case USNOB_SURVEY_SERC_I:
	case USNOB_SURVEY_SERC_I_OR_POSS_II_N:
		return 'N'; // infrared (715 or 730-900 mm)
	default:
		return '\0';
	}
}

int usnob_get_blue_mag(usnob_entry* entry, float* mag) {
    float sum = 0.0;
    int n = 0;
    // Dumbass mag averaging.
    if (usnob_is_observation_valid(entry->obs + OBS_BLUE1)) {
        sum += entry->obs[OBS_BLUE1].mag;
        n++;
    }
    if (usnob_is_observation_valid(entry->obs + OBS_BLUE2)) {
        sum += entry->obs[OBS_BLUE2].mag;
        n++;
    }
    if (n == 0)
        return -1;
    *mag = sum / (double)n;
    return 0;
}
int usnob_get_red_mag(usnob_entry* entry, float* mag) {
    float sum = 0.0;
    int n = 0;
    // Dumbass mag averaging.
    if (usnob_is_observation_valid(entry->obs + OBS_RED1)) {
        sum += entry->obs[OBS_RED1].mag;
        n++;
    }
    if (usnob_is_observation_valid(entry->obs + OBS_RED2)) {
        sum += entry->obs[OBS_RED2].mag;
        n++;
    }
    if (n == 0)
        return -1;
    *mag = sum / (double)n;
    return 0;
}
int usnob_get_infrared_mag(usnob_entry* entry, float* mag) {
    if (usnob_is_observation_valid(entry->obs + OBS_N)) {
        *mag = entry->obs[OBS_N].mag;
        return 0;
    }
    return -1;
}


anbool usnob_is_usnob_star(usnob_entry* entry) {
	return (entry->ndetections >= 2);
}

anbool usnob_is_observation_valid(struct observation* obs) {
	return (obs->field > 0);
}

anbool usnob_is_band_blue(unsigned char band) {
	return (band == 'O' || band == 'J');
}

anbool usnob_is_band_red(unsigned char band) {
	return (band == 'E' || band == 'F');
}

anbool usnob_is_band_ir(unsigned char band) {
	return (band == 'N');
}

anbool usnob_is_observation_blue(struct observation* obs) {
	return usnob_is_band_blue(usnob_get_survey_band(obs->survey));
}

anbool usnob_is_observation_red(struct observation* obs) {
	return usnob_is_band_red(usnob_get_survey_band(obs->survey));
}

anbool usnob_is_observation_ir(struct observation* obs) {
	return usnob_is_band_ir(usnob_get_survey_band(obs->survey));
}

int usnob_get_slice(usnob_entry* entry) {
	return (entry->usnob_id >> 24) & 0xFF;
}

int usnob_get_index(usnob_entry* entry) {
	return (entry->usnob_id & 0x00ffffff);
}

int usnob_parse_entry(unsigned char* line, usnob_entry* usnob) {
    int obs;
	int A, S, P, i, j, M, R, Q, y, x, k, e, v, u;
	uint32_t ival;
    uint32_t* uline;

    uline = (uint32_t*)line;

	// bytes 0-3: uint, RA in units of 0.01 arcsec.
	ival = u32_letoh(uline[0]);
	if (ival > (100*60*60*360)) {
		fprintf(stderr, "USNOB: RA should be in [0, %u), but got %u.\n",
				100*60*60*360, ival);
		assert(ival <= (100*60*60*360));
		return -1;
	}
	usnob->ra = arcsec2deg(ival * 0.01);

	// bytes 4-7: uint, SPD (south polar distance) in units of 0.01 arcsec.
	ival = u32_letoh(uline[1]);
	assert(ival <= (100*60*60*180));
	// DEC = south polar distance - 90 degrees
	usnob->dec = arcsec2deg(ival * 0.01) - 90.0;

	// bytes 8-11: uint, packed in base-10:
	//    iPSSSSAAAA
	ival = u32_letoh(uline[2]);
	A = (ival % 10000);
	ival     /= 10000;
	S = (ival % 10000);
	ival     /= 10000;
	P = (ival % 10);
	ival     /= 10;
	i = (ival % 10);

	// A: mu_RA, in units of 0.002 arcsec per year, offset by
	//    -10 arcsec per year.
    // (rewrite in this form to avoid cancellation error for zero)
	usnob->pm_ra = 0.002 * (A - 5000); // -10.0 + (0.002 * A);

	// S: mu_SPD, in units of 0.002 arcsec per year, offset by
	//    -10 arcsec per year.
	// This is a derivative of SPD which is equal to a derivative of DEC.
	usnob->pm_dec = 0.002 * (S - 5000); //-10.0 + (0.002 * S);

	// P: total mu probability, in units of 0.1.
	usnob->pm_prob = 0.1 * P;

	// i: motion catalog flag: 0=no, 1=yes.
	assert((i == 0) || (i == 1));
	usnob->motion_catalog = i;

	// bytes 12-15: uint, packed in base-10:
	//     jMRQyyyxxx
	ival = u32_letoh(uline[3]);
	x = (ival % 1000);
	ival     /= 1000;
	y = (ival % 1000);
	ival     /= 1000;
	Q = (ival % 10);
	ival     /= 10;
	R = (ival % 10);
	ival     /= 10;
	M = (ival % 10);
	ival     /= 10;
	j = (ival % 10);

	// x: sigma_mu_RA, in units of 0.001 arcsec per year.
	usnob->sigma_pm_ra = 0.001 * x;

	// y: sigma_mu_SPD, in units of 0.001 arcsec per year.
	// Again, a derivative of SPD = derivate of DEC.
	usnob->sigma_pm_dec = 0.001 * y;

	// Q: sigma_RA_fit, in units of 0.1 arcsec.
	usnob->sigma_ra_fit = arcsec2deg(0.1 * Q);

	// R: sigma_SPD_fit, in units of 0.1 arcsec.
	usnob->sigma_dec_fit = arcsec2deg(0.1 * R);

	// M: number of detections; in [2, 5] for USNOB stars,
	//   1 for rejected USNOB stars,
	//   0 for Tycho-2 stars.
	usnob->ndetections = M;

	// j: diffraction spike flag: 0=no, 1=yes.
	assert((j == 0) || (j == 1));
	usnob->diffraction_spike = j;

	// bytes 16-19: uint, packed in base-10:
	//     keeevvvuuu
	ival = u32_letoh(uline[4]);
	u = (ival % 1000);
	ival     /= 1000;
	v = (ival % 1000);
	ival     /= 1000;
	e = (ival % 1000);
	ival     /= 1000;
	k = (ival % 10);

	// u: sigma_RA, in units of 0.001 arcsec.
	usnob->sigma_ra = arcsec2deg(0.001 * u);

	// v: sigma_SPD, in units of 0.001 arcsec.
	// Again, equal to sigma_DEC.
	usnob->sigma_dec = arcsec2deg(0.001 * v);

	// e: mean epoch, in 0.1 yr, offset by -1950.
	usnob->epoch = 1950.0 + 0.1 * e;

	// k: YS4.0 correlation flag: 0=no, 1=yes.
    // (if M==0, it's a Tycho star and this has a different meaning)
	assert((M == 0) || (k == 0) || (k == 1));
	usnob->ys4 = (k == 1) ? 1 : 0;

	for (obs=0; obs<5; obs++) {
		int G, S, F, m, C, r, R;

		// bytes 20-23, 24-27, ...: uint, packed in base-10:
		//     GGSFFFmmmm
		ival = u32_letoh(uline[5 + obs]);
		m = (ival % 10000);
		ival     /= 10000;
		F = (ival % 1000);
		ival     /= 1000;
		S = (ival % 10);
		ival     /= 10;
		G = (ival % 100);

		// m: magnitude, in units of 0.01 mag.
		usnob->obs[obs].mag = 0.01 * m;

		// F: field number in the original survey; 1-937.
		if (M >= 2)
			assert(F <= 937);
		usnob->obs[obs].field = F;

		// S: survey number of original survey.
		usnob->obs[obs].survey = S;

		// G: star-galaxy estimate.  0=galaxy, 11=star.
		if (M >= 2) {
            // Hmm, this is triggered by USNOB10/033/b0337.cat
            // byte offset 10282000, observation 2, and MANY others.
			/*
             assert(G <= 11 || G == 19);
             if ((G > 11) && (G != 19)) {
             fprintf(stderr, "USNOB: star/galaxy estimate should be in {[0, 11], 19}, but found %u.\n", G);
             }
             */
		}
		usnob->obs[obs].star_galaxy = G;

		// bytes 40-43, 44-47, ...: uint, packed in base-10:
		//     CrrrrRRRR
		ival = u32_letoh(uline[10 + obs]);
		R = (ival % 10000);
		ival     /= 10000;
		r = (ival % 10000);
		ival     /= 10000;
		C = (ival % 10);

        if (M >= 2 && F == 0) {
            // empty observation.
            usnob->obs[obs].xi_resid  = 0.0;
            usnob->obs[obs].eta_resid = 0.0;
        } else {
            // R: xi residual, in units of 0.01 arcsec, offset by -50 arcsec.
            usnob->obs[obs].xi_resid = arcsec2deg(0.01 * (R - 5000)); // (-50.0 + 0.01 * R);
            // r: eta residual, in units of 0.01 arcesc, offset by -50.
            usnob->obs[obs].eta_resid = arcsec2deg(0.01 * (r - 5000)); // (-50.0 + 0.01 * r);
        }

		// C: source of photometric calibration.
		usnob->obs[obs].calibration = C;

		// bytes 60-63, 64-67, ...: uint
		// index into PMM scan file.
		ival = u32_letoh(uline[15 + obs]);
		assert(ival <= 9999999);
		usnob->obs[obs].pmmscan = ival;
	}

	return 0;
}