Changeset 38030

Timestamp:

Mar 28, 2015, 5:32:49 AM (11 years ago)

Author:

eugene

Message:

skip non-gpc1 chips when applying KH correction; in galaxy shape error analysis, measure stdev of theta for a domain of theta - pi : theta + pi; add flags for cases where galaxy shape error analysis fails; add header keyword to note errors are in major,minor,theta not PAR[Sxx,Sxy,Syy]

Location:

branches/eam_branches/ipp-20150326/psastro/src

Files:

• psastroConvert.c (modified) (1 diff)
• psastroGalaxyShapeErrors.c (modified) (5 diffs)

branches/eam_branches/ipp-20150326/psastro/src/psastroConvert.c

@@ -335 +335 @@
   pmChip *chip = readout->parent->parent;
   char *chipName = psMetadataLookupStr (NULL, chip->concepts, "CHIP.NAME");
+
+  // skip stacks
+  if (!strcmp(chipName, "SkyChip")) {
+    psLogMsg ("psastro.correctKH", PS_LOG_DETAIL, "skipping KH correction: not a gpc1 chip\n");
+    return true;
+  }
+
+  // only try to address gpc1 chips (should probably check the camera)
+  if (strncmp(chipName, "XY", 2)) {
+    psLogMsg ("psastro.correctKH", PS_LOG_DETAIL, "skipping KH correction: not a gpc1 chip\n");
+    return true;
+  }
 
   psAssert (strlen(chipName) == 4, "error in chip name");

branches/eam_branches/ipp-20150326/psastro/src/psastroGalaxyShapeErrors.c

@@ -12 +12 @@
 
 # include "psastroInternal.h"
+
+double psastroNormalizeAngleToMidpoint (double angle, double midpoint);
 
 bool psastroGalaxyShapeErrors (psMetadata *recipe, pmReadout *readout) {
@@ -71 +73 @@
             if (!(isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXY]) && isfinite(PAR[PM_PAR_SYY]))) {
               // set a mask bit
+              model->flags |= PM_MODEL_STATUS_NAN_SHAPE;
               continue;
             }
             if (!(isfinite(dPAR[PM_PAR_SXX]) && isfinite(dPAR[PM_PAR_SXY]) && isfinite(dPAR[PM_PAR_SYY]))) {
               // set a (different) mask bit
+              model->flags |= PM_MODEL_STATUS_NAN_SHAPE_ERR;
               continue;
             }
@@ -89 +93 @@
             bool useReff = model->class->useReff;
             psEllipseAxes axes;
+
+            float SxxRef = PAR[PM_PAR_SXX];
+            float SxyRef = PAR[PM_PAR_SXY];
+            float SyyRef = PAR[PM_PAR_SYY];
+            if (!(isfinite(SxxRef) && isfinite(SxyRef) && isfinite(SyyRef))) {
+              // skip objects which are already poor ellipses
+              model->flags |= PM_MODEL_STATUS_NAN_SHAPE;
+              continue;
+            }
+            pmModelParamsToAxes (&axes, SxxRef, SxyRef, SyyRef, useReff);
+            if (!(isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta))) {
+              // skip objects which are already poor ellipses
+              model->flags |= PM_MODEL_STATUS_NAN_SHAPE;
+              continue;
+            }
+            double ThetaRef = axes.theta;
+
             for (int k = 0; k < nSample; k++) {
 
@@ -101 +122 @@
                 if (!(isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta))) continue;
 
+                // theta should be in range ThetaRef - PI : ThetaRef + PI
+                double ThetaNorm = psastroNormalizeAngleToMidpoint (axes.theta, ThetaRef);
+
                 psVectorAppend (majorValues, axes.major);
                 psVectorAppend (minorValues, axes.minor);
-                psVectorAppend (thetaValues, axes.theta);
+                psVectorAppend (thetaValues, ThetaNorm);
             }
             if (majorValues->n == 0) {
                 psWarning ("failed mc error search");
+                model->flags |= PM_MODEL_STATUS_FAIL_SHAPE_ERR;
+                continue;
             }
 
@@ -129 +155 @@
     psFree (thetaValues);
 
+    psMetadata *header = psMetadataLookupMetadata (&status, readout->analysis, "PSASTRO.HEADER");
+    psMetadataAddStr (header, PS_LIST_TAIL, "GALAXY.SHAPE.ERRORS", PS_META_REPLACE, "error analysis mode", "MONTE CARLO");
+
     psLogMsg ("psastro.galaxy", PS_LOG_INFO, "monte carlo shape errors for %d of %d objects: %f sec\n", nObjects, (int) sources->n, psTimerMark ("psastro.galaxy.shape.errors"));
     return true;
 }
 
+// these angles are in radians
+double psastroNormalizeAngleToMidpoint (double angle, double midpoint) {
+
+    double x = cos(angle);
+    double y = sin(angle);
+
+    double result = atan2 (y, x);
+    if (result < midpoint - M_PI) result += 2.0*M_PI;
+    if (result > midpoint + M_PI) result -= 2.0*M_PI;
+
+    return result;
+}
+