# include "psastroInternal.h"

// this function loads the header WCS astrometry terms into the fpa terms and applies the
// astrometry to the detected objects.

// this function assumes the initial astrometry arrives in the form of WCS keywords in the
// headers corresponding to the chips.

bool psastroAstromGuess (int *nStars, pmConfig *config) {

    bool newFPA = true;
    bool status = false;
    double RAmin  = +FLT_MAX;
    double RAmax  = -FLT_MAX;
    double DECmin = +FLT_MAX;
    double DECmax = -FLT_MAX;

    double RAminSky = NAN;
    double RAmaxSky = NAN;

    pmChip *chip = NULL;
    pmCell *cell = NULL;
    pmReadout *readout = NULL;

    *nStars = 0;

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

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

    // physical pixel scale in microns per pixel
    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    if (!status) {
	psError(PS_ERR_IO, true, "Failed to lookup pixel scale"); 
	return false; 
    } 

    pmFPAview *view = pmFPAviewAlloc (0);
    pmFPA *fpa = input->fpa;

    // load mosaic-level astrometry?
    bool bilevelAstrometry = false;
    pmHDU *phu = pmFPAviewThisPHU (view, fpa);
    if (phu) {
      char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
      if (ctype) {
	bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
      }
    }
    if (bilevelAstrometry) {
      pmAstromReadBilevelMosaic (fpa, phu->header);
    } 

    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 || !chip->data_exists) { continue; }

        // read WCS data from the corresponding header
        pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
	if (bilevelAstrometry) {
	  pmAstromReadBilevelChip (chip, hdu->header); 
	} else {
	  pmAstromReadWCS (fpa, chip, hdu->header, pixelScale);
	}

        if (newFPA) {
            newFPA = false;
	    while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
	    while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
            RAminSky = fpa->toSky->R - M_PI;
            RAmaxSky = fpa->toSky->R + M_PI;
        }

        // apply the new WCS guess data to all of the data in the readouts
        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
            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                if (! readout->data_exists) { continue; }

                psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
                if (rawstars == NULL) { continue; }

		*nStars += rawstars->n;
                for (int i = 0; i < rawstars->n; i++) {
                    pmAstromObj *raw = rawstars->data[i];

                    psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
                    psPlaneTransformApply (raw->TP, fpa->toTPA, raw->FP);
                    psDeproject (raw->sky, raw->TP, fpa->toSky);

                    // rationalize ra to sky range centered on boresite
                    while (raw->sky->r < RAminSky) raw->sky->r += 2.0*M_PI;
                    while (raw->sky->r > RAmaxSky) raw->sky->r -= 2.0*M_PI;

                    RAmin = PS_MIN (raw->sky->r, RAmin);
                    RAmax = PS_MAX (raw->sky->r, RAmax);

                    DECmin = PS_MIN (raw->sky->d, DECmin);
                    DECmax = PS_MAX (raw->sky->d, DECmax);
                }

		// dump or plot the resulting projected positions
		if (psTraceGetLevel("psastro.dump") > 0) {
		    psastroDumpRawstars (rawstars, fpa, chip);
		}

		if (psTraceGetLevel("psastro.plot") > 0) {
		    psastroPlotRawstars (rawstars, fpa, chip);
		}
            }
        }
    }

    // how many total sources are available to us?
    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
    if (*nStars == 0) {
	psLogMsg ("psastro", 2, "no sources available for astrometry\n");
	psFree (view);
	return true;
    }

    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
	      DEG_RAD*RAmin, DEG_RAD*DECmin, 
	      DEG_RAD*RAmax, DEG_RAD*DECmax);

    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MAX",  PS_META_REPLACE, "", RAmax);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", DECmin);
    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", DECmax);

    psFree (view);
    return true;
}

/* coordinate frame hierachy
   pixels (on a given readout)
   cell
   chip
   FP (focal plane)
   TP (tangent plane)
   sky (ra, dec)
*/
