Index: /branches/eam_branches/ipp-20150326/psastro/src/psastroConvert.c
===================================================================
--- /branches/eam_branches/ipp-20150326/psastro/src/psastroConvert.c	(revision 38029)
+++ /branches/eam_branches/ipp-20150326/psastro/src/psastroConvert.c	(revision 38030)
@@ -335,4 +335,16 @@
   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");
Index: /branches/eam_branches/ipp-20150326/psastro/src/psastroGalaxyShapeErrors.c
===================================================================
--- /branches/eam_branches/ipp-20150326/psastro/src/psastroGalaxyShapeErrors.c	(revision 38029)
+++ /branches/eam_branches/ipp-20150326/psastro/src/psastroGalaxyShapeErrors.c	(revision 38030)
@@ -12,4 +12,6 @@
 
 # include "psastroInternal.h"
+
+double psastroNormalizeAngleToMidpoint (double angle, double midpoint);
 
 bool psastroGalaxyShapeErrors (psMetadata *recipe, pmReadout *readout) {
@@ -71,8 +73,10 @@
 	    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,4 +93,21 @@
 	    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,10 +122,15 @@
 		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,6 +155,22 @@
     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;
+}
+
