# include "dvoshell.h"

// to add a new field to the extractions:
// 1) add the selection for the field below
// 2) add the field to the ParseMeasureField in dbFields.c
// 3) add the field to the measure enum list in dvoshell.h

/* time concepts */
static time_t TimeReference;
static int TimeFormat;

static time_t TimeRefPM;

static CoordTransform *celestial_to_galactic = NULL;
static CoordTransform *celestial_to_ecliptic = NULL;

static int haveGalacticAve = FALSE;
static double GLON_AVE = 0.0;
static double GLAT_AVE = 0.0;

static int haveEclipticAve = FALSE;
static double ELON_AVE = 0.0;
static double ELAT_AVE = 0.0;

static int haveGalacticMeas = FALSE;
static double GLON_MEAS = 0.0;
static double GLAT_MEAS = 0.0;

static int haveEclipticMeas = FALSE;
static double ELON_MEAS = 0.0;
static double ELAT_MEAS = 0.0;

void dbExtractMeasuresInit () {
  TimeRefPM = ohana_date_to_sec ("2000/01/01");
  GetTimeFormat (&TimeReference, &TimeFormat);
}

// define a locally-static transform
int dbExtractMeasuresInitTransform (CoordTransformSystem target) {

  // galactic transform is kept forever
  if (target == COORD_GALACTIC) {
    if (celestial_to_galactic != NULL) return (TRUE);
    celestial_to_galactic = InitTransform (COORD_CELESTIAL, target);
    return (TRUE);
  }

  // ecliptic transform must be updated (is weakly time-dependent)
  if (target == COORD_ECLIPTIC) {
    if (celestial_to_ecliptic != NULL) {
      free (celestial_to_ecliptic);
    }
    celestial_to_ecliptic = InitTransform (COORD_CELESTIAL, target);
    return (TRUE);
  }
  return (FALSE);
}

int dbExtractMeasuresInitAve () {
  haveGalacticAve = FALSE;
  haveEclipticAve = FALSE;
  return (TRUE);
}

int dbExtractMeasuresInitMeas () {
  haveGalacticMeas = FALSE;
  haveEclipticMeas = FALSE;
  return (TRUE);
}

