# include "psastroInternal.h"
# define NONLIN_TOL 0.001 /* tolerance in pixels */

XXX add a function to supply for all entries, as opposed to fixing the poor astrometry ones.
bool psastroEnforceChips (pmConfig *config, psMetadata *recipe) {

  bool status;

  bool fixChips = psMetadataLookupBool (&status, config->arguments, "PSASTRO.FIX.CHIPS");
  if (!status) {
      fixChips = psMetadataLookupBool (&status, recipe, "PSASTRO.FIX.CHIPS");
  }
  if (!fixChips) return true;

  // identify reference astrometry table
  // if not defined, correction was not requested; skip step
  pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.MODEL");
  if (!astrom) return true;

  // 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;
  }

  // make sure the astrometry model is loaded
  // de-activate all files except PSASTRO.MODEL
  // XXX do I need to supply the input->fpa->concepts to the astrom->fpa->concepts?
  psFree (astrom->fpa->concepts);
  astrom->fpa->concepts = psMemIncrRefCounter (input->fpa->concepts);
  pmFPAfileActivate (config->files, false, NULL);
  pmFPAfileActivate (config->files, true, "PSASTRO.MODEL");

  pmFPAview *view = pmFPAviewAlloc (0);

  // files associated with the science image
  if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
      psError (PS_ERR_IO, false, "Can't load the astrometry model file");
      return false;
  }

  // loop over all chips, replace input astrometry elements with those from astrom 
  pmChip *obsChip = NULL;
  while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
    psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }

    // set the chip astrometry using the astrom file
    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);

    psFree (obsChip->toFPA);
    psFree (obsChip->fromFPA);

    // XXX can we do a straight copy?  I think so, but we may need to apply a scaling factor based
    // on the toSky or toTPA values
    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);

    // XXX test to write out adjusted header
    pmAstromWriteBilevelChip (obsChip->hdu->header, obsChip, NONLIN_TOL);
  }

  psFree (input->fpa->toSky);
  psFree (input->fpa->toTPA);
  psFree (input->fpa->fromTPA);
  input->fpa->toSky   = psMemIncrRefCounter (astrom->fpa->toSky);
  input->fpa->toTPA   = psMemIncrRefCounter (astrom->fpa->toTPA);
  input->fpa->fromTPA = psMemIncrRefCounter (astrom->fpa->fromTPA);

  psMetadata *updates = psMetadataAlloc();
  pmAstromWriteBilevelMosaic (updates, input->fpa, NONLIN_TOL);
  psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", updates);
  psFree (updates);

  psFree (view);
  return true;
}

# if (0)

  // physical pixel scale in microns per pixel
  double pixelTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.TOLERANCE");
  if (!status) {
    psError(PS_ERR_IO, true, "Failed to lookup pixel tolerance"); 
    psFree (view);
    return false; 
  } 
  double angleTol = psMetadataLookupF32 (&status, recipe, "PSASTRO.ANGLE.TOLERANCE");
  if (!status) {
    psError(PS_ERR_IO, true, "Failed to lookup angle tolerance"); 
    psFree (view);
    return false; 
  } 

  pmChip *obsChip = NULL;

  // loop over all chips
  while ((obsChip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
    psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, obsChip->file_exists, obsChip->process);
    if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }

    // set the chip astrometry using the astrom file
    pmChip *refChip = pmFPAviewThisChip (view, astrom->fpa);

    // bad Astrometry test:  ref pixel or angle outside nominal

    psPlane refPixel = {0.0, 0.0, 0.0, 0.0};
    psPlane obsCoord, refCoord;

    // find location of 0,0 pixel in focal plane coords for this chip
    psPlaneTransformApply (&obsCoord, obsChip->toFPA, &refPixel);

    // find location of 0,0 pixel in focal plane coords for ref chip
    psPlaneTransformApply (&refCoord, refChip->toFPA, &refPixel);
    
    psPlane offPixel = {0.0, 0.0, 0.0, 0.0};
    psPlane obsOffPt, refOffPt;

    // find location of 0,0 pixel in focal plane coords for this chip
    psPlaneTransformApply (&obsOffPt, obsChip->toFPA, &offPixel);

    // find location of 0,0 pixel in focal plane coords for ref chip
    psPlaneTransformApply (&refOffPt, refChip->toFPA, &offPixel);
    
    double obsAngle = atan2 (obsOffPt.y - obsCoord.y, obsOffPt.x - obsCoord.x);
    double refAngle = atan2 (refOffPt.y - refCoord.y, refOffPt.x - refCoord.x);

    bool badAstrom = false;
    badAstrom |= fabs(obsCoord.x - refCoord.x) > pixelTol;
    badAstrom |= fabs(obsCoord.y - refCoord.y) > pixelTol;
    badAstrom |= fabs(obsAngle   - refAngle)   > angleTol;

    if (!badAstrom) continue;

    psFree (obsChip->toFPA);
    psFree (obsChip->fromFPA);

    // XXX can we do a straight copy?  I think so, but we may need to apply a scaling factor based
    // on the toSky or toTPA values
    obsChip->toFPA   = psMemIncrRefCounter (refChip->toFPA);
    obsChip->fromFPA = psMemIncrRefCounter (refChip->fromFPA);
  }
  psFree (view);
  return true;
}

// XXX for this function to work, I need to:
// 1 *) define badAstrom (ref pixel too far from nominal?) (ref angle too far off?)
// 2) load the astrom fpa in the pmFPAfileIO stage 
// 3) allow for this operation to be optional
# endif
