/** @file psastroMosaicZeroPoint.c
 *
 *  @brief 
 *
 *  @ingroup libpsastro
 *
 *  @author IfA
 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2009-02-07 02:03:34 $
 *  Copyright 2009 Institute for Astronomy, University of Hawaii
 */

# include "psastroInternal.h"

bool psastroZeroPoint (pmConfig *config) {

    float zeropt, exptime;
    pmChip *chip = NULL;
    pmCell *cell = NULL;
    pmReadout *readout = NULL;

    // select the current recipe
    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    if (!recipe) {
	psError(PSASTRO_ERR_CONFIG, false, "Can't find PSASTRO recipe!\n");
	return false;
    }

    // select the input data sources
    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    if (!input) {
	psError(PSASTRO_ERR_CONFIG, false, "Can't find input data!\n");
	return false;
    }
    pmFPA *fpa = input->fpa;

    pmFPAview *view = pmFPAviewAlloc (0);

    // given the existing per-chip astrometry, determine matches between raw and ref stars
    // is this needed? yes, if we didn't do SingleChip astrometry first
    if (!psastroMosaicSetMatch (fpa, recipe, 0)) {
        psError(PSASTRO_ERR_UNKNOWN, false, "failed to match raw and ref stars for mosaic");
        return false;
    }

    // XXX eventually: look up the photcode info from the recipe, use the ZEROPT and A_0 terms
    // to determine transparency = Mref - Mraw - 2.5*log10(exptime) - ZEROPT - A_0*Color_ref
    // to get there, we need: a) getstar to write out the color-term magnitudes and b)  the
    // pmAstromObj to carry a color value as well.

    // really error-out here?  or just skip?
    if (!psastroZeroPointFromRecipe (&zeropt, &exptime, fpa, recipe)) {
	psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe");
	return false;
    }

    // this loop selects the matched stars for all chips
    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
        if (!chip->process || !chip->file_exists) { continue; }
	if (!chip->toFPA) { continue; }

	while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
            if (!cell->process || !cell->file_exists) { continue; }

	    // process each of the readouts
	    // XXX there can only be one readout per chip, right?
	    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
		if (! readout->data_exists) { continue; }

		// calculate dMag for the matched stars
		psastroZeroPointReadout (readout, zeropt, exptime);

	    }
	}
    }

    psFree (view);
    return true;
}

/**
 * we measure <dMag> and \sigma_dMag and write them to the header
 */
bool psastroZeroPointReadout(pmReadout *readout, float zeropt, float exptime) {

    bool status;

    // select the raw objects for this readout
    psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.RAWSTARS");
    if (rawstars == NULL) return false;

    // select the raw objects for this readout
    psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS");
    if (refstars == NULL) return false;

    psArray *matches = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.MATCH");
    if (matches == NULL) return false;

    psVector *dMag  = psVectorAllocEmpty (100, PS_TYPE_F32);

    int Npts = 0;
    for (int i = 0; i < matches->n; i++) {

      pmAstromMatch *match = matches->data[i];

      pmAstromObj *raw = rawstars->data[match->raw];

      pmAstromObj *ref = refstars->data[match->ref];

      dMag->data.F32[Npts] = ref->Mag - raw->Mag - 2.5*log10(exptime);
      psVectorExtend (dMag, 100, 1);
      Npts++;
    }

    psTrace ("psModules.astrom", 4, "Npts: %d\n", Npts);

    if (Npts < 3) {
      fprintf (stderr, "zero point NaN +/- NaN\n");
      psFree (dMag);
      return false;
    }

    // stats structure and mask for use in measuring the clipping statistic
    // this analysis has too few data points to use the robust median method
    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    psVectorStats (stats, dMag, NULL, NULL, 0);
    fprintf (stderr, "zero point %f +/- %f using %d stars; transparency %f\n", stats->clippedMean, stats->clippedStdev, Npts, zeropt - stats->clippedMean);

    psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");

    psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_OBS", PS_META_REPLACE, "measured zero point", stats->clippedMean);
    psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_REF", PS_META_REPLACE, "measured zero point", zeropt);
    psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_ERR", PS_META_REPLACE, "measured zero point", stats->clippedStdev);
    psMetadataAddF32 (header, PS_LIST_TAIL, "ZPT_OFF", PS_META_REPLACE, "measured zero point", zeropt - stats->clippedMean);

    psFree (dMag);
    psFree (stats);
    return true;
}

# define ESCAPE(MSG) { \
  psLogMsg ("psastro", PS_LOG_INFO, MSG); \
  return false; }

bool psastroZeroPointFromRecipe (float *zeropt, float *exptime, pmFPA *fpa, psMetadata *recipe) {

    bool status;

    // select the filter; default to fixed photcode and mag limit otherwise
    char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID");
    if (!status) ESCAPE ("missing FPA.FILTER in concepts");

    *exptime = psMetadataLookupF32 (&status, fpa->concepts, "FPA.EXPOSURE");
    if (!status) ESCAPE ("missing FPA.EXPOSURE in concepts");

    // we need to select the PHOTCODE.DATA folder that matches our filter
    psMetadataItem *item = psMetadataLookup (recipe, "PHOTCODE.DATA");
    if (!item) ESCAPE ("PHOTCODE.DATA folders missing");
    if (item->type != PS_DATA_METADATA_MULTI) ESCAPE ("PHOTCODE.DATA not a multi");

    // PHOTCODE.DATA is a multi of metadata items 
    psListIterator *iter = psListIteratorAlloc(item->data.list, PS_LIST_HEAD, false);

    psMetadataItem *refItem = NULL;
    while ((refItem = psListGetAndIncrement (iter))) {
	if (refItem->type != PS_DATA_METADATA) ESCAPE ("PHOTCODE.DATA entry is not a metadata folder");
    
	char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER");
	if (!status) {
	    psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
	    continue;
	}

	// does this entry match the current filter?
	if (strcmp (refFilter, filter)) continue;

	psLogMsg ("psastro", PS_LOG_DETAIL, "PHOTCODE.DATA found for filter %s", filter);

	*zeropt = psMetadataLookupF32 (&status, refItem->data.md, "ZEROPT");
	if (!status) {
	    psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER");
	    continue;
	}
	psFree (iter);
	return true;
    }
    psFree (iter);
    return false;
}