/* return measure.field based on the selection */
double dbExtractMeasures (Average *average, SecFilt *secfilt, Measure *measure, dbField *field) {

  int Nsec;
  double value;
  double ra, dec, x, y, dT;

  Coords *mosaic;
  PhotCode *equiv;

  value = NAN;

  switch (field->ID) {
    case MEAS_MAG: /* magnitudes are already determined above */
      equiv = GetPhotcodeEquivbyCode (measure[0].photcode);

      // we return the magnitude for this measure if:
      if (field->photcode->type == PHOT_MAG) goto valid_photcode;
      if ((field->photcode->type == PHOT_REF) && (measure[0].photcode == field->photcode->code)) goto valid_photcode;
      if ((field->photcode->type == PHOT_DEP) && (measure[0].photcode == field->photcode->code)) goto valid_photcode;

      if ((equiv != NULL) && (field->photcode->type == PHOT_SEC) && (equiv[0].code == field->photcode->code)) goto valid_photcode;
      break;

  valid_photcode:
      switch (field->magMode) {
	case MAG_INST:
	  value = PhotInst (measure);  
	  break;
	case MAG_CAT:
	  value = PhotCat  (measure); 
	  break;
	case MAG_SYS:
	  value = PhotSys  (measure, average, secfilt); 
	  break;
	case MAG_REL:
	  value = PhotRel  (measure, average, secfilt); 
	  break;
	case MAG_CAL:
	  value = PhotCal  (measure, average, secfilt, measure, equiv); 
	  break;
	case MAG_AVE:
	  value = PhotAve  (equiv, average, secfilt); 
	  break;
	case MAG_REF:
	  value = PhotRef  (equiv, average, secfilt, measure); 
	  break;
	case MAG_ERR:
	  if (field->photcode->type == PHOT_SEC) {
	    value = PhotdM  (equiv, average, secfilt); 
	  } else {
	    value = measure[0].dM;
	  }
	  break;
	case MAG_CHISQ:
	  value = PhotXm  (equiv, average, secfilt); 
	  break;
	case MAG_NCODE:
	  if (equiv == NULL) {
	    value = NAN;
	    break;
	  }
	  Nsec = GetPhotcodeNsec (equiv->code);
	  if (Nsec == -1) break;
	  value = secfilt[Nsec].Ncode;
	  break;
	case MAG_NPHOT:
	  if (equiv == NULL) {
	    value = NAN;
	    break;
	  }
	  Nsec = GetPhotcodeNsec (equiv->code);
	  if (Nsec == -1) break;
	  value = secfilt[Nsec].Nused;
	  break;
      }
      break;
    case MEAS_RA: /* OK */
      value = average[0].R - measure[0].dR / 3600.0;
      break;
    case MEAS_DEC: /* OK */
      value = average[0].D - measure[0].dD / 3600.0;
      break;
    case MEAS_RA_AVE: /* OK */
      value = average[0].R;
      break;
    case MEAS_DEC_AVE: /* OK */
      value = average[0].D;
      break;

    case MEAS_GLON:
      if (!haveGalacticMeas) {
	ApplyTransform (&GLON_MEAS, &GLAT_MEAS, average[0].R - measure[0].dR / 3600.0, average[0].D - measure[0].dD / 3600.0, celestial_to_galactic);
	haveGalacticMeas = TRUE;
      }
      value = GLON_MEAS;
      break;
    case MEAS_GLAT:
      if (!haveGalacticMeas) {
	ApplyTransform (&GLON_MEAS, &GLAT_MEAS, average[0].R - measure[0].dR / 3600.0, average[0].D - measure[0].dD / 3600.0, celestial_to_galactic);
	haveGalacticMeas = TRUE;
      }
      value = GLAT_MEAS;
      break;
    case MEAS_ELON:
      if (!haveEclipticMeas) {
	ApplyTransform (&ELON_MEAS, &ELAT_MEAS, average[0].R - measure[0].dR / 3600.0, average[0].D - measure[0].dD / 3600.0, celestial_to_ecliptic);
	haveEclipticMeas = TRUE;
      }
      value = ELON_MEAS;
      break;
    case MEAS_ELAT:
      if (!haveEclipticMeas) {
	ApplyTransform (&ELON_MEAS, &ELAT_MEAS, average[0].R - measure[0].dR / 3600.0, average[0].D - measure[0].dD / 3600.0, celestial_to_ecliptic);
	haveEclipticMeas = TRUE;
      }
      value = ELAT_MEAS;
      break;

    case MEAS_GLON_AVE:
      if (!haveGalacticAve) {
	ApplyTransform (&GLON_AVE, &GLAT_AVE, average[0].R, average[0].D, celestial_to_galactic);
	haveGalacticAve = TRUE;
      }
      value = GLON_AVE;
      break;
    case MEAS_GLAT_AVE:
      if (!haveGalacticAve) {
	ApplyTransform (&GLON_AVE, &GLAT_AVE, average[0].R, average[0].D, celestial_to_galactic);
	haveGalacticAve = TRUE;
      }
      value = GLAT_AVE;
      break;
    case MEAS_ELON_AVE:
      if (!haveEclipticAve) {
	ApplyTransform (&ELON_AVE, &ELAT_AVE, average[0].R, average[0].D, celestial_to_ecliptic);
	haveEclipticAve = TRUE;
      }
      value = ELON_AVE;
      break;
    case MEAS_ELAT_AVE:
      if (!haveEclipticAve) {
	ApplyTransform (&ELON_AVE, &ELAT_AVE, average[0].R, average[0].D, celestial_to_ecliptic);
	haveEclipticAve = TRUE;
      }
      value = ELAT_AVE;
      break;

    case MEAS_RA_AVE_ERR: /* OK */
      value = average[0].dR;
      break;
    case MEAS_DEC_AVE_ERR: /* OK */
      value = average[0].dD;
      break;
    case MEAS_U_RA: /* OK */
      value = average[0].uR;
      break;
    case MEAS_U_DEC: /* OK */
      value = average[0].uD;
      break;
    case MEAS_U_RA_ERR: /* OK */
      value = average[0].duR;
      break;
    case MEAS_U_DEC_ERR: /* OK */
      value = average[0].duD;
      break;
    case MEAS_PAR: /* OK */
      value = average[0].R;
      break;
    case MEAS_PAR_ERR: /* OK */
      value = average[0].D;
      break;
    case MEAS_XP: /* OK */
      value = average[0].Xp;
      break;
    case MEAS_NMEAS: /* OK */
      value = average[0].Nmeasure;
      break;
    case MEAS_NMISS: /* OK */
      value = average[0].Nmissing;
      break;
    case MEAS_OBJFLAGS: /* OK */
      value = average[0].code;
      break;
    // note that these represent the ra displacement relative to the average, not 
    // the error.
    case MEAS_RA_OFFSET: /* OK */
      value = measure[0].dR;
      break;
    case MEAS_DEC_OFFSET: /* OK */
      value = measure[0].dD;
      break;
    case MEAS_RA_FIT_OFFSET: /* OK */
      dT = (measure[0].t - TimeRefPM) / (86400*365.25);
      value = average[0].uR * dT + measure[0].dR;
      break;
    case MEAS_DEC_FIT_OFFSET: /* OK */
      dT = (measure[0].t - TimeRefPM) / (86400*365.25);
      value = average[0].uD * dT + measure[0].dD;
      break;
    case MEAS_RA_OFFSET_ERR: /* OK */
      value = NAN;
      break;
    case MEAS_DEC_OFFSET_ERR: /* OK */
      value = NAN;
      break;
    case MEAS_AIRMASS: /* OK */
      value = measure[0].airmass;
      break;
    case MEAS_AZ: /* OK */
      value = measure[0].az;
      break;
    case MEAS_EXPTIME: /* OK */
      value = pow (10.0, measure[0].dt * 0.4);
      break;
    case MEAS_PHOTCODE: /* OK */
      value = measure[0].photcode;
      break;
    case MEAS_TIME: /* OK */
      value = TimeValue (measure[0].t, TimeReference, TimeFormat);
      break;
    case MEAS_FWHM: /* OK */
      value = 0.01*(measure[0].FWx + measure[0].FWy) / 2.0;
      break;
    case MEAS_FWHM_MAJ: /* OK */
      value = 0.01*measure[0].FWx;
      break;
    case MEAS_FWHM_MIN: /* OK */
      value = 0.01*measure[0].FWy;
      break;
    case MEAS_THETA: /* OK */
      value = measure[0].theta;
      break;
    case MEAS_DOPHOT: /* OK */
      value = measure[0].dophot;
      break;
    case MEAS_DB_FLAGS: /* ? */
      value = measure[0].dbFlags;
      break;
    case MEAS_PHOT_FLAGS: /* ? */
      value = measure[0].photFlags;
      break;
    case MEAS_XCCD: /* OK */
/* I need to perform this conversion for ELIXIR and LONEOS formats on load */      
# if 1
      value = measure[0].Xccd;
# else
      { 
	Image *image;
	ra  = average[0].R - measure[0].dR / 3600.0;
	dec = average[0].D - measure[0].dD / 3600.0;
	image = MatchImage (measure[0].t, measure[0].photcode);
	if (image == NULL) break;
	RD_to_XY (&x, &y, ra, dec, &image[0].coords);
	value = x;
      }
# endif
      break;
    case MEAS_YCCD: /* OK */
/* I need to perform this conversion for ELIXIR and LONEOS formats on load */      
# if 1
      value = measure[0].Yccd;
# else
      {
	Image *image;
	ra  = average[0].R - measure[0].dR / 3600.0;
	dec = average[0].D - measure[0].dD / 3600.0;
	image = MatchImage (measure[0].t, measure[0].photcode);
	if (image == NULL) break;
	RD_to_XY (&x, &y, ra, dec, &image[0].coords);
	value = y;
      }
# endif
      break;
    case MEAS_XMOSAIC: /* OK */
      ra  = average[0].R - measure[0].dR / 3600.0;
      dec = average[0].D - measure[0].dD / 3600.0;
      mosaic = MatchMosaic (measure[0].t, measure[0].photcode);
      if (mosaic == NULL) break;
      RD_to_XY (&x, &y, ra, dec, mosaic);
      value = x;
      break;
    case MEAS_YMOSAIC: /* OK */
      ra  = average[0].R - measure[0].dR / 3600.0;
      dec = average[0].D - measure[0].dD / 3600.0;
      mosaic = MatchMosaic (measure[0].t, measure[0].photcode);
      if (mosaic == NULL) break;
      RD_to_XY (&x, &y, ra, dec, mosaic);
      value = y;
      break;

    case MEAS_SKY: /* OK */
      value = measure[0].Sky;
      break;
    case MEAS_dSKY: /* OK */
      value = measure[0].dSky;
      break;
    case MEAS_DET_ID: /* OK */
      value = measure[0].detID;
      break;
    case MEAS_OBJ_ID: /* OK */
      value = measure[0].averef;
      break;
    case MEAS_IMAGE_ID: /* OK */
      value = measure[0].imageID;
      break;
    case MEAS_PSF_QF: /* OK */
      value = measure[0].qPSF;
      break;
    case MEAS_PSF_CHISQ: /* OK */
      value = measure[0].psfChisq;
      break;
    case MEAS_CR_NSIGMA: /* OK */
      value = measure[0].crNsigma;
      break;
    case MEAS_EXT_NSIGMA: /* OK */
      value = measure[0].extNsigma;
      break;
    case MEAS_STARGAL: /* OK */
      value = measure[0].stargal;
      break;
  }
  return (value);
}

/** the mosaic entries do not use the registered mosaic found 
    by MatchImage (via FindMosaicForImage).  Rather, they use
    a coordinate frame saved by SetImageSelection 
**/
