Changeset 36837 for trunk/psastro/src/psastroConvert.c

Timestamp:

Jun 8, 2014, 9:27:56 AM (12 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-20140423: add Koppehoefer correction; make sure psastro failures are noted clearl in the output header; do not load the model if we choose skipastro

File:

• trunk/psastro/src/psastroConvert.c (modified) (7 diffs)

trunk/psastro/src/psastroConvert.c

@@ -16 +16 @@
 static psArray *chooseStars(psArray *inStars, char *listName, psArray *sources, psVector *index, int nMax, float iMagMin, float iMagMax, pmSourceMode skip);
 
-bool psastroConvertFPA (pmFPA *fpa, psMetadata *recipe) {
+bool psastroConvertFPA (pmConfig *config, pmFPA *fpa, psMetadata *recipe) {
 
     pmChip *chip;
@@ -35 +35 @@
                 if (! readout->data_exists) { continue; }
 
-                psastroConvertReadout (readout, recipe);
+                if (!psastroConvertReadout (config, view, readout, recipe)) {
+                  psError(PSASTRO_ERR_CONFIG, true, "problem converting readout\n");
+                  return false;
+                }
             }
         }
@@ -44 +47 @@
 
 // pass/apply the WCS information?
-bool psastroConvertReadout (pmReadout *readout, psMetadata *recipe) {
+bool psastroConvertReadout (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe) {
 
     bool status;
@@ -57 +60 @@
     // convert the pmSource objects into pmAstromObj objects (drop !STAR and SATSTAR?)
     psArray *inStars = pmSourceToAstromObj (sources);
+
+    // apply Koppenhoefer correction if needed
+    if (!psastroCorrectKH (config, view, readout, recipe, inStars)){
+      psError(PSASTRO_ERR_CONFIG, true, "failed to correct Koppenhoefer Effect\n");
+      return false;
+    }
 
     // sort in ascending magnitude order
@@ -183 +192 @@
         psF32 *dPAR = model->dparams->data.F32;
 
+        // Note from EAM: I measured the KH correction from the nightly science database
+        // The KH effect was active until 2011/05/11 when the camera voltages were modified
+        // Between 2011/09/06 and 2013/04/30, there was a bug in pmPSF_ModelToAxes with the 
+        // result that the reported values were too large by a factor of sqrt(2).  Fortunately,
+        // this did not affect the data in the nightly science DVO used to measure the
+        // effect (all affected data was processed BEFORE the bug was introduced), and it
+        // does not affect any of the PV2 or PV3 data processed AFTER the bug was fixed.
+        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->type);
+
         pmAstromObj *obj = pmAstromObjAlloc ();
 
@@ -192 +210 @@
         obj->Mag = source->psfMag;
         obj->dMag = source->psfMagErr;
+        obj->SBinst = source->psfMag + 5.0*log10(axes.major);
 
         // XXX do we have the information giving the readout and cell offset?
@@ -221 +240 @@
 }
 
+static float pltScale[] = {
++0.0,
+-0.255347177386,
+-0.255585625172,
+-0.255733260512,
+-0.255747352242,
+-0.255611001253,
+-0.255311933756,
+ +0.0,
+ +0.0,
+ +0.0,
+-0.255296201766,
+-0.255688636720,
+-0.256161080003,
+-0.256464074850,
+-0.256324860930,
+-0.256147363186,
+-0.255734743476,
+-0.255422729909,
+ +0.0,
+ +0.0,
+-0.255508323163,
+-0.256142899513,
+-0.255493895173,
+-0.256997376740,
+-0.256935591102,
+-0.256644588679,
+-0.256123321533,
+-0.255515242517,
+ +0.0,
+ +0.0,
+-0.255633033872,
+-0.256356916189,
+-0.256929253668,
+-0.257250442505,
+-0.257260403216,
+-0.256962024927,
+-0.256387320280,
+-0.255681164920,
+ +0.0,
+ +0.0,
+-0.255671992302,
+-0.256371312678,
+-0.256936249495,
+-0.257237892330,
+-0.257248004675,
+-0.256911942035,
+-0.256365648448,
+-0.255632205397,
+ +0.0,
+ +0.0,
+-0.255455072403,
+-0.256080205411,
+-0.256570428520,
+-0.256908982217,
+-0.256899240971,
+-0.256607972383,
+-0.256107793808,
+-0.255482232273,
+ +0.0,
+ +0.0,
+-0.255237901688,
+-0.255669494152,
+-0.256101093233,
+-0.256372770309,
+-0.256373448133,
+-0.256118197262,
+-0.255689391643,
+-0.255279636681,
+ +0.0,
+ +0.0,
+ +0.0,
+-0.255286350489,
+-0.255517200172,
+-0.255689459741,
+-0.255668706447,
+-0.255539349556,
+-0.255363308907,
+};
+
+bool psastroCorrectKH (pmConfig *config, pmFPAview *view, pmReadout *readout, psMetadata *recipe, psArray *inStars) {
+
+  bool status;
+
+  psAssert (readout, "missing readout");
+  psAssert (readout->parent, "missing cell");
+  psAssert (readout->parent->parent, "missing chip");
+
+  // should we stay or should we go?
+  bool apply = psMetadataLookupBool (&status, recipe, "KH.CORRECT.APPLY.EXP");
+  if (!apply) return true;
+
+  // get the chip ID from the chip name:
+  pmChip *chip = readout->parent->parent;
+  char *chipName = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+
+  psAssert (strlen(chipName) == 4, "error in chip name");
+  psAssert (chipName[0] == 'X', "error in chip name");
+  psAssert (chipName[1] == 'Y', "error in chip name");
+  
+  int chipID = atoi (&chipName[2]);
+
+  // XXX hardwired list of chips to correct
+  // some notes:
+  // XY24 surpisingly does not need correction
+  // XY27 has poor astrometry with a weird correction for all mags
+  // XY36 surpisingly does not need correction
+  // XY67 has poor astrometry with a weird correction for all mags
+
+  // correct this chip? (raise an error for unexpected chipID values)
+  switch (chipID) {
+    case  1: // do not correct
+    case  2: // do not correct
+    case  3: // do not correct
+    case 10: // do not correct
+    case 11: // do not correct
+    case 12: // do not correct
+    case 13: // do not correct
+    case 20: // do not correct
+    case 21: // do not correct
+    case 22: // do not correct
+    case 23: // do not correct
+    case 24: // do not correct
+    case 27: // do not correct
+    case 30: // do not correct
+    case 31: // do not correct
+    case 32: // do not correct
+    case 33: // do not correct
+    case 36: // do not correct
+    case 44: // do not correct
+    case 45: // do not correct
+    case 46: // do not correct
+    case 47: // do not correct
+    case 54: // do not correct
+    case 55: // do not correct
+    case 56: // do not correct
+    case 57: // do not correct
+    case 64: // do not correct
+    case 65: // do not correct
+    case 66: // do not correct
+    case 67: // do not correct
+    case 74: // do not correct
+    case 75: // do not correct
+    case 76: // do not correct
+      return true;
+    case  4: // correct
+    case  5: // correct
+    case  6: // correct
+    case 14: // correct
+    case 15: // correct
+    case 16: // correct
+    case 17: // correct
+    case 25: // correct
+    case 26: // correct
+    case 34: // correct
+    case 35: // correct
+    case 37: // correct
+    case 40: // correct
+    case 41: // correct
+    case 42: // correct
+    case 43: // correct
+    case 50: // correct
+    case 51: // correct
+    case 52: // correct
+    case 53: // correct
+    case 60: // correct
+    case 61: // correct
+    case 62: // correct
+    case 63: // correct
+    case 71: // correct
+    case 72: // correct
+    case 73: // correct
+      break;
+    default:
+      psAbort ("chipID is invalid");
+  }
+
+  // grab the KH correction file
+  pmFPAfile *KHfile = psMetadataLookupPtr (NULL, config->files, "PSASTRO.KH.CORRECT");
+  if (!KHfile) {
+    psError(PM_ERR_CONFIG, false, "KH correction file not found");
+    return false;
+  }
+
+  // grab the corresponding chip
+  pmChip *KHchip = pmFPAviewThisChip (view, KHfile->fpa);
+  psAssert (KHchip, "found KH file, but not chip?");
+
+  // grab the correction spline
+  KHcorrectData *spline = psMetadataLookupPtr (&status, KHchip->analysis, "KH.CORRECT");
+  if (!spline) {
+    psError(PM_ERR_CONFIG, false, "KH correction not found");
+    return false;
+  }
+
+  // the hard-wired array of plate-scales per chip above come from 
+  // a single smf.  in there, a negative plate scale means parity is flipped
+  // on the sky.
+  float myPltScale = fabs(pltScale[chipID]);
+
+  // apply the correction to the detections
+  for (int i = 0; i < inStars->n; i++) {
+    pmAstromObj *obj = inStars->data[i];
+
+    float Xraw = obj->pix->x;
+    float SBinst = obj->SBinst;
+
+    float dX = KHcorrectApply (spline, SBinst);
+
+    float Xfix = Xraw + dX / myPltScale;
+
+    // note that we carry around pix, cell, chip but the real analysis only operates on chip
+    // if in the future we add transformations between chip->cell->pix, then we will need to
+    // make these consistent as well.
+    obj->pix->x = Xfix;
+    obj->cell->x = Xfix;
+    obj->chip->x = Xfix;
+  }
+
+  return true;
+
+}
+