Changeset 41892 for trunk/psModules

Timestamp:

Nov 4, 2021, 6:05:18 PM (5 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-dev-20210817 (pmPattern updates, new inverse transform extra orders api, forward transform uses ORD, pmSourceIO_CMF.c.in conversion to pmFitsTableNew)

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/astrom/pmAstrometryDistortion.c (modified) (3 diffs)
• src/astrom/pmAstrometryModel.c (modified) (9 diffs)
• src/astrom/pmAstrometryObjects.c (modified) (7 diffs)
• src/astrom/pmAstrometryObjects.h (modified) (1 diff)
• src/astrom/pmAstrometryRegions.h (modified) (1 diff)
• src/astrom/pmAstrometryUtils.c (modified) (5 diffs)
• src/astrom/pmAstrometryUtils.h (modified) (1 diff)
• src/astrom/pmAstrometryWCS.c (modified) (14 diffs)
• src/camera/pmFPAConstruct.c (modified) (1 diff)
• src/camera/pmFPAfile.c (modified) (2 diffs)
• src/camera/pmFPAfile.h (modified) (1 diff)
• src/camera/pmFPAfileFitsIO.c (modified) (1 diff)
• src/camera/pmFPAfileIO.c (modified) (10 diffs)
• src/camera/pmFPAfileIO.h (modified) (1 diff)
• src/detrend/pmDetrendDB.c (modified) (1 diff)
• src/detrend/pmDetrendDB.h (modified) (1 diff)
• src/detrend/pmPattern.c (modified) (10 diffs)
• src/detrend/pmPatternIO.c (modified) (2 diffs)
• src/detrend/pmPatternIO.h (modified) (1 diff)
• src/imcombine/pmStack.c (modified) (5 diffs)
• src/imcombine/pmStack.h (modified) (1 diff)
• src/objects/pmSourceIO_CMF.c.in (modified) (32 diffs)
• src/objects/pmSourceIO_Ghosts.c (modified) (6 diffs)

trunk/psModules/src/astrom/pmAstrometryDistortion.c

@@ -29 +29 @@
 #include "pmAstrometryRegions.h"
 #include "pmAstrometryDistortion.h"
+#include "pmAstrometryUtils.h"
 #include "pmKapaPlots.h"
 
@@ -183 +184 @@
     PS_ASSERT_PTR_NON_NULL(fpa, false);
     PS_ASSERT_ARRAY_NON_NULL(gradients, false);
+
+    int ExtraOrders = pmAstrometryGetExtraOrders();
 
     psPolynomial2D *localX = NULL;
@@ -340 +343 @@
     psRegion *region = pmAstromFPAExtent (fpa);
 
-    // XXX psFree (fpa->fromTPA);
-    // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(fpa->toTPA->x->nX+4, fpa->toTPA->x->nY+4);
-    // XXX fpa->fromTPA = psPlaneTransformInvert(myPT, fpa->toTPA, *region, 50);
-    // XXX psFree (myPT);
-
     // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear)
     psFree (fpa->fromTPA);
-    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 100, 6);
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, *region, 50, ExtraOrders);
     psFree (region);
 

trunk/psModules/src/astrom/pmAstrometryModel.c

@@ -461 +461 @@
     bool status;
 
+    int ExtraOrders = pmAstrometryGetExtraOrders();
+
     // set FITS cursor
     if (!psFitsMoveExtName (file->fits, "CHIPS")) {
@@ -498 +500 @@
         int nY = psMetadataLookupS32(&status, row, "NYORDER"); REQUIRE (status, "missing NYORDER");
         if (chip->toFPA == NULL) {
-            chip->toFPA = psPlaneTransformAlloc(nX, nY);
+            chip->toFPA = psPlaneTransformAlloc(nX, nY, PS_POLYNOMIAL_ORD); // chip->fpa uses ordinary poly
         } else {
             REQUIRE (chip->toFPA->x->nX == nX, "mismatch in chip order");
@@ -525 +527 @@
         // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear)
         psFree (chip->fromFPA);
-        chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6);
-
+        chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders);
         psFree (region);
     }
@@ -539 +540 @@
 
     bool status;
+
+    int ExtraOrders = pmAstrometryGetExtraOrders();
 
     if (!psFitsMoveExtName (file->fits, "FP")) {
@@ -565 +568 @@
         if (file->fpa->toTPA == NULL) {
             // allocate the new transformation
-            file->fpa->toTPA = psPlaneTransformAlloc(nX, nY);
+            file->fpa->toTPA = psPlaneTransformAlloc(nX, nY, PS_POLYNOMIAL_ORD); // fpa->tpa uses ORD
         } else {
             REQUIRE (file->fpa->toTPA->x->nX == nX, "mismatch in chip order");
@@ -585 +588 @@
     psRegion *region = pmAstromFPAExtent (file->fpa);
 
-    // XXX psFree (file->fpa->fromTPA);
-    // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(file->fpa->toTPA->x->nX+4, file->fpa->toTPA->x->nY+4);
-    // XXX file->fpa->fromTPA = psPlaneTransformInvert(myPT, file->fpa->toTPA, *region, 50);
-    // XXX psFree (myPT);
-
     // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear)
     psFree (file->fpa->fromTPA);
-    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 100, 6);
+    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50, ExtraOrders);
 
     psFree (model);
@@ -689 +687 @@
     bool status;
 
+    int ExtraOrders = pmAstrometryGetExtraOrders();
+
     // these externally supplied values are used to set the final transformation terms
     double RA  = psMetadataLookupF64 (&status, concepts, "FPA.RA"); REQUIRE (status, "missing FPA.RA");
@@ -728 +728 @@
     psLogMsg ("psModules.astrom", 4, "Position Angle: %f, Model Position Angle Zero Point: %f\n", POS, PosZero);
 
-//    psPlaneTransform *fromTPA = psPlaneTransformRotate (NULL, file->fpa->fromTPA, rotatorParity * (PosZero - POS));
-//    psFree (file->fpa->fromTPA);
-//    file->fpa->fromTPA = fromTPA;
-
     psPlaneTransform *toTPA = psPlaneTransformRotate (NULL, file->fpa->toTPA, rotatorParity * (POS - PosZero));
     psFree (file->fpa->toTPA);
@@ -739 +735 @@
     psRegion *region = pmAstromFPAExtent (file->fpa);
 
-    // XXX psFree (file->fpa->fromTPA);
-    // XXX psPlaneTransform *myPT = psPlaneTransformAlloc(file->fpa->toTPA->x->nX+4, file->fpa->toTPA->x->nY+4);
-    // XXX file->fpa->fromTPA = psPlaneTransformInvert(myPT, file->fpa->toTPA, *region, 50);
-    // XXX psFree (myPT);
-
     psFree (file->fpa->fromTPA);
-    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region,100, 10);
+    file->fpa->fromTPA = psPlaneTransformInvert(NULL, file->fpa->toTPA, *region, 50, ExtraOrders);
 
     psFree (region);

trunk/psModules/src/astrom/pmAstrometryObjects.c

@@ -40 +40 @@
 // XXX this is defined in pmPSFtry.h, which makes no sense
 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value);
 
 #define PM_ASTROMETRYOBJECTS_DEBUG 1
@@ -339 +340 @@
     psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
 
+    // we measure the stdev of the median residual in NxN bins.
+    // use only valid (not NAN) measurements
+    psVector *xPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *xResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+
     for (int i = 0; i < match->n; i++) {
         if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
@@ -344 +352 @@
         pmAstromObj *rawStar = raw->data[pair->raw];
         if (!isfinite(rawStar->dMag)) continue;
+
+        bool isValid = true;
+        isValid = isValid & isfinite (x->data.F32[i]);
+        isValid = isValid & isfinite (y->data.F32[i]);
+        isValid = isValid & isfinite (xRes->data.F32[i]);
+        isValid = isValid & isfinite (yRes->data.F32[i]);
+        if (isValid) {
+          psVectorAppend (xPosValid, x->data.F32[i]);
+          psVectorAppend (yPosValid, y->data.F32[i]);
+          psVectorAppend (xResValid, xRes->data.F32[i]);
+          psVectorAppend (yResValid, yRes->data.F32[i]);
+        }
+
+        // for the systematic error, use only high S/N stars
         if (rawStar->dMag > 0.02) continue;
 
@@ -369 +391 @@
       results->dXsys   = results->dYsys   = NAN;
       results->dXrange = results->dYrange = NAN;
+      results->dXstdev = results->dYstdev = NAN;
     } else {
       results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05);
@@ -375 +398 @@
       results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev);
       results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
-    }
-
-    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange);
-    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange);
+
+      results->dXstdev = pmAstrom2DSystematics (xPosValid, yPosValid, xResValid);
+      results->dYstdev = pmAstrom2DSystematics (xPosValid, yPosValid, yResValid);
+    }
+
+    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f, dXstdev: %f\n", results->dXsys, results->dXrange, results->dXstdev);
+    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f, dYstdev: %f\n", results->dYsys, results->dYrange, results->dYstdev);
 
     psFree (xErr);
@@ -386 +412 @@
     psFree (xResGood);
     psFree (yResGood);
+    psFree (xPosValid);
+    psFree (yPosValid);
+    psFree (xResValid);
+    psFree (yResValid);
 
     psFree (x);
@@ -395 +425 @@
 
     return (results);
+}
+
+# define VAL_COUNT 9.0
+# define MIN_COUNT 5.0
+
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value) {
+
+    // pre-filter the values to ensure no NANs, other invalid
+    if (xPos->n < VAL_COUNT*2*2) {
+      return NAN;
+    }
+
+    // generate a grid covering the full range of x,y and measure the median value in each
+    // grid cell.  calculate the stdev of those median values.
+
+    // VAL_COUNT (9) is the (min) goal density, but a single cell may have only MIN_COUNT (7)
+    // we have N points.  require a min of 9 pts per cell (configurable?).  grid is square.
+    // Ncell*Ncell*9 = Npts, Ncell = MIN(sqrt(Npts/9), 5)
+
+    int Ncell = PS_MIN(sqrt((float)(xPos->n / VAL_COUNT)), 2);  // have to at least have a 2x2 grid
+
+    // find the range of x,y values
+    float xMin = 1e9, xMax = -1e9, yMin = 1e9, yMax = 1e9;
+    for (int i = 0; i < xPos->n; i++) {
+        xMin = PS_MIN(xPos->data.F32[i],xMin);
+        xMax = PS_MAX(xPos->data.F32[i],xMax);
+        yMin = PS_MIN(yPos->data.F32[i],yMin);
+        yMax = PS_MAX(yPos->data.F32[i],yMax);
+    }
+
+    float xStep = Ncell / (xMax - xMin);
+    float yStep = Ncell / (yMax - yMin);
+  
+    if (isnan(xStep)) {
+      return NAN;
+    }
+    if (isnan(yStep)) {
+      return NAN;
+    }
+
+    psArray *xBin = psArrayAlloc (Ncell);
+    for (int ix = 0; ix < Ncell; ix++) {
+        psArray *yBin = psArrayAlloc (Ncell); 
+        xBin->data[ix] = yBin;
+        for (int iy = 0; iy < Ncell; iy++) {
+            yBin->data[iy] = psVectorAllocEmpty(128, PS_TYPE_F32);
+        }
+    }    
+
+    // xValue = ix/xStep + xMin -> ix = (xValue - xMin) * xStep;
+
+    for (int i = 0; i < xPos->n; i++) {
+        int ix = PS_MIN(PS_MAX(0, (xPos->data.F32[i] - xMin) * xStep), Ncell - 1);
+        int iy = PS_MIN(PS_MAX(0, (yPos->data.F32[i] - yMin) * yStep), Ncell - 1);
+    
+        psArray *yBin = xBin->data[ix];
+        psVector *Bin = yBin->data[iy];
+
+        psVectorAppend(Bin, value->data.F32[i]);
+    }
+
+    psVector *sample = psVectorAllocEmpty (Ncell*Ncell, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    // calculate the median for each vector and save on a vector
+    for (int ix = 0; ix < Ncell; ix++) {
+        psArray *yBin = xBin->data[ix];
+        for (int iy = 0; iy < Ncell; iy++) {
+            psVector *Bin = yBin->data[iy];
+            if (Bin->n < MIN_COUNT) { continue; }
+
+            // psVectorStats resets stats so we can call this repeatedly
+            psVectorStats (stats, Bin, NULL, NULL, 0); 
+            if (isfinite(stats->sampleMedian)) {
+                psVectorAppend (sample, stats->sampleMedian);
+            }
+        }
+    }    
+  
+    stats->options = PS_STAT_SAMPLE_STDEV;
+    psVectorStats (stats, sample, NULL, NULL, 0); 
+    float result = stats->sampleStdev;
+
+    if (!isfinite(stats->sampleStdev)) {
+      fprintf (stderr, "*** bad solution ***\n");
+    }
+
+    // NOTE: the elements of xBin are freed automatically, which extends down to the vectors
+    psFree (xBin); 
+    psFree (stats);
+    psFree (sample);
+  
+    return result;
 }
 

trunk/psModules/src/astrom/pmAstrometryObjects.h

@@ -102 +102 @@
     double  dXrange;                    ///< 10% - 90% range X residuals (unmasked, high S/N)
     double  dYrange;                    ///< 10% - 90% range Y residuals (unmasked, high S/N)
+    double  dXstdev;                    ///< stdev of median residual in NxN bins
+    double  dYstdev;                    ///< stdev of median residual in NxN bins
 }
 pmAstromFitResults;

trunk/psModules/src/astrom/pmAstrometryRegions.h

@@ -10 +10 @@
 #ifndef PM_ASTROMETRY_REGIONS_H
 #define PM_ASTROMETRY_REGIONS_H
+
+// used by pmAstrometryDistortion.c, pmAstrometryWCS.c, pmAstrometryModel.c
+// # define EXTRA_ORDERS 6
 
 /// @addtogroup Astrometry

trunk/psModules/src/astrom/pmAstrometryUtils.c

@@ -23 +23 @@
 #include "pmAstrometryUtils.h"
 
+static int transform_extra_orders = 3;
+
+int pmAstrometryGetExtraOrders (void) {
+  return transform_extra_orders;
+}
+
+int pmAstrometrySetExtraOrders (int orders) {
+  transform_extra_orders = orders;
+  return transform_extra_orders;
+}
+
 // this is used by the test output block
 static int Nout = 0;
@@ -120 +131 @@
 
 // convert a transformation L(x,y) to L'(x-xo,y-yo)
+// is this used for an upward (e.g., chip->fpa) or a downward (fpa->chip) transform?
 psPlaneTransform *psPlaneTransformSetCenter (psPlaneTransform *output, psPlaneTransform *input, double Xo, double Yo)
 {
 
-    // validate fit order
+    // validate fit type:
+    // polynomial in input transforms must match type -- and for now be ORD
+    psAssert (input->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (input->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
 
     if (output == NULL) {
-        output = psPlaneTransformAlloc(input->x->nX, input->x->nY);
+        output = psPlaneTransformAlloc(input->x->nX, input->x->nY, PS_POLYNOMIAL_ORD);
     }
 
@@ -169 +184 @@
 
 // rotate a transformation L(x,y) by theta
-psPlaneTransform *psPlaneTransformRotate (psPlaneTransform *output, psPlaneTransform *input, double theta)
+psPlaneTransform *psPlaneTransformRotate  (psPlaneTransform *output, psPlaneTransform *input, double theta)
 {
     /* given the polynomial transformations:
@@ -181 +196 @@
 
     if (output == NULL) {
-        output = psPlaneTransformAlloc(input->x->nX, input->x->nY);
+        // generate a new transform using the same order and type as the input
+        output = psPlaneTransformAlloc(input->x->nX, input->x->nY, input->x->type);
     }
 
@@ -216 +232 @@
 
     // all coeffs and masks initially set to 0
-    transform = psPlaneTransformAlloc (order, order);
+    transform = psPlaneTransformAlloc (order, order, PS_POLYNOMIAL_ORD);
 
     for (int i = 0; i <= order; i++) {

trunk/psModules/src/astrom/pmAstrometryUtils.h

@@ -25 +25 @@
 bool psPlaneDistortIsDiagonal (psPlaneDistort *distort);
 
+int pmAstrometryGetExtraOrders (void);
+int pmAstrometrySetExtraOrders (int orders);
+
 /// @}
 #endif

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -501 +501 @@
     psPlaneTransform *toFPA;
 
+    int ExtraOrders = pmAstrometryGetExtraOrders();
+
     // create transformation with 0,0 reference pixel and units of degrees/pixel
     toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
@@ -552 +554 @@
         // apply the exiting fromTPA transformation to make the new toFPA consistent with the toTPA layter
         // XXX this only works if toTPA is at most a linear transformation
-        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY);
+        psPlaneTransform *toFPAnew = psPlaneTransformAlloc(toFPA->x->nX, toFPA->x->nY, PS_POLYNOMIAL_ORD);
         for (int i = 0; i <= toFPA->x->nX; i++) {
             for (int j = 0; j <= toFPA->x->nY; j++) {
@@ -608 +610 @@
     // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear)
     psFree (chip->fromFPA);
-    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6);
+    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders);
     psFree (region);
 
@@ -648 +650 @@
      */
 
+    int ExtraOrders = pmAstrometryGetExtraOrders();
+
     psFree (chip->toFPA);
     if ((fabs(wcs->crpix1) > 0.01) || (fabs(wcs->crpix2) > 0.01)) {
       chip->toFPA = psPlaneTransformSetCenter (NULL, wcs->trans, -wcs->crpix1, -wcs->crpix2);
     } else {
-      chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY);
+      chip->toFPA = psPlaneTransformAlloc(wcs->trans->x->nX, wcs->trans->x->nY, PS_POLYNOMIAL_ORD);
 
       // copy the toFPA x,y, transformations to the wcs version
@@ -668 +672 @@
     // as of r40806, psPlaneTransformInvert supplies the extra order (if non-linear)
     psFree (chip->fromFPA);
-    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 100, 6);
+    chip->fromFPA = psPlaneTransformInvert(NULL, chip->toFPA, *region, 50, ExtraOrders);
     psFree (region);
 
@@ -700 +704 @@
     // the region defines the FPA pixels covered by the tranformation
     psFree (fpa->fromTPA);
-    int additional_orders = 4;  // This is the number of orders that should be added.
+    int additional_orders = pmAstrometryGetExtraOrders();  // This is the number of orders that should be added.
     bool status = false;
     int config_additional_orders = psMetadataLookupS32(&status,fpa->analysis, "ADDITIONAL_WCS_ORDERS");
@@ -706 +710 @@
       additional_orders = config_additional_orders; 
     }
-    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 100, additional_orders);
+    fpa->fromTPA = psPlaneTransformInvert(NULL, fpa->toTPA, region, 50, additional_orders);
     return true;
 }
@@ -739 +743 @@
     // XXX require fpa->toTPA->nX == 1
 
-    psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY);
+    psPlaneTransform *toTPA = psPlaneTransformAlloc(chip->toFPA->x->nX, chip->toFPA->x->nY, PS_POLYNOMIAL_ORD);
 
     for (int i = 0; i <= toTPA->x->nX; i++) {
@@ -958 +962 @@
     }
 
-    psPlaneTransform *newTrans = psPlaneTransformAlloc(1, 1);
+    psPlaneTransform *newTrans = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD);
 
     if (!psPlaneTransformFit(newTrans, src, dst, 0, 0)) {
@@ -997 +1001 @@
     PS_ASSERT_PTR_NON_NULL(inChip, NULL);
 
+    int ExtraOrders = pmAstrometryGetExtraOrders();
+
     if (outFPA == NULL) {
         outFPA = inFPA;
@@ -1032 +1038 @@
 
     // NOTE: the extraOrders value (4) should be ignored since outToFPA is specified to be linear
-    psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 100, 6);
+    psPlaneTransform *outFromFPA = psPlaneTransformInvert(NULL, outToFPA, *outputBounds, 50, ExtraOrders);
     if (!outFromFPA) {
         psFree(outToFPA);
@@ -1116 +1122 @@
     }
 
-    psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1);
+    psPlaneTransform *newToFPA = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD);
     newToFPA->x->coeffMask[1][1] = 1;
     newToFPA->y->coeffMask[1][1] = 1;
@@ -1148 +1154 @@
     psFree(dst);
 
-    // this is a linear transformation
+    // this is a linear transformation, no extra orders are needed
     psPlaneTransform *newFromFPA = psPlaneTransformInvert(NULL, newToFPA, *bounds, 1, 0);
     if (!newFromFPA) {
@@ -1187 +1193 @@
     psMemSetDeallocator(wcs, (psFreeFunc) pmAstromWCSFree);
 
-    wcs->trans = psPlaneTransformAlloc (nXorder, nYorder);
+    // note: WCS transforms are always defined as chip to sky
+    wcs->trans = psPlaneTransformAlloc (nXorder, nYorder, PS_POLYNOMIAL_ORD);
     wcs->toSky = NULL;
     wcs->wcsCDkeys = 0;

trunk/psModules/src/camera/pmFPAConstruct.c

@@ -853 +853 @@
     } else {
         const char *content = getContent(fileInfo, phdu->header, contents); // The chip type
+        if (!content) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to find CONTENT entry in header");
+            return false;
+        }
 
         int chipNum = -1;               // Chip number

trunk/psModules/src/camera/pmFPAfile.c

@@ -551 +551 @@
     if (!strcasecmp(type, "KH.CORRECT"))     {
         return PM_FPA_FILE_KH_CORRECT;
+    }
+    if (!strcasecmp(type, "PATTERN.ROW.AMP"))     {
+        return PM_FPA_FILE_PATTERN_ROW_AMP;
     }
     if (!strcasecmp(type, "SUBKERNEL"))     {
@@ -606 +609 @@
       case PM_FPA_FILE_KH_CORRECT:
         return ("KH.CORRECT");
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
+        return ("PATTERN.ROW.AMP");
       case PM_FPA_FILE_SUBKERNEL:
         return ("SUBKERNEL");

trunk/psModules/src/camera/pmFPAfile.h

@@ -52 +52 @@
     PM_FPA_FILE_SRCTEXT,
     PM_FPA_FILE_PATTERN,
+    PM_FPA_FILE_PATTERN_ROW_AMP,
     PM_FPA_FILE_LINEARITY,
     PM_FPA_FILE_EXPNUM,

trunk/psModules/src/camera/pmFPAfileFitsIO.c

@@ -152 +152 @@
       case PM_FPA_FILE_ASTROM_REFSTARS: 
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
         {
           pmHDU *hdu = pmFPAviewThisHDU(view, fpa);

trunk/psModules/src/camera/pmFPAfileIO.c

@@ -145 +145 @@
 }
 
+// list all defined pmFPAfiles
+bool pmFPAfileIOList (pmConfig *config)
+{
+    PS_ASSERT_PTR_NON_NULL(config, false);
+    PS_ASSERT_PTR_NON_NULL(config->files, false);
+
+    psMetadata *files = config->files;
+
+    // attempt to perform all create, read, write, close operations
+    psMetadataItem *item = NULL;
+    psMetadataIterator *iter = psMetadataIteratorAlloc (files, PS_LIST_HEAD, NULL);
+    while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
+        pmFPAfile *file = item->data.V;
+        
+        fprintf (stderr, "%s : %d %d %d : %d %d %d %d\n", file->name, file->type, file->mode, file->state,
+                 file->fileLevel, file->dataLevel, file->freeLevel, file->mosaicLevel);
+    }
+    psFree (iter);
+    return true;
+}
+
 // read the file, if necessary and possible
 bool pmFPAfileRead(pmFPAfile *file, const pmFPAview *view, pmConfig *config)
@@ -231 +252 @@
       case PM_FPA_FILE_PATTERN:
         status = pmPatternRead(view, file, config);
+        break;
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
+        status = pmPatternRowAmpRead(view, file, config);
         break;
       case PM_FPA_FILE_SX:
@@ -342 +366 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
       case PM_FPA_FILE_JPEG:
       case PM_FPA_FILE_KAPA:
@@ -359 +384 @@
 
     if (file->state & PM_FPA_STATE_INACTIVE) {
-        psTrace("psModules.camera", 6, "skip write for %s, files is inactive", file->name);
+        psTrace("psModules.camera", 6, "skip write for %s, file is inactive", file->name);
         return true;
     }
@@ -430 +455 @@
       return true;
     }
+    if (file->type == PM_FPA_FILE_PATTERN_ROW_AMP) {
+      psTrace("psModules.camera", 6, "skip write for %s, no write function defined", file->name);
+      return true;
+    }
 
     // open the file if not yet opened
@@ -526 +555 @@
       case PM_FPA_FILE_KH_CORRECT:
         psError(PS_ERR_IO, true, "cannot write type KH.CORRECT (%s)", file->name);
+        break;
+
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
+        psError(PS_ERR_IO, true, "cannot write type PATTERN.ROW.AMP (%s)", file->name);
         break;
 
@@ -603 +636 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
       case PM_FPA_FILE_LINEARITY:
       case PM_FPA_FILE_EXPNUM:
@@ -681 +715 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
       case PM_FPA_FILE_EXPNUM:
         psTrace ("psModules.camera", 6, "NOT freeing %s (%s) : save for further analysis\n", file->filename, file->name);
@@ -844 +879 @@
       case PM_FPA_FILE_ASTROM_REFSTARS:
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
       case PM_FPA_FILE_LINEARITY:
       case PM_FPA_FILE_EXPNUM:
@@ -1048 +1084 @@
       case PM_FPA_FILE_ASTROM_MODEL:
       case PM_FPA_FILE_KH_CORRECT:
+      case PM_FPA_FILE_PATTERN_ROW_AMP:
       case PM_FPA_FILE_SX:
       case PM_FPA_FILE_RAW:

trunk/psModules/src/camera/pmFPAfileIO.h

@@ -55 +55 @@
 bool pmFPAfileReadPHU (pmFPAfile *file, const pmFPAview *view, pmConfig *config);
 
+bool pmFPAfileIOList (pmConfig *config);
+
 /// @}
 # endif

trunk/psModules/src/detrend/pmDetrendDB.c

@@ -112 +112 @@
         DETREND_STRING_CASE(AUXMASK);
         DETREND_STRING_CASE(KH_CORRECT);
+        DETREND_STRING_CASE(PATTERN_ROW_AMP);
     default:
         return NULL;

trunk/psModules/src/detrend/pmDetrendDB.h

@@ -41 +41 @@
     PM_DETREND_TYPE_AUXMASK,
     PM_DETREND_TYPE_KH_CORRECT,
+    PM_DETREND_TYPE_PATTERN_ROW_AMP,
 } pmDetrendType;
 

trunk/psModules/src/detrend/pmPattern.c

@@ -9 +9 @@
 #define PATTERN_ROW_BKG_FIX 1
 
+/* some in-line notes:
+
+   patternMaskRow sets the data value to NAN, inconsistent with new plan
+   
+   here is the outline of pmPatternRow
+
+   * at this point, we have already done overscan subtraction, right?
+
+   * measure stats on the full cell (MEDIAN, STDEV)
+   ** subsample?
+   ** if it fails, it masks the entire cell and set the value to NAN
+
+   * calculate an upper and lower threshold (median +/- T * sigma)
+   * define a normalized x-coordinate ('index') : 
+   ** see note below on chebys
+
+   * each row is treated independently
+   * pixels are masked for the fit if they are out-of-range 
+     or if they are already masked
+
+   ** note that the clipping threshold will be larger if there 
+      are pixels which have astronomical structures
+      
+      a possible better option would be to set the threshold based on the median
+      and a sigma calculated from Poisson stats (do we know the gain?)
+      
+   ** fit is allowed to proceed if even N+1 pixels exist, which is clearly too low
+   
+   ** Remaining pixels are fitted with clip-fit 
+
+   ** solution is subtracted from the data
+   (this is implemented with psPolynomial1DEvalVector)
+   perhaps faster if we fixed the order to 2 and hardwired the result
+   
+   * after each row is fitted, the intercept (A value) is fitted
+   as a function of the y-coordinate and the result is subtracted
+
+   * the slope value is also fitted as a function of the column 
+     and added back in -- I'm not sure I understand this step.
+
+   *****************
+
+   ** what we calculate are related to chebychevs (domain is -1 : +1)
+   *** T0(x) = 1
+   *** T1(x) = x
+   *** T2(x) = 2x^2 - 1
+
+   *** we calculate y = A + Bx + Cx^2
+
+   a_0 + a_1 x + a_2 (2x^2 - 1) = A + B + Cx^2
+
+   a_1       = B
+   a_0 - a_2 = A
+   2 a_2     = C
+
+   a_0       = A + C/2
+   a_1       = B
+   a_2       = C/2
+
+   *****************
+   
+   I have 3 goals in re-working the code:
+   
+   1) improve overall speed
+   2) improve reliability of the fit
+   3) skip fit if we can
+
+   Let's assume the signal in the cell is light + bias drift
+
+   The bias drift has an amplitude of ~5 - 10 DN
+
+   That makes a detectable source with ~N * a few counts (multiple pixels in a row)
+   
+   So, the effective flux is ~10 * 5 = 50 DN
+   for which sky level is this value - 3 sigma?
+   
+   50 / sqrt(sky sigma^2 * effective area)
+
+   area ~ 5pixels, sky sigma^2 = sky
+
+   10 * Npix / sqrt(sky * Npix) = 3
+
+   Npix = sky * (S/N)^2 / (peak^2)
+   
+   sky = Npix * peak^2 / SN^2
+
+   if (sky < Npix * peak^2 / SN^2), we should skip:
+   
+   Npix ~ 5
+   peak ~ 10
+   SN ~ 3 (or even less)
+
+   sky < 5 * 100 / 9 = 55 or so
+
+
+   To address these in order:
+   
+   1) speed: 
+   * the analysis threaded is not threaded: thread across cells
+   
+   2) 
+
+ */
 
 // Mask a row as bad
@@ -33 +136 @@
 }
 
+// Comparison and swap functions for sorting values directly
+#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B])
+#define SORT_SWAP(TYPE,A,B) {                   \
+    if (A != B) { \
+        TYPE temp = sampleArray[A];                     \
+        sampleArray[A] = sampleArray[B]; \
+        sampleArray[B] = temp; \
+    } \
+}
+
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // Measurement and application
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                          psStatsOptions clipMean, psStatsOptions clipStdev,
+                          psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                        psStatsOptions clipMean, psStatsOptions clipStdev,
+                        psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+// XXX allow user choice of binned vs unbinned analysis?
 bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
+                  psStatsOptions clipMean, psStatsOptions clipStdev,
+                  psImageMaskType maskVal, psImageMaskType maskBad) {
+
+  bool status = false;
+  if (true) {
+    status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+  } else {
+    status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+  }
+  return status;
+}
+
+// USE_BACKGROUND_STDEV: if TRUE, the analysis will use the measured robust stdev to clip the out-of-range pixles
+// if FALSE, the stdev will be estimated based on the Poisson statistics of the median (sky level).
+# define USE_BACKGROUND_STDEV 0
+
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
                   psStatsOptions clipMean, psStatsOptions clipStdev,
                   psImageMaskType maskVal, psImageMaskType maskBad)
@@ -48 +189 @@
     PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
 
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
     psImage *image = ro->image;         // Image to correct
     psImage *mask = ro->mask;           // Mask for image
@@ -55 +203 @@
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
     if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
-        psWarning("Unable to calculate statistics on readout; masking entire readout.");
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
         psErrorClear();
         psFree(stats);
         psFree(rng);
-        psImageInit(image, NAN);
-        if (mask) {
-            psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
-        }
-        if (ro->variance) {
-            psImageInit(image, NAN);
-        }
         return true;
     }
+
+# if (USE_BACKGROUND_STDEV) 
     float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
     float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
     float background = stats->robustMedian;
+# else
+    // the signal we are looking for is a small variation on top of the background.  if
+    // the background is uniform with only read noise + sky noise, then the pixel-to-pixel
+    // stdev should only be due to known noise sources and predictable.  If the
+    // pixel-to-pixel variations are from other features, then those variations will
+    // probably dominate the row-by-row bias variations.
+
+    // instead of using the image pixel statistics to measure the stdev, lets assume only
+    // dark noise plus poisson sky noise.  we are not carrying in the read noise, but it is
+    // fairly modest for GPC1 (~10 DN)
+
+    // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have
+    // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17
+
+    // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2)
+    // so smaller by a factor of 1.4 than what we predict, which is not very large
+
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
+    // XXX retrieve noise and gain from 'concepts' if possible
+# define READNOISE 10 /* arbitrary number */
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); 
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
+    float background = stats->robustMedian;
+
+    // if the noise from the background is too large 
+    float significance = nominalAmplitude / sigma;
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      return true; 
+    }
+    fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+# endif    
+    
     psFree(stats);
     psFree(rng);
@@ -109 +294 @@
     psVectorInit(yaxisMask, 0);
 #endif
+
+    // we really need more than order + 1 points (= 4).
+    // this should be tunable, but let's try 5 - 10%
+    int validNmin = numCols * 0.1;
+
     for (int y = 0; y < numRows; y++) {
         psVectorInit(clipMask, 0);
         data = psImageRow(data, image, y);
         int num = 0;                    // Number of good pixels
+
+        // if the unmasked pixels only span a small range in x then we cannot fit the
+        // 2nd order polynomial variations very well.  Require a minimum fractional range
+        float validXmin = +1;
+        float validXmax = -1;
+
+        // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY)
+        // (x % 3) ||
         for (int x = 0; x < numCols; x++) {
-            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
-                data->data.F32[x] < lower || data->data.F32[x] > upper) {
-                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                data->data.F32[x] < lower || data->data.F32[x] > upper) {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
             } else {
                 clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0;
                 num++;
+                validXmin = PS_MIN(indices->data.F32[x], validXmin);
+                validXmax = PS_MAX(indices->data.F32[x], validXmax);
             }
         }
-        if (num < order + 1) {
+
+        // XXX how much time is spent in the fitting
+        if (num < validNmin) {
             // Not enough points to fit
             patternMaskRow(ro, y, maskBad);
@@ -131 +333 @@
             continue;
         }
+        // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise?
         if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) {
             psWarning("Unable to fit polynomial to row %d", y);
@@ -215 +418 @@
           for (int x = 0; x < numCols; x++) {
             image->data.F32[y][x] += solution->data.F32[y];
-            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]);
+            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]);
           }
           corr->data.F64[y][0]  -= solution->data.F32[y];
@@ -241 +444 @@
           for (int x = 0; x < numCols; x++) {
             image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
-            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]);
+            // XXX EAM : this was [x] which is wrong
+            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]);
           }
           corr->data.F64[y][1]  -= solution->data.F32[y] ;
@@ -262 +466 @@
 
     psFree(indices);
+    psFree(clip);
+    psFree(clipMask);
+    psFree(poly);
+    psFree(data);
+
+    return true;
+}
+
+# define NPIX 15
+
+// bin by NPIX in the x-direction to reduce the number of calculations needed to measure
+// the pattern
+bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                  psStatsOptions clipMean, psStatsOptions clipStdev,
+                  psImageMaskType maskVal, psImageMaskType maskBad)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+    PS_ASSERT_INT_NONNEGATIVE(order, false);
+    PS_ASSERT_INT_NONNEGATIVE(iter, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
+
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
+    psImage *image = ro->image;         // Image to correct
+    psImage *mask = ro->mask;           // Mask for image
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
+        psErrorClear();
+        psFree(stats);
+        psFree(rng);
+        
+        // EAM 20211011 : we used to mask cells which fail the above, but this seems excessive
+        // psImageInit(image, NAN);
+        // if (mask) {
+        //     psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
+        // }
+        // if (ro->variance) {
+        //     psImageInit(image, NAN);
+        // }
+        return true;
+    }
+
+    // if USE_BACKGROUND_STDEV is TRUE, the observed standard deviation is used to set the
+    // thresholds.  this is going to be an overestimate if there is any structure in the
+    // image.  If FALSE, the thresholds are set based on poisson stats for the background
+    // level.  We assume the gain is 1, so this is an overestimate if the gain is > 1
+
+# if (USE_BACKGROUND_STDEV) 
+    float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
+    float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
+    float background = stats->robustMedian;
+# else
+    // the signal we are looking for is a small variation on top of the background.  if
+    // the background is uniform with only read noise + sky noise, then the pixel-to-pixel
+    // stdev should only be due to known noise sources and predictable.  If the
+    // pixel-to-pixel variations are from other features, then those variations will
+    // probably dominate the row-by-row bias variations.
+
+    // instead of using the image pixel statistics to measure the stdev, lets assume only
+    // dark noise plus poisson sky noise.  we are not carrying in the read noise, but it is
+    // fairly modest for GPC1 (~10 DN)
+
+    // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have
+    // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17
+
+    // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2)
+    // so smaller by a factor of 1.4 than what we predict, which is not very large
+
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
+    // retrieve noise and gain from 'concepts' if possible
+# define READNOISE 10 /* arbitrary number */
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); 
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
+    float background = stats->robustMedian;
+
+    // if the noise from the background is too large 
+    float significance = nominalAmplitude / sigma;
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      psFree(stats);
+      psFree(rng);
+      return true; 
+    }
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+# endif    
+
+    psFree(stats);
+    psFree(rng);
+
+    // the vector 'indices' maps the x-coordinate to a range [-1:1].  the element number (i) of indices
+    // related to the x-coordinate (column number) by x = (i + 0.5) * NPIX
+
+    int nSamples = numCols / NPIX;
+
+    psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fit solutions
+    psVector *xFit = psVectorAlloc(nSamples, PS_TYPE_F32); // x-coordinate for fitting
+    psVector *yFit = psVectorAlloc(nSamples, PS_TYPE_F32); // flux values for fitting
+
+    // xFit elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX
+
+    float norm = 2.0 / (float)numCols;  // Normalisation for indices
+    for (int sample = 0; sample < nSamples; sample ++) {
+        int x = (sample + 0.5)*NPIX;
+        xFit->data.F32[sample] = x * norm - 1.0;
+    }
+    for (int x = 0; x < numCols; x ++) {
+        indices->data.F32[x] = x * norm - 1.0;
+    }
+
+    psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics
+    clip->clipIter = iter;
+    clip->clipSigma = rej;
+    psVector *clipMask = psVectorAlloc(nSamples, PS_TYPE_VECTOR_MASK); // Mask for clipping
+    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit
+    psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit
+
+    psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied
+    psImageInit(corr, NAN);
+
+    // CZW: 2011-11-30
+    // Define the vectors to hold the "x" and "y" slope trends.
+    // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term
+    // of the PATTERN.ROW fit between individual rows across the cell.  Similarly, the 1-st
+    // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we
+    // are fitting with PATTERN.ROW in the first place).  However, the thing we're trying to
+    // fix with PATTERN.ROW is the detector level bias wiggles.  These should be overlaid on
+    // the true sky level.  Therefore, simply applying the PATTERN.ROW correction will
+    // introduce cell-to-cell sky variations as these two trends are removed.  To avoid this,
+    // We store the 0th and 1st order values used for each row, and then fit a polynomial to
+    // these results.  By re-adding these systematic trends back, we can remove the row-to-row
+    // variations without improperly removing the real sky trend.
+    psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term
+    psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit
+    psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term
+    psVectorInit(yaxisMask, 0);
+
+    // validNmin is the minimum number of samples needed to measure the trend.
+    // this should be tunable, but let's try 5 - 10%
+    int validNmin = PS_MAX (nSamples * 0.1, order + 2);
+
+    for (int y = 0; y < numRows; y++) {
+        psVectorInit(clipMask, 0);
+        data = psImageRow(data, image, y);
+        int num = 0;                    // Number of good pixels
+
+        // if the unmasked pixels only span a small range in x then we cannot fit the
+        // 2nd order polynomial variations very well.  Require a minimum fractional range
+        float validXmin = +1;
+        float validXmax = -1;
+
+        for (int sample = 0; sample < nSamples; sample ++) {
+
+            // store valid samples in the array to be sorted
+            float sampleArray[NPIX];
+            int seq = 0;
+            for (int j = 0; j < NPIX; j++) {
+                int pix = sample  * NPIX + j; // real pixel elements in x-dir
+                psAssert (pix >= 0, "invalid pix value");
+                psAssert (pix < numCols, "invalid pix value");
+                if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][pix] & maskVal)) continue;
+                if (data->data.F32[pix] < lower || data->data.F32[pix] > upper) continue;
+                sampleArray[seq] = data->data.F32[pix]; // store the value to be sorted
+                seq ++;
+            } 
+            if (seq < 1) {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF;
+                yFit->data.F32[sample] = NAN;
+                continue;
+            }
+            // note that we are treating the x-coordinate as the center
+            // of the binned pixel group, even if some or most pixels have
+            // been masked.  compared to the amplitude of the slope, this
+            // error is small
+            clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0;
+            validXmin = PS_MIN(xFit->data.F32[sample], validXmin);
+            validXmax = PS_MAX(xFit->data.F32[sample], validXmax);
+            num++;
+
+            // PSSORT operates on sampleArray (see define of macro SORT_SWAP above)
+            PSSORT (seq, SORT_COMPARE, SORT_SWAP, float);
+
+            int midPt = 0.5 * seq;
+            if (seq % 2 == 1) { psAssert (midPt >= 0, "invalid midPt"); }
+            if (seq % 2 == 0) { psAssert (midPt >= 1, "invalid midPt"); }
+            psAssert (midPt < NPIX, "invalid midPt");
+
+            float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]);
+            yFit->data.F32[sample] = medValue;
+        }
+
+        // If not enough points are valid, skip
+        if (num < validNmin) {
+            // Ignore this row in our subsequent fits, because the fit failed.
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+        // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise?
+        if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, yFit, NULL, xFit)) {
+            psWarning("Unable to fit polynomial to row %d", y);
+            psErrorClear();
+            // Ignore this row in our subsequent fits, because the fit failed.
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+        // Store the results we found for this row.
+        yaxisData->data.F32[y] = poly->coeff[0];
+        xaxisData->data.F32[y] = poly->coeff[1];
+        psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]);
+
+        memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
+        if (!solution) {
+            psWarning("Unable to evaluate polynomial for row %d", y);
+            psErrorClear();
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+            continue;
+        }
+
+        psAssert (solution->n == numCols, "oops");
+        for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] -= solution->data.F32[x];
+            psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]);
+        }
+        psFree(solution);
+    }
+
+    // Put the global trends back that were removed by the PATTERN.ROW correction.
+    // Set up the indices for the polynomial
+    psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32);
+    norm = 2.0 / (float)numRows;
+    for (int y = 0; y < numRows; y++) {
+      yaxisIndices->data.F32[y] = y * norm - 1.0;
+      psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]);
+    }
+
+    // Fit the trend of the constant term, producing the y-axis global trend
+    psStatsInit(clip);
+    psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
+    if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) {
+      psWarning("Unable to fit polynomial to y-axis trend");
+      psErrorClear();
+      // If we've failed, we need to do something, so add back in the background level, and
+      // expect that the final image will have background mismatches.
+      for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+          image->data.F32[y][x] += background;
+        }
+        corr->data.F64[y][0]  -= background;
+      }
+    } else {
+      psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices);
+      if (!solution) {
+        psWarning("Unable to evaluate polynomial");
+        psErrorClear();
+        // If we've failed, we need to do something, so add back in the background level, and
+        // expect that the final image will have background mismatches.
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += background;
+          }
+          corr->data.F64[y][0]  -= background;
+        }
+      } else {
+        psAssert (solution->n == numRows, "oops");
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += solution->data.F32[y];
+            // XXX EAM : this was [x], which is wrong
+            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]);
+          }
+          corr->data.F64[y][0]  -= solution->data.F32[y];
+        }
+      }
+      psFree(solution);
+    }      
+
+    // Fit the trend of the linear term, producing the x-axis global trend
+    // We can use the same mask vector, as the same rows failed the row-fit earlier.
+    psStatsInit(clip);
+    psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit.
+    if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) {
+      psWarning("Unable to fit polynomial to x-axis trend");
+      psErrorClear();
+    } else {
+      psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices);
+      if (!solution) {
+        psWarning("Unable to evaluate polynomial");
+        psErrorClear();
+      } else {
+        psAssert (solution->n == numRows, "oops");
+        for (int y = 0; y < numRows; y++) {
+          for (int x = 0; x < numCols; x++) {
+            image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x];
+            // XXX EAM : this was set to [x] which is wrong (numCols > numRows)
+            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]);
+          }
+          corr->data.F64[y][1]  -= solution->data.F32[y] ;
+        }
+      }
+      psFree(solution);
+    }
+    psFree(yaxisPoly);
+    psFree(xaxisPoly);
+    psFree(yaxisIndices);
+    psFree(yaxisMask);
+    psFree(yaxisData);
+    psFree(xaxisData);
+    
+    psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
+                       "Pattern row correction", corr);
+    psFree(corr);
+
+    psFree(indices);
+    psFree(xFit);
+    psFree(yFit);
     psFree(clip);
     psFree(clipMask);
@@ -1316 +1855 @@
     return true;
 }
+

trunk/psModules/src/detrend/pmPatternIO.c

@@ -8 +8 @@
 #include "pmFPAview.h"
 #include "pmFPAfile.h"
+#include "pmFPAUtils.h"
 #include "pmFPAfileFitsIO.h"
 #include "pmFPAHeader.h"
@@ -452 +453 @@
     return pmReadoutReadPattern(readout, file->fits);
 }
+
+/**************** PatternRowAmp(litude) I/O *************************/
+
+bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
+    {
+        // read the full model in one pass: require the level to be FPA
+        if (view->chip != -1) {
+            psError(PS_ERR_IO, false, "Pattern Row Amplitude must be read at the FPA level");
+            return false;
+        }
+
+        if (!pmPatternRowAmpReadFPA (file)) {
+            psError(PS_ERR_IO, false, "Failed to read Pattern Row Amplitude for fpa");
+            return false;
+        }
+        return true;
+    }
+
+// read in all chip-level Pattern Row Amplitude data for this FPA
+bool pmPatternRowAmpReadFPA (pmFPAfile *file) {
+
+    if (!pmPatternRowAmpReadChips (file)) {
+        psError(PS_ERR_IO, false, "Failed to read Pattern Row Amplitude for chips");
+        return false;
+    }
+
+    return true;
+}
+
+// read the set of tables, one for each chip
+bool pmPatternRowAmpReadChips (pmFPAfile *file) {
+
+    bool haveData, status;
+
+    // loop over the extensions
+    // for each extension, use the extname (eg, XY01.ptn) to assign to a chip
+
+    // move to the start of the file
+    haveData = psFitsMoveExtNum (file->fits, 1, false);
+    if (!haveData) {
+        psError(PS_ERR_IO, false, "Failed to read even the first extension?");
+        return false;
+    }
+
+    int nGood = 0;
+    int nTotal = 0;
+    while (haveData) {
+
+        // load the header
+        psMetadata *header = psFitsReadHeader(NULL, file->fits); // The FITS header
+        if (!header) psAbort("cannot read model header");
+
+        // load the full model in one shot
+        psArray *model = psFitsReadTable (file->fits);
+        if (!model) psAbort("cannot read model");
+        
+        // determine the chip:
+        char *extname = psMetadataLookupStr (&status, header, "EXTNAME");
+        psLogMsg ("psModules.detrend", 8, "read %ld rows from Pattern Row Amplitude file, extname %s\n", model->n, extname);
+        nTotal += model->n;
+
+        // I expect to find a name of the form: chipName.ptn (eg, XY01.ptn)
+        // where chipName like 'XY01'
+        psAssert (strlen(extname) == 8, "invalid extension %s", extname);
+        psAssert (extname[5] == 'p', "invalid extension %s", extname);
+        psAssert (extname[6] == 't', "invalid extension %s", extname);
+        psAssert (extname[7] == 'n', "invalid extension %s", extname);
+
+        char chipName[5];
+        strncpy (chipName, extname, 4);
+        chipName[4] = 0;
+
+        pmChip *chip = pmConceptsChipFromName (file->fpa, chipName);
+        if (!chip) psAbort ("invalid chip?");
+
+        // parse the model entries
+        for (int i = 0; i < model->n; i++) {
+            psMetadata *row = model->data[i];
+            psAssert (row, "missing model row");
+
+            char *cellName = psMetadataLookupStr(&status, row, "CELL_NAME");
+            if (!cellName) continue;
+
+            float amplitude = psMetadataLookupF32(&status, row, "VALUE_MEDIAN");
+
+            int cellNumber = pmChipFindCell (chip, cellName);
+            psAssert ((cellNumber >=0) && (cellNumber < chip->cells->n), "invalid cell number");
+
+            pmCell *cell = chip->cells->data[cellNumber];
+            if (!cell) continue;
+
+            psAssert (cell->analysis, "oops");
+
+            psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
+            nGood ++;
+        }
+
+        psFree (model);
+        psFree (header);
+
+        // move to the next extension
+        haveData = psFitsMoveExtNum (file->fits, 1, true);
+    }
+    psLogMsg ("psModules.detrend", 4, "read %d of %d rows from Pattern Row Amplitude file\n", nGood, nTotal);
+
+    return true;
+}

trunk/psModules/src/detrend/pmPatternIO.h

@@ -35 +35 @@
     );
 
+
+bool pmPatternRowAmpRead (const pmFPAview *view, pmFPAfile *file, const pmConfig *config);
+bool pmPatternRowAmpReadFPA (pmFPAfile *file);
+bool pmPatternRowAmpReadChips (pmFPAfile *file);
+
 #endif

trunk/psModules/src/imcombine/pmStack.c

@@ -756 +756 @@
         // *goodMask &= mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn]; // save the mask bits still used
         // check for set bits and increment counter as appropriate
+        // count the number of times a given mask bit is set in the input pixels.
+        // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+        // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
         { 
             psImageMaskType value = 0x0001;
@@ -961 +964 @@
 
     // KMM values;
-    float Punimodal = 1.0,KMMmean = NAN,KMMsigma = NAN,KMMpi = NAN;
+    float Punimodal = 1.0, KMMmean = NAN, KMMsigma = NAN, KMMpi = NAN;
     int KMM_MINIMUM_INPUTS = 6;
     bool useKMM = false;
@@ -1510 +1513 @@
 }
 
+# define MIN_GOOD_PERCENTILE 2
+bool pmStackCombineByPercentile(
+    pmReadout *combined,
+    psArray *stackData,
+    psF64 minRange,
+    psF64 maxRange,
+    psImageMaskType badMaskBits,        // treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
+    psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+) {
+    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
+    int xSize, ySize;                   // Size of the output image
+
+    // we need to copy the readouts to their own array for validation
+    psArray *stackReadouts = psArrayAlloc(stackData->n);
+    for (int i = 0; i < stackData->n; i++) {
+        pmStackData *data = stackData->data[i]; // Stack data for this input
+        pmReadout *ro = data->readout;  // Readout of interest
+        stackReadouts->data[i] = psMemIncrRefCounter(ro);
+    }
+
+    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, stackReadouts)) {
+        psError(psErrorCodeLast(), false, "Input stack is not valid.");
+        psFree(stackReadouts);
+        return false;
+    }
+    psFree(stackReadouts);
+
+    psVector *pixelData = psVectorAlloc(stackData->n, PS_TYPE_F32);
+    psVector *pixelMask = psVectorAlloc(stackData->n, PS_TYPE_VECTOR_MASK);
+
+    for (int y = minInputRows; y < maxInputRows; y++) {
+        for (int x = minInputCols; x < maxInputCols; x++) {
+
+            int nGood = 0; 
+            for (int i = 0; i < stackData->n; i++) {
+
+                pmStackData *data = stackData->data[i]; // Stack data for this input
+                pmReadout *ro = data->readout;  // Readout of interest
+                psImage *image = ro->image;
+                psImage *mask = ro->mask;
+
+                int xIn = x - data->readout->col0;
+                int yIn = y - data->readout->row0; // Coordinates on input readout
+
+                if (!isfinite(image->data.F32[yIn][xIn])) continue;
+                if (fabs(image->data.F32[yIn][xIn]) > 1e5) continue;
+                // XXX need to test the input mask as well here
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & badMaskBits);
+
+                // XXX save the count of suspect bits
+
+# if (0)
+                // count the number of times a given mask bit is set in the input pixels.
+                // NOTE: since we have explicitly skipped the pixels with any bad bits, these are only
+                // the suspect bits (nGoodBits is a bit of a misnomer: it is more like 'nSuspectBitsForGoodInputs'
+                if (1) { 
+                  psImageMaskType value = 0x0001;
+                  for (int nbit = 0; nbit < 16; nbit ++) {
+                    if (mask->data.PS_TYPE_IMAGE_MASK_DATA[yIn][xIn] & value) {
+                      // XXX nGoodBits[nbit] ++;
+                    }
+                    value <<= 1;
+                  }
+                }
+# endif
+
+                pixelData->data.F32[nGood] = image->data.F32[yIn][xIn];
+                nGood ++;
+            }
+            pixelData->n = nGood;
+            psVectorSortInPlace (pixelData);
+
+            if (nGood < MIN_GOOD_PERCENTILE) {
+                combined->image->data.F32[y][x] = NAN;
+                combined->variance->data.F32[y][x] = NAN;
+                combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 1;
+                continue;
+            }
+
+# if (0)
+# define SUSPECT_FRACTION 0.65
+            // set the mask bits if nGoodBits[i] > f*numGood
+            if (0) {
+              psImageMaskType value = 0x0001;
+              for (int nbit = 0; nbit < 16; nbit ++) {
+                if (nGoodBits[nbit] > SUSPECT_FRACTION*numGood) {
+                  *goodMask |= value;
+                }
+                value <<= 1;
+              }
+            }
+# endif
+
+            int Ns = MIN(MAX(0, minRange * nGood),nGood); // e.g., 0.1 * 50 = 5
+            int Ne = MIN(MAX(0, maxRange * nGood),nGood); // e.g., 0.9 * 50 = 45 
+            int Npt = Ne - Ns;
+
+            float sum = 0.0;
+            for (int n = Ns; n < Ne; n++) {
+                sum += pixelData->data.F32[n];
+            }
+            float mean = sum / (float) Npt;
+
+            float varSum = 0.0;
+            for (int n = Ns; n < Ne; n++) {
+                varSum += SQ(pixelData->data.F32[n] - mean);
+            }
+            // variance on the mean (stdev / sqrt(N))^2
+            float varValue = varSum / (Npt - 1) / Npt;
+
+            // XXX this is probably not the correct variance to save
+            // the predicted variance should be the 1 / sum (1/var)
+
+            combined->image->data.F32[y][x] = mean;
+            combined->variance->data.F32[y][x] = varValue;
+            combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+        }
+    }
+  
+    psFree(pixelData);
+    psFree(pixelMask);
+
+    return (true);
+}
+
 /// Stack input images
 bool pmStackCombine(
@@ -1644 +1773 @@
         expweight = expmaps->variance;
     }
-
-/*     // Pre-reject inputs using KMM bimodality test. */
-/*     if (1)  { */
-/* /\*       KMMRejectUnpopular(input,x,y); *\/ */
-/*       rejection = true; */
-/*     } */
    
     // Set up rejection list
@@ -1660 +1783 @@
     for (int y = minInputRows; y < maxInputRows; y++) {
         for (int x = minInputCols; x < maxInputCols; x++) {
+            // this is only for TESTING:
             CHECKPIX(x, y, "Combining pixel %d,%d: %x %x %f %f %f %f %d %d %d\n", x, y, badMaskBits, blankMaskBits, iter, rej, sys, olympic, useVariance, safe, rejection);
-
-/*          // Pre-reject inputs using KMM bimodality test. */
-/*        if (1)  { */
-/*          KMMRejectUnpopular(input,x,y); */
-/* /\*      rejection = true; *\/ */
-/*        } */
-/*        else { */
-/*          KMMRejectBright(input,x,y); */
-/*        } */
-
-#ifdef TESTING
-            if (PS_SQR(x - TEST_X) + PS_SQR(y - TEST_Y) <= PS_SQR(TEST_RADIUS)) {
-                fprintf (stderr, "combine\n");
-            }
-# endif
 
             psVector *reject = NULL; // Images to reject for this pixel

trunk/psModules/src/imcombine/pmStack.h

@@ -47 +47 @@
                                 );
 
+bool pmStackCombineByPercentile(
+    pmReadout *combined,
+    psArray *input,
+    psF64 minRange,
+    psF64 maxRange,
+    psImageMaskType badMaskBits,        // treat these bits as 'bad'
+    psImageMaskType suspectMaskBits,    // treat these bits as 'suspect'
+    psImageMaskType blankMaskBits       // use this mask value for pixels missing input data (distinguish between Ninput = 0 and Ngood = 0?)
+  );
+
 /// Stack input images
 bool pmStackCombine(pmReadout *combined,///< Combined readout (output)

trunk/psModules/src/objects/pmSourceIO_CMF.c.in

@@ -57 +57 @@
 // followed by a zero-size matrix, followed by the table data
 
-bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
+/** this file is used to generate c-code for a number of different formats.  the perl program
+ ** 'mksource.pl' reads this template file and generates an output file for one of the named
+ ** formats, e.g., PS1_V1 or PS1_DV4 (see list in mksource.pl).  
+
+ ** Any instances of the word @CMFMODE@ are replaced by the named version of the format.
+
+ ** Any line which starts with a command of the form: @COMMAND[,COMMAND...]@ is kept or removed
+ ** from the output c-code depending on the commands.  These may be: ALL (lines are kept for
+ ** all formats) or the name of a format.  A line with just a plain format will be kept only
+ ** for that format.  The format name may be optionally proceded by >, <, <=, >=, ! or =.  In
+ ** these cases, the line is kept if the format is greater / lesser / not the same, etc as the
+ ** listed format.  The sequence of formats is defined for series (V, SV, DV) in mksource.pl so
+ ** that, e.g., >PS1_V3 would include e.g., PS1_V4 and PS1_V5, but not PS1_DV4.  The number in
+ ** a series may be replaced with ? to stand for the full set.  The format 
+ ** commands may be grouped on a line separated by commas.  Note that only ! (not) is
+ ** restrictive: only this rule will reduce the number of matches.  So for example,
+ ** @ALL,<PS1_DV3@ will have the effect of @ALL@ since any format >= PS1_DV3 will be matched by
+ ** ALL.  But @ALL,!PS1_DV3@ will generate a line for any format except PS1_DV3.
+
+ **/ 
+
+bool pmSourcesWrite_CMF_@CMFMODE@_New (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
 {
+    // fprintf (stderr, "writing with %s\n", __func__);
     PS_ASSERT_PTR_NON_NULL(fits, false);
     PS_ASSERT_PTR_NON_NULL(sources, false);
     PS_ASSERT_PTR_NON_NULL(extname, false);
-
-    psArray *table;
-    psMetadata *row;
 
     pmChip *chip = readout->parent->parent;
@@ -73 +92 @@
         pmSource *source = (pmSource *) sources->data[0];
         if (source->seq == -1) {
-            // let's write these out in S/N order
-            sources = psArraySort (sources, pmSourceSortByFlux);
+            // let's write these out in S/N order
+            sources = psArraySort (sources, pmSourceSortByFlux);
         } else {
-            sources = psArraySort (sources, pmSourceSortBySeq);
-        }
-    }
-
-    table = psArrayAllocEmpty (sources->n);
+            sources = psArraySort (sources, pmSourceSortBySeq);
+        }
+    }
 
     float magOffset; 
@@ -88 +105 @@
     pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
 
+    // before we can generate the table structure, we need to know if certain fields were measured.
+    // this information is not available a priori, so we need to check all sources
+    
+    bool haveLensingOBJsmear = false;
+    bool haveLensingOBJshear = false;
+    bool haveLensingPSFsmear = false;
+    bool haveLensingPSFshear = false;
+    bool haveLensingPSF      = false;
+
+    for (int i = 0; (i < sources->n) && !haveLensingOBJsmear && !haveLensingOBJshear && !haveLensingPSFsmear && !haveLensingPSFshear && !haveLensingPSF; i++) {
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        if (source->lensingOBJ && source->lensingOBJ->smear) haveLensingOBJsmear = true;
+        if (source->lensingOBJ && source->lensingOBJ->shear) haveLensingOBJshear = true;
+        if (source->lensingPSF && source->lensingPSF->smear) haveLensingPSFsmear = true;
+        if (source->lensingPSF && source->lensingPSF->shear) haveLensingPSFshear = true;
+        if (source->lensingPSF                             ) haveLensingPSF      = true;
+    }
+
+    /************ generate the table columns *****************/
+
+    // before we allocate the table, generate an empty array of table columns and generate them
+    psArray *tableColumns = psArrayAllocEmpty (100);
+
+    // add the named / typed columns to the collection of columns
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "IPP_IDET",         PS_DATA_U32); // "IPP detection identifier index"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "X_PSF",            PS_DATA_F32); // "PSF x coordinate"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "Y_PSF",            PS_DATA_F32); // "PSF y coordinate"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "X_PSF_SIG",        PS_DATA_F32); // "Sigma in PSF x coordinate"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "Y_PSF_SIG",        PS_DATA_F32); // "Sigma in PSF y coordinate"
+
+    // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
+    @PS1_V1@                    psFitsTableColumnAdd (tableColumns, "RA_PSF",           PS_DATA_F32); // "PSF RA coordinate (degrees)"
+    @PS1_V1@                    psFitsTableColumnAdd (tableColumns, "DEC_PSF",          PS_DATA_F32); // "PSF DEC coordinate (degrees)"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "POSANGLE",         PS_DATA_F32); // "position angle at source (degrees)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PLTSCALE",         PS_DATA_F32); // "plate scale at source (arcsec/pixel)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG",     PS_DATA_F32); // "PSF fit instrumental magnitude"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_INST_MAG_SIG", PS_DATA_F32); // "Sigma of PSF instrumental magnitude"
+
+    @ALL,!PS1_V1,!PS1_V2@       psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX",    PS_DATA_F32); // "PSF fit instrumental flux (counts)"
+    @ALL,!PS1_V1,!PS1_V2@       psFitsTableColumnAdd (tableColumns, "PSF_INST_FLUX_SIG", PS_DATA_F32); // "Sigma of PSF instrumental flux"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "AP_MAG",           PS_DATA_F32); // "magnitude in standard aperture"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "AP_MAG_RAW",       PS_DATA_F32); // "magnitude in reported aperture"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "AP_MAG_RADIUS",    PS_DATA_F32); // "radius used for aperture mags"
+    @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX",          PS_DATA_F32); // "instrumental flux in standard aperture"
+    @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableColumnAdd (tableColumns, "AP_FLUX_SIG",      PS_DATA_F32); // "aperture flux error"
+    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "AP_NPIX",          PS_DATA_S32); // "aperture unmasked pixels"
+
+    @<PS1_V3,PS1_SV1,PS1_DV?@   psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG",      PS_DATA_F32); // "PSF Magnitude using supplied calibration"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "CAL_PSF_MAG_SIG",  PS_DATA_F32); // "measured scatter of zero point calibration"
+
+    // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
+    @ALL,!PS1_V1@               psFitsTableColumnAdd (tableColumns, "RA_PSF",           PS_DATA_F64); // "PSF RA coordinate (degrees)"
+    @ALL,!PS1_V1@               psFitsTableColumnAdd (tableColumns, "DEC_PSF",          PS_DATA_F64); // "PSF DEC coordinate (degrees)"
+
+    @>=PS1_V3,>PS1_SV1@         psFitsTableColumnAdd (tableColumns, "PEAK_FLUX_AS_MAG", PS_DATA_F32); // "Peak flux expressed as magnitude"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "SKY",              PS_DATA_F32); // "Sky level"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "SKY_SIGMA",        PS_DATA_F32); // "Sigma of sky level"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_CHISQ",        PS_DATA_F32); // "Chisq of PSF-fit"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "CR_NSIGMA",        PS_DATA_F32); // "Nsigma deviations from PSF to CF"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "EXT_NSIGMA",       PS_DATA_F32); // "Nsigma deviations from PSF to EXT"
+
+    // PSF shape parameters:
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_MAJOR",        PS_DATA_F32); // "PSF width (major axis)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_MINOR",        PS_DATA_F32); // "PSF width (minor axis)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_THETA",        PS_DATA_F32); // "PSF orientation angle"
+    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_CORE",         PS_DATA_F32); // "k term if defined"
+    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MAJ",     PS_DATA_F32); // "PSF FWHM (major axis)"
+    @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableColumnAdd (tableColumns, "PSF_FWHM_MIN",     PS_DATA_F32); // "PSF FWHM (minor axis)"
+
+    // psf data quality
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_QF",           PS_DATA_F32); // "PSF coverage/quality factor (bad)"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "PSF_QF_PERFECT",   PS_DATA_F32); // "PSF coverage/quality factor (poor)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_NDOF",         PS_DATA_S32); // "degrees of freedom"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PSF_NPIX",         PS_DATA_S32); // "number of pixels in fit"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_XX",       PS_DATA_F32); // "second moments (X^2)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_XY",       PS_DATA_F32); // "second moments (X*Y)"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "MOMENTS_YY",       PS_DATA_F32); // "second moments (Y*Y)"
+    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M3C",      PS_DATA_F32); // "third momemt cos theta"
+    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M3S",      PS_DATA_F32); // "third momemt sin theta"
+    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M4C",      PS_DATA_F32); // "fourth momemt cos theta"
+    @>PS1_V2,PS1_SV?@           psFitsTableColumnAdd (tableColumns, "MOMENTS_M4S",      PS_DATA_F32); // "fourth momemt sin theta"
+  
+    // Lensing parameters are only written if they are measured
+    if (haveLensingOBJsmear) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_OBJ",  PS_DATA_F32); // "smear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_OBJ",   PS_DATA_F32); // "smear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_OBJ",   PS_DATA_F32); // "smear polarizability element (objects)"
+    }
+    if (haveLensingOBJshear) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_OBJ",  PS_DATA_F32); // "shear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_OBJ",   PS_DATA_F32); // "shear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_OBJ",   PS_DATA_F32); // "shear polarizability element (objects)"
+    }
+    if (false) {
+        // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_OBJ",      PS_DATA_F32); // "shear polarizability element (objects)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_OBJ",      PS_DATA_F32); // "shear polarizability element (objects)"
+    }
+
+    if (haveLensingPSFsmear) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SM_PSF",  PS_DATA_F32); // "smear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SM_PSF",   PS_DATA_F32); // "smear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SM_PSF",   PS_DATA_F32); // "smear polarizability element (PSFs)"
+    }
+    if (haveLensingPSFshear) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X11_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X12_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_X22_SH_PSF",  PS_DATA_F32); // "shear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_SH_PSF",   PS_DATA_F32); // "shear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_SH_PSF",   PS_DATA_F32); // "shear polarizability element (PSFs)"
+    }
+    if (haveLensingPSF) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E1_PSF",      PS_DATA_F32); // "shear polarizability element (PSFs)"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "LENS_E2_PSF",      PS_DATA_F32); // "shear polarizability element (PSFs)"
+    }
+
+    // if lensing params exist also include the backmapped chipID and chip coordinates
+    if (haveLensingPSF) {
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM",     PS_DATA_S16); // "id of warp input chip"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X",       PS_DATA_S16); // "x coord in warp input chip"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y",       PS_DATA_S16); // "y coord in warp input chip"
+        @>PS1_V4@               psFitsTableColumnAdd (tableColumns, "PADDING3",         PS_DATA_S16); // "more padding"
+    }
+
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "MOMENTS_R1",       PS_DATA_F32); // "first radial moment"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "MOMENTS_RH",       PS_DATA_F32); // "half radial moment"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX",        PS_DATA_F32); // "Kron Flux (in 2.5 R1)"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_ERR",    PS_DATA_F32); // "Kron Flux Error"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_INNER",  PS_DATA_F32); // "Kron Flux (in 1.0 R1)"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "KRON_FLUX_OUTER",  PS_DATA_F32); // "Kron Flux (in 4.0 R1)"
+
+    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_RAD",    PS_DATA_F32); // "Radius where object hits sky"
+    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_FLUX",   PS_DATA_F32); // "Flux / pix where object hits sky"
+    @>PS1_V3@                   psFitsTableColumnAdd (tableColumns, "SKY_LIMIT_SLOPE",  PS_DATA_F32); // "d(Flux/pix)/dRadius where object hits sky"
+
+    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_NUM",     PS_DATA_S16); // "id of warp input chip"
+    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_X",       PS_DATA_S16); // "x coord in warp input chip"
+    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "SRC_CHIP_Y",       PS_DATA_S16); // "y coord in warp input chip"
+    @>PS1_DV4@                  psFitsTableColumnAdd (tableColumns, "PADDING3",         PS_DATA_S16); // "more padding"
+
+    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NPOS",        PS_DATA_S32); // "nPos (n pix > 3 sigma)"
+    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_FRATIO",      PS_DATA_F32); // "fPos / (fPos + fNeg)"
+    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_BAD",  PS_DATA_F32); // "nPos / (nPos + nNeg)"
+    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_MASK", PS_DATA_F32); // "nPos / (nPos + nMask)"
+    @PS1_DV?@                   psFitsTableColumnAdd (tableColumns, "DIFF_NRATIO_ALL",  PS_DATA_F32); // "nPos / (nGood + nMask + nBad)"
+
+    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_R_P",         PS_DATA_F32); // "distance to positive match source"
+    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_SN_P",        PS_DATA_F32); // "signal-to-noise of pos match src"
+    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_R_M",         PS_DATA_F32); // "distance to negative match source"
+    @>PS1_DV1@                  psFitsTableColumnAdd (tableColumns, "DIFF_SN_M",        PS_DATA_F32); // "signal-to-noise of neg match src"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "FLAGS",            PS_DATA_U32); // "psphot analysis flags"
+    @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableColumnAdd (tableColumns, "FLAGS2",           PS_DATA_U32); // "psphot analysis flags"
+    @PS1_V3,PS1_SV2@            psFitsTableColumnAdd (tableColumns, "PADDING2",         PS_DATA_S32); // "more padding"
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "N_FRAMES",         PS_DATA_U16); // "Number of frames overlapping source center"
+
+    @ALL@                       psFitsTableColumnAdd (tableColumns, "PADDING",          PS_DATA_S16); // "more padding"
+
+    // note that this definition is inconsistent with the definition in
+    // Ohana/src/libautocode/def.  This version creates a table with data not
+    // properly aligned with the 8-byte boundaries. The structure defined by
+    // libautocode does this, but has a different order of elements (and adds
+    // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
+    // I'll leave it. But in future a PS1_SV2 should be forced to match
+    // libautocode. Note that addstar knows to detect the alternate version of
+    // PS1_SV1 and correctly interpret its fields.
+
+    // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
+    // subtracted
+
+    // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
+    // edge; 3) any pixels in the 3x3 peak region are masked;
+
+    // generate an FITS table using the array of columns defined above, with space for all rows
+    psFitsTable *table = psFitsTableCreate (tableColumns, sources->n);
+
+    /************ write the data to the table *****************/
+
     // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
     // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
     for (int i = 0; i < sources->n; i++) {
-        // this is the source associated with this image
+        // this is the source associated with this image
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // If source->seq is -1, source was generated in this analysis.  If source->seq is
@@ -102 +313 @@
         // generated on Alloc, and would thus be wrong for read in sources.
         if (source->seq == -1) {
-            source->seq = i;
-        }
-
-        // set the 'best' values for various output fields:
-        pmSourceOutputs outputs;
-        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
-
-        pmSourceOutputsMoments moments;
-        pmSourceOutputsSetMoments (&moments, source);
-
-        @PS1_DV?@ pmSourceDiffStats diffStats;
-        @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
-        @PS1_DV?@ if (source->diffStats) {
-        @PS1_DV?@     diffStats = *source->diffStats;
-        @PS1_DV?@ }
-
-        row = psMetadataAlloc ();
+            source->seq = i;
+        }
+
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
+
+        @PS1_DV?@ pmSourceDiffStats diffStats;
+        @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
+        @PS1_DV?@ if (source->diffStats) { diffStats = *source->diffStats;}
 
         // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
         // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment"
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
-
-        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
-        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
-        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
-
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
-
-        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
-        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
-
-        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
-        @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
-        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              source->apRadius);
-        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
-        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
-        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX",          PS_DATA_S32, "aperture unmasked pixels",                   source->apNpixels);
-
-        @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
-
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
-        
-        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
-        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
-        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
-
-        @>=PS1_V3,>PS1_SV1@       psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
-
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
+        @ALL@                       psFitsTableSetU32 (table, i, "IPP_IDET",          source->seq);                     // "IPP detection identifier index",            
+        @ALL@                       psFitsTableSetF32 (table, i, "X_PSF",             outputs.xPos);                    // "PSF x coordinate",                          
+        @ALL@                       psFitsTableSetF32 (table, i, "Y_PSF",             outputs.yPos);                    // "PSF y coordinate",                          
+        @ALL@                       psFitsTableSetF32 (table, i, "X_PSF_SIG",         outputs.xErr);                    // "Sigma in PSF x coordinate",                 
+        @ALL@                       psFitsTableSetF32 (table, i, "Y_PSF_SIG",         outputs.yErr);                    // "Sigma in PSF y coordinate",                 
+
+        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision                           
+        @PS1_V1@                    psFitsTableSetF32 (table, i, "RA_PSF",            outputs.ra);                      // "PSF RA coordinate (degrees)",                
+        @PS1_V1@                    psFitsTableSetF32 (table, i, "DEC_PSF",           outputs.dec);                     // "PSF DEC coordinate (degrees)",               
+
+        @ALL@                       psFitsTableSetF32 (table, i, "POSANGLE",          outputs.posAngle);                // "position angle at source (degrees)",        
+        @ALL@                       psFitsTableSetF32 (table, i, "PLTSCALE",          outputs.pltScale);                // "plate scale at source (arcsec/pixel)",      
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_INST_MAG",      source->psfMag);                  // "PSF fit instrumental magnitude",            
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_INST_MAG_SIG",  source->psfMagErr);               // "Sigma of PSF instrumental magnitude",       
+
+        @ALL,!PS1_V1,!PS1_V2@       psFitsTableSetF32 (table, i, "PSF_INST_FLUX",     source->psfFlux);                 // "PSF fit instrumental flux (counts)",        
+        @ALL,!PS1_V1,!PS1_V2@       psFitsTableSetF32 (table, i, "PSF_INST_FLUX_SIG", source->psfFluxErr);              // "Sigma of PSF instrumental flux",            
+
+        @ALL@                       psFitsTableSetF32 (table, i, "AP_MAG",            source->apMag);                   // "magnitude in standard aperture",            
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "AP_MAG_RAW",        source->apMagRaw);                // "magnitude in reported aperture",            
+        @ALL@                       psFitsTableSetF32 (table, i, "AP_MAG_RADIUS",     source->apRadius);                // "radius used for aperture mags",             
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX",           source->apFlux);                  // "instrumental flux in standard aperture",    
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psFitsTableSetF32 (table, i, "AP_FLUX_SIG",       source->apFluxErr);               // "aperture flux error",                       
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetS32 (table, i, "AP_NPIX",           source->apNpixels);               // "aperture unmasked pixels",                  
+
+        @<PS1_V3,PS1_SV1,PS1_DV?@   psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG",  outputs.peakMag);                 // "Peak flux expressed as magnitude"
+
+        @ALL@                       psFitsTableSetF32 (table, i, "CAL_PSF_MAG",       outputs.calMag);                  // "PSF Magnitude using supplied calibration",  
+        @ALL@                       psFitsTableSetF32 (table, i, "CAL_PSF_MAG_SIG",   zeroptErr);                       // "measured scatter of zero point calibration",
+
+        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...                                      
+        @ALL,!PS1_V1@               psFitsTableSetF64 (table, i, "RA_PSF",            outputs.ra);                      // "PSF RA coordinate (degrees)",               
+        @ALL,!PS1_V1@               psFitsTableSetF64 (table, i, "DEC_PSF",           outputs.dec);                     // "PSF DEC coordinate (degrees)",              
+
+        @>=PS1_V3,>PS1_SV1@         psFitsTableSetF32 (table, i, "PEAK_FLUX_AS_MAG",  outputs.peakMag);                 // "Peak flux expressed as magnitude",          
+        @ALL@                       psFitsTableSetF32 (table, i, "SKY",               source->sky);                     // "Sky level",                                 
+        @ALL@                       psFitsTableSetF32 (table, i, "SKY_SIGMA",         source->skyErr);                  // "Sigma of sky level",                        
+
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_CHISQ",         outputs.chisq);                   // "Chisq of PSF-fit",                          
+        @ALL@                       psFitsTableSetF32 (table, i, "CR_NSIGMA",         source->crNsigma);                // "Nsigma deviations from PSF to CF",          
+        @ALL@                       psFitsTableSetF32 (table, i, "EXT_NSIGMA",        source->extNsigma);               // "Nsigma deviations from PSF to EXT",         
 
         // PSF shape parameters:
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_MAJOR",         outputs.psfMajor);                // "PSF width (major axis)",                    
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_MINOR",         outputs.psfMinor);                // "PSF width (minor axis)",                    
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_THETA",         outputs.psfTheta);                // "PSF orientation angle",                     
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_CORE",          outputs.psfCore);                 // "k term if defined",                         
+
+        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MAJ",      outputs.psfMajorFWHM);            // "PSF FWHM (major axis)",                     
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psFitsTableSetF32 (table, i, "PSF_FWHM_MIN",      outputs.psfMinorFWHM);            // "PSF FWHM (minor axis)",                     
+
+        // psf data quality
+        @ALL@                       psFitsTableSetF32 (table, i, "PSF_QF",            source->pixWeightNotBad);         // "PSF coverage/quality factor (bad)",         
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "PSF_QF_PERFECT",    source->pixWeightNotPoor);        // "PSF coverage/quality factor (poor)",        
+        @ALL@                       psFitsTableSetS32 (table, i, "PSF_NDOF",          outputs.nDOF);                    // "degrees of freedom",                        
+        @ALL@                       psFitsTableSetS32 (table, i, "PSF_NPIX",          outputs.nPix);                    // "number of pixels in fit",                   
+
+        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_XX",        moments.Mxx);                     // "second moments (X^2)",                      
+        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_XY",        moments.Mxy);                     // "second moments (X*Y)",                      
+        @ALL@                       psFitsTableSetF32 (table, i, "MOMENTS_YY",        moments.Myy);                     // "second moments (Y*Y)",                      
+
+        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M3C",       moments.M_c3);                    // "third momemt cos theta",                    
+        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M3S",       moments.M_s3);                    // "third momemt sin theta",                    
+        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M4C",       moments.M_c4);                    // "fourth momemt cos theta",                   
+        @>PS1_V2,PS1_SV?@           psFitsTableSetF32 (table, i, "MOMENTS_M4S",       moments.M_s4);                    // "fourth momemt sin theta",                   
+
+        // Lensing parameters:
+        if (source->lensingOBJ && source->lensingOBJ->smear) {
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SM_OBJ",   source->lensingOBJ->smear->X11);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SM_OBJ",   source->lensingOBJ->smear->X12);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SM_OBJ",   source->lensingOBJ->smear->X22);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SM_OBJ",    source->lensingOBJ->smear->e1);   // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SM_OBJ",    source->lensingOBJ->smear->e2);   // "smear polarizability element (objects)",     
+        }
+        if (source->lensingOBJ && source->lensingOBJ->shear) {
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SH_OBJ",   source->lensingOBJ->shear->X11);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SH_OBJ",   source->lensingOBJ->shear->X12);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SH_OBJ",   source->lensingOBJ->shear->X22);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SH_OBJ",    source->lensingOBJ->shear->e1);   // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SH_OBJ",    source->lensingOBJ->shear->e2);   // "shear polarizability element (objects)",     
+        }
+        if (false && source->lensingOBJ) {
+            // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_PSF",       source->lensingOBJ->e1);          // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_PSF",       source->lensingOBJ->e2);          // "shear polarizability element (objects)",     
+        }
+
+        if (source->lensingPSF && source->lensingPSF->smear) {
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SM_PSF",   source->lensingPSF->smear->X11);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SM_PSF",   source->lensingPSF->smear->X12);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SM_PSF",   source->lensingPSF->smear->X22);  // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SM_PSF",    source->lensingPSF->smear->e1);   // "smear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SM_PSF",    source->lensingPSF->smear->e2);   // "smear polarizability element (objects)",     
+        }
+        if (source->lensingPSF && source->lensingPSF->shear) {
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X11_SH_PSF",   source->lensingPSF->shear->X11);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X12_SH_PSF",   source->lensingPSF->shear->X12);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_X22_SH_PSF",   source->lensingPSF->shear->X22);  // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_SH_PSF",    source->lensingPSF->shear->e1);   // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_SH_PSF",    source->lensingPSF->shear->e2);   // "shear polarizability element (objects)",     
+        }
+        if (source->lensingPSF) {
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E1_PSF",       source->lensingPSF->e1);          // "shear polarizability element (objects)",     
+            @>PS1_V4@               psFitsTableSetF32 (table, i, "LENS_E2_PSF",       source->lensingPSF->e2);          // "shear polarizability element (objects)",     
+        }
+
+        // if lensing params exist also include the backmapped chipID and chip coordinates
+        if (source->lensingPSF && source->lensingPSF->shear) {
+            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_NUM",      source->chipNum);                 // "id of warp input chip",       
+            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_X",        source->chipX);                   // "x coord in warp input chip",  
+            @>PS1_V4@               psFitsTableSetS16 (table, i, "SRC_CHIP_Y",        source->chipY);                   // "y coord in warp input chip",  
+            @>PS1_V4@               psFitsTableSetS16 (table, i, "PADDING3",          0);                               // "more padding",                
+        }
+
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "MOMENTS_R1",        moments.Mrf);                     // "first radial moment",                        
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "MOMENTS_RH",        moments.Mrh);                     // "half radial moment",                         
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX",         moments.Krf);                     // "Kron Flux (in 2.5 R1)",                      
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_ERR",     moments.dKrf);                    // "Kron Flux Error",                            
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_INNER",   moments.Kinner);                  // "Kron Flux (in 1.0 R1)",                      
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psFitsTableSetF32 (table, i, "KRON_FLUX_OUTER",   moments.Kouter);                  // "Kron Flux (in 4.0 R1)",                      
+        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_RAD",     source->skyRadius);               // "Radius where object hits sky",               
+        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_FLUX",    source->skyFlux);                 // "Flux / pix where object hits sky",           
+        @>PS1_V3@                   psFitsTableSetF32 (table, i, "SKY_LIMIT_SLOPE",   source->skySlope);                // "d(Flux/pix)/dRadius where object hits sky",  
+
+        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_NUM",      source->chipNum);                 // "id of warp input chip",       
+        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_X",        source->chipX);                   // "x coord in warp input chip",  
+        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "SRC_CHIP_Y",        source->chipY);                   // "y coord in warp input chip",  
+        @>PS1_DV4@                  psFitsTableSetS16 (table, i, "PADDING3",          0);                               // "more padding",                
+
+        @PS1_DV?@                   psFitsTableSetS32 (table, i, "DIFF_NPOS",         diffStats.nGood);                 // "nPos (n pix > 3 sigma)",                     
+        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_FRATIO",       diffStats.fRatio);                // "fPos / (fPos + fNeg)",                       
+        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_BAD",   diffStats.nRatioBad);             // "nPos / (nPos + nNeg)",                       
+        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_MASK",  diffStats.nRatioMask);            // "nPos / (nPos + nMask)",                      
+        @PS1_DV?@                   psFitsTableSetF32 (table, i, "DIFF_NRATIO_ALL",   diffStats.nRatioAll);             // "nPos / (nGood + nMask + nBad)",              
+
+        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_R_P",          diffStats.Rp);                    // "distance to positive match source",  
+        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_SN_P",         diffStats.SNp);                   // "signal-to-noise of pos match src",   
+        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_R_M",          diffStats.Rm);                    // "distance to negative match source",  
+        @>PS1_DV1@                  psFitsTableSetF32 (table, i, "DIFF_SN_M",         diffStats.SNm);                   // "signal-to-noise of neg match src",   
+
+        @ALL@                      psFitsTableSetU32 (table, i, "FLAGS",              source->mode);                    // "psphot analysis flags",                      
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ psFitsTableSetU32 (table, i, "FLAGS2",             source->mode2);                   // "psphot analysis flags",                      
+        @PS1_V3,PS1_SV2@           psFitsTableSetS32 (table, i, "PADDING2",           0);                               // "padding",                                    
+
+        @ALL@                      psFitsTableSetU16 (table, i, "N_FRAMES",           source->nFrames);                 // "Number of frames overlapping source center", 
+        @ALL@                      psFitsTableSetS16 (table, i, "PADDING",            0);                               // "padding",                                    
+    }
+
+    // XXX why do we make a copy here to be supplemented with the masks?  why not do this in the calling function?
+    psMetadata *header = psMetadataCopy(NULL, tableHeader);
+    pmSourceMasksHeader(header);
+
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTableNew(fits, header, table, extname)) {
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+    psFree(table);
+    psFree(header);
+
+    return true;
+}
+
+bool pmSourcesWrite_CMF_@CMFMODE@_Old (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
+{
+    // fprintf (stderr, "writing with %s\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+    PS_ASSERT_PTR_NON_NULL(extname, false);
+
+    psArray *table;
+    psMetadata *row;
+
+    pmChip *chip = readout->parent->parent;
+
+    // if the sequence is defined, write these in seq order; otherwise
+    // write them in S/N order:
+    if (sources->n > 0) {
+        pmSource *source = (pmSource *) sources->data[0];
+        if (source->seq == -1) {
+            // let's write these out in S/N order
+            sources = psArraySort (sources, pmSourceSortByFlux);
+        } else {
+            sources = psArraySort (sources, pmSourceSortBySeq);
+        }
+    }
+
+    table = psArrayAllocEmpty (sources->n);
+
+    float magOffset; 
+    float zeroptErr; 
+    float fwhmMajor; 
+    float fwhmMinor;
+    pmSourceOutputsCommonValues (&magOffset, &zeroptErr, &fwhmMajor, &fwhmMinor, readout, imageHeader);
+
+    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
+    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
+    for (int i = 0; i < sources->n; i++) {
+        // this is the source associated with this image
+        pmSource *thisSource = sources->data[i];
+
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        // If source->seq is -1, source was generated in this analysis.  If source->seq is
+        // not -1, source was read from elsewhere: in the latter case, preserve the source
+        // ID.  source.seq is used instead of source.id since the latter is a const
+        // generated on Alloc, and would thus be wrong for read in sources.
+        if (source->seq == -1) {
+            source->seq = i;
+        }
+
+        // set the 'best' values for various output fields:
+        pmSourceOutputs outputs;
+        pmSourceOutputsSetValues (&outputs, source, chip, fwhmMajor, fwhmMinor, magOffset);
+
+        pmSourceOutputsMoments moments;
+        pmSourceOutputsSetMoments (&moments, source);
+
+        @PS1_DV?@ pmSourceDiffStats diffStats;
+        @PS1_DV?@ pmSourceDiffStatsInit(&diffStats);
+        @PS1_DV?@ if (source->diffStats) {
+        @PS1_DV?@     diffStats = *source->diffStats;
+        @PS1_DV?@ }
+
+        row = psMetadataAlloc ();
+
+        // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
+        // This set of psMetadataAdd Entries marks the "----" "Start of the PSF segment"
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           outputs.xPos);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           outputs.yPos);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  outputs.xErr);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  outputs.yErr);
+
+        // NOTE: pre-PS1_V2, we only reported RA & DEC in floats for reference, not precision
+        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F32, "PSF RA coordinate (degrees)",                outputs.ra);
+        @PS1_V1@                  psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F32, "PSF DEC coordinate (degrees)",               outputs.dec);
+
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         outputs.posAngle);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       outputs.pltScale);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        source->psfMagErr);
+
+        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental flux (counts)",         source->psfFlux);
+        @ALL,!PS1_V1,!PS1_V2@     psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental flux",             source->psfFluxErr);
+
+        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in reported aperture",             source->apMagRaw);
+        @ALL@                       psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              source->apRadius);
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",     source->apFlux);
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",                        source->apFluxErr);
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "AP_NPIX",          PS_DATA_S32, "aperture unmasked pixels",                   source->apNpixels);
+
+        @<PS1_V3,PS1_SV1,PS1_DV?@ psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   outputs.calMag);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
+        
+        // NOTE: RA & DEC (both double) need to be on an 8-byte boundary...
+        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                outputs.ra);
+        @ALL,!PS1_V1@             psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               outputs.dec);
+
+        @>=PS1_V3,>PS1_SV1@       psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           outputs.peakMag);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
+
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           outputs.chisq);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
+
+        // PSF shape parameters:
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     outputs.psfMajor);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     outputs.psfMinor);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      outputs.psfTheta);
         @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_CORE",         PS_DATA_F32, "k term if defined",                          outputs.psfCore);
 
-        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
+        // I use a look-up table and linear interpolation to map PSF_MAJOR,PSF_MINOR + PSF_CORE to FWHM values
         @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MAJ",        PS_DATA_F32, "PSF FWHM (major axis)",                   outputs.psfMajorFWHM);
         @>PS1_V4,>PS1_SV2,>PS1_DV3@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_FWHM_MIN",        PS_DATA_F32, "PSF FWHM (minor axis)",                   outputs.psfMinorFWHM);
 
         // psf data quality
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
-        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
-
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
-
-        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
-        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
-        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
-        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
-
-        // Lensing parameters:
-        if (source->lensingOBJ && source->lensingOBJ->smear) {
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
-        }
-        if (source->lensingOBJ && source->lensingOBJ->shear) {
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1); 
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2); 
-        }
-        if (false && source->lensingOBJ) {
-          // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e1); 
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e2); 
-        }
-
-        if (source->lensingPSF && source->lensingPSF->smear) {
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1); 
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2); 
-        }
-        if (source->lensingPSF && source->lensingPSF->shear) {
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1); 
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2); 
-        }
-        if (source->lensingPSF) {
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e1); 
-          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e2); 
-        }
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor (bad)",          source->pixWeightNotBad);
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",         source->pixWeightNotPoor);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         outputs.nDOF);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    outputs.nPix);
+
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                       moments.Mxx);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                       moments.Mxy);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                       moments.Myy);
+
+        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3C",      PS_DATA_F32, "third momemt cos theta",                     moments.M_c3);
+        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M3S",      PS_DATA_F32, "third momemt sin theta",                     moments.M_s3);
+        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4C",      PS_DATA_F32, "fourth momemt cos theta",                    moments.M_c4);
+        @>PS1_V2,PS1_SV?@         psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_M4S",      PS_DATA_F32, "fourth momemt sin theta",                    moments.M_s4);
+
+        // Lensing parameters:
+        if (source->lensingOBJ && source->lensingOBJ->smear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_OBJ",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e1);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_OBJ",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingOBJ->smear->e2);
+        }
+        if (source->lensingOBJ && source->lensingOBJ->shear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_OBJ",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_OBJ",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->shear->e2); 
+        }
+        if (false && source->lensingOBJ) {
+          // do not bother to save these as they are equivalent to Mxx,Mxy,Myy above
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingOBJ->e2); 
+        }
+
+        if (source->lensingPSF && source->lensingPSF->smear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SM_PSF",  PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SM_PSF",   PS_DATA_F32, "smear polarizability element (objects)",     source->lensingPSF->smear->e2); 
+        }
+        if (source->lensingPSF && source->lensingPSF->shear) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X11_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X11);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X12_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X12);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_X22_SH_PSF",  PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->X22);
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_SH_PSF",   PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->shear->e2); 
+        }
+        if (source->lensingPSF) {
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E1_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e1); 
+          @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "LENS_E2_PSF",      PS_DATA_F32, "shear polarizability element (objects)",     source->lensingPSF->e2); 
+        }
 
         // if lensing params exist also include the backmapped chipID and chip coordinates
-        if (source->lensingPSF && source->lensingPSF->shear) {
+        if (source->lensingPSF && source->lensingPSF->shear) {
             @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",   PS_DATA_S16, "id of warp input chip",     source->chipNum); 
             @>PS1_V4@               psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",    PS_DATA_S16, "x coord in warp input chip",     source->chipX); 
@@ -244 +683 @@
         @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",                      moments.Kouter);
 
-        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
-        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
-        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
-
-        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",     PS_DATA_S16, "id of warp input chip",                      source->chipNum); 
-        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",       PS_DATA_S16, "x coord in warp input chip",                 source->chipX); 
-        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y",       PS_DATA_S16, "y coord in warp input chip",                 source->chipY); 
-        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "PADDING3",         PS_DATA_S16, "more padding", 0);
-
-        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
-        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
-        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
-        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
-        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
-
-        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
-        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",           diffStats.SNp);
-        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",          diffStats.Rm);
-        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",           diffStats.SNm);
-
-        @ALL@                      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
-        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
-        @PS1_V3,PS1_SV2@           psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
-        @PS1_SV?@
-
-          // note that this definition is inconsistent with the definition in
-          // Ohana/src/libautocode/def.  This version creates a table with data not
-          // properly aligned with the 8-byte boundaries. The structure defined by
-          // libautocode does this, but has a different order of elements (and adds
-          // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
-          // I'll leave it. But in future a PS1_SV2 should be forced to match
-          // libautocode. Note that addstar knows to detect the alternate version of
-          // PS1_SV1 and correctly interpret its fields.
-
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
-        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
+        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_RAD",    PS_DATA_F32, "Radius where object hits sky",               source->skyRadius);
+        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_FLUX",   PS_DATA_F32, "Flux / pix where object hits sky",           source->skyFlux);
+        @>PS1_V3@                 psMetadataAdd (row, PS_LIST_TAIL, "SKY_LIMIT_SLOPE",  PS_DATA_F32, "d(Flux/pix)/dRadius where object hits sky",  source->skySlope);
+
+        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_NUM",     PS_DATA_S16, "id of warp input chip",                      source->chipNum); 
+        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_X",       PS_DATA_S16, "x coord in warp input chip",                 source->chipX); 
+        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "SRC_CHIP_Y",       PS_DATA_S16, "y coord in warp input chip",                 source->chipY); 
+        @>PS1_DV4@                psMetadataAdd (row, PS_LIST_TAIL, "PADDING3",         PS_DATA_S16, "more padding", 0);
+
+        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                     diffStats.nGood);
+        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                       diffStats.fRatio);
+        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                       diffStats.nRatioBad);
+        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                      diffStats.nRatioMask);
+        @PS1_DV?@                 psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",              diffStats.nRatioAll);
+
+        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",          diffStats.Rp);
+        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",           diffStats.SNp);
+        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",          diffStats.Rm);
+        @>PS1_DV1@                        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",           diffStats.SNm);
+
+        @ALL@                      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags",                      source->mode2);
+        @PS1_V3,PS1_SV2@           psMetadataAdd (row, PS_LIST_TAIL, "PADDING2",         PS_DATA_S32, "more padding", 0);
+        @PS1_SV?@
+
+          // note that this definition is inconsistent with the definition in
+          // Ohana/src/libautocode/def.  This version creates a table with data not
+          // properly aligned with the 8-byte boundaries. The structure defined by
+          // libautocode does this, but has a different order of elements (and adds
+          // padding2 to fix things). We have generated may files with PS1_SV1 as is, so
+          // I'll leave it. But in future a PS1_SV2 should be forced to match
+          // libautocode. Note that addstar knows to detect the alternate version of
+          // PS1_SV1 and correctly interpret its fields.
+
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", source->nFrames);
+        @ALL@                     psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
 
         psArrayAdd (table, 100, row);
@@ -320 +759 @@
 }
 
+bool pmSourcesWrite_CMF_@CMFMODE@ (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname, psMetadata *recipe)
+{
+  // bool status = pmSourcesWrite_CMF_@CMFMODE@_Old (fits, readout, sources, imageHeader, tableHeader, extname, recipe);
+     bool status = pmSourcesWrite_CMF_@CMFMODE@_New (fits, readout, sources, imageHeader, tableHeader, extname, recipe);
+    return status;
+}
+
 // read in a readout from the fits file
-psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header)
+psArray *pmSourcesRead_CMF_@CMFMODE@_New (psFits *fits, psMetadata *header)
 {
+    // fprintf (stderr, "reading with %s\n", __func__);
+
     PS_ASSERT_PTR_NON_NULL(fits, false);
     PS_ASSERT_PTR_NON_NULL(header, false);
@@ -355 +803 @@
     psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
 
-    // convert the table to the pmSource entriesa
+    // reads full table into memory
+    psFitsTable *table = psFitsReadTableNew (fits);
+
+    // convert the table to the pmSource entries
+    for (int i = 0; i < numSources; i++) {
+        pmSource *source = pmSourceAlloc ();
+        pmModel *model = pmModelAlloc (modelType);
+        source->modelPSF  = model;
+        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
+
+        // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
+        PAR = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
+
+        @ALL@                       source->seq                 = psFitsTableGetU32 (&status, table, i, "IPP_IDET");
+        @ALL@                       PAR[PM_PAR_XPOS]            = psFitsTableGetF32 (&status, table, i, "X_PSF");
+        @ALL@                       PAR[PM_PAR_YPOS]            = psFitsTableGetF32 (&status, table, i, "Y_PSF");
+        @ALL@                       dPAR[PM_PAR_XPOS]           = psFitsTableGetF32 (&status, table, i, "X_PSF_SIG");
+        @ALL@                       dPAR[PM_PAR_YPOS]           = psFitsTableGetF32 (&status, table, i, "Y_PSF_SIG");
+        @ALL@                       axes.major                  = psFitsTableGetF32 (&status, table, i, "PSF_MAJOR");
+        @ALL@                       axes.minor                  = psFitsTableGetF32 (&status, table, i, "PSF_MINOR");
+        @ALL@                       axes.theta                  = psFitsTableGetF32 (&status, table, i, "PSF_THETA");
+        @ALL@                       axes.theta                  = axes.theta * PS_RAD_DEG;
+        
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7]           = psFitsTableGetF32 (&status, table, i, "PSF_CORE");
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ } 
+
+        @ALL@                       PAR[PM_PAR_SKY]             = psFitsTableGetF32 (&status, table, i, "SKY");
+        @ALL@                       dPAR[PM_PAR_SKY]            = psFitsTableGetF32 (&status, table, i, "SKY_SIGMA");
+        @ALL@                       source->sky                 = PAR[PM_PAR_SKY];
+        @ALL@                       source->skyErr              = dPAR[PM_PAR_SKY];
+
+        // XXX use these to determine PAR[PM_PAR_I0]?
+        @ALL@                       source->psfMag              = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG");
+        @ALL@                       source->psfMagErr           = psFitsTableGetF32 (&status, table, i, "PSF_INST_MAG_SIG");
+        @ALL@                       source->apMag               = psFitsTableGetF32 (&status, table, i, "AP_MAG");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  source->apMagRaw            = psFitsTableGetF32 (&status, table, i, "AP_MAG_RAW");
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux              = psFitsTableGetF32 (&status, table, i, "AP_FLUX");
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr           = psFitsTableGetF32 (&status, table, i, "AP_FLUX_SIG");
+
+        // XXX use these to determine PAR[PM_PAR_I0] if they exist?
+        // XXX add these to PS1_SV1?
+        @>PS1_V2,PS1_SV?,PS1_DV?@   source->psfFlux             = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX");
+        @>PS1_V2,PS1_SV?,PS1_DV?@   source->psfFluxErr          = psFitsTableGetF32 (&status, table, i, "PSF_INST_FLUX_SIG");
+
+        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
+        @ALL@     PAR[PM_PAR_I0]                                = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
+        @ALL@     dPAR[PM_PAR_I0]                               = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->psfMagErr : NAN;
+
+        pmPSF_AxesToModel (PAR, axes, model->class->useReff);
+
+        @ALL@     float peakMag                                 = psFitsTableGetF32 (&status, table, i, "PEAK_FLUX_AS_MAG");
+        @ALL@     float peakFlux                                = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN;
+
+        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
+        @ALL@     source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
+        @ALL@     source->peak->rawFlux = peakFlux;
+        @ALL@     source->peak->smoothFlux = peakFlux;
+        @ALL@     source->peak->xf                              = PAR[PM_PAR_XPOS]; // more accurate position
+        @ALL@     source->peak->yf                              = PAR[PM_PAR_YPOS]; // more accurate position
+        @ALL@     source->peak->dx                              = dPAR[PM_PAR_XPOS];
+        @ALL@     source->peak->dy                              = dPAR[PM_PAR_YPOS];
+
+        @ALL@     source->pixWeightNotBad                       = psFitsTableGetF32 (&status, table, i, "PSF_QF");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor     = psFitsTableGetF32 (&status, table, i, "PSF_QF_PERFECT");
+        @ALL@     source->crNsigma                              = psFitsTableGetF32 (&status, table, i, "CR_NSIGMA");
+        @ALL@     source->extNsigma                             = psFitsTableGetF32 (&status, table, i, "EXT_NSIGMA");
+        @ALL@     source->apRadius                              = psFitsTableGetF32 (&status, table, i, "AP_MAG_RADIUS");
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels           = psFitsTableGetS32 (&status, table, i, "AP_NPIX");
+
+        // note that some older versions used PSF_PROBABILITY: this was not well defined.
+        @ALL@     model->chisq                                  = psFitsTableGetF32 (&status, table, i, "PSF_CHISQ");
+        @ALL@     model->nDOF                                   = psFitsTableGetS32 (&status, table, i, "PSF_NDOF");
+        @ALL@     model->nPix                                   = psFitsTableGetS32 (&status, table, i, "PSF_NPIX");
+
+        @ALL@     source->moments = pmMomentsAlloc ();
+        @ALL@     source->moments->Mx = source->peak->xf; // we don't have both Mx,My and xf,yf in the cmf
+        @ALL@     source->moments->My = source->peak->yf; // we don't have both Mx,My and xf,yf in the cmf
+
+        @ALL@     source->moments->Mxx                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_XX");
+        @ALL@     source->moments->Mxy                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_XY");
+        @ALL@     source->moments->Myy                          = psFitsTableGetF32 (&status, table, i, "MOMENTS_YY");
+
+        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
+        // we are storing enough information so the output will be consistent with the input
+        @>PS1_V2,PS1_SV?@ source->moments->Mxxx         = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3C");
+        @>PS1_V2,PS1_SV?@ source->moments->Mxxy         =  0.00;
+        @>PS1_V2,PS1_SV?@ source->moments->Mxyy         =  0.00;
+        @>PS1_V2,PS1_SV?@ source->moments->Myyy         = -1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M3S");
+        @>PS1_V2,PS1_SV?@ source->moments->Mxxxx        = +1.00 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4C");
+        @>PS1_V2,PS1_SV?@ source->moments->Mxxxy        =  0.00;
+        @>PS1_V2,PS1_SV?@ source->moments->Mxxyy        =  0.00;
+        @>PS1_V2,PS1_SV?@ source->moments->Mxyyy        = -0.25 * psFitsTableGetF32 (&status, table, i, "MOMENTS_M4S");
+        @>PS1_V2,PS1_SV?@ source->moments->Myyyy        =  0.00;
+
+        // Lensing parameters (on read if PS1_V5+)
+        if (haveLensOBJ) {
+            source->lensingOBJ = pmSourceLensingAlloc ();
+            source->lensingOBJ->smear = pmLensingParsAlloc();
+            source->lensingOBJ->shear = pmLensingParsAlloc();
+
+            @>PS1_V4@ source->lensingOBJ->smear->X11            = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_OBJ");
+            @>PS1_V4@ source->lensingOBJ->smear->X12            = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_OBJ");
+            @>PS1_V4@ source->lensingOBJ->smear->X22            = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_OBJ");
+            @>PS1_V4@ source->lensingOBJ->smear->e1             = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_OBJ");
+            @>PS1_V4@ source->lensingOBJ->smear->e2             = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_OBJ");
+            @>PS1_V4@ source->lensingOBJ->shear->X11            = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_OBJ");
+            @>PS1_V4@ source->lensingOBJ->shear->X12            = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_OBJ");
+            @>PS1_V4@ source->lensingOBJ->shear->X22            = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_OBJ");
+            @>PS1_V4@ source->lensingOBJ->shear->e1             = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_OBJ");
+            @>PS1_V4@ source->lensingOBJ->shear->e2             = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_OBJ");
+        }
+
+        @>PS1_V4@ source->chipNum                               = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_NUM");
+        @>PS1_V4@ source->chipX                                 = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_X");
+        @>PS1_V4@ source->chipY                                 = psFitsTableGetS16 (&status, table, i, "SRC_CHIP_Y");
+
+        if (haveLensPSF) {
+            source->lensingPSF = pmSourceLensingAlloc ();
+            source->lensingPSF->smear = pmLensingParsAlloc();
+            source->lensingPSF->shear = pmLensingParsAlloc();
+            
+            @>PS1_V4@            source->lensingPSF->smear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SM_PSF");
+            @>PS1_V4@            source->lensingPSF->smear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SM_PSF");
+            @>PS1_V4@            source->lensingPSF->smear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SM_PSF");
+            @>PS1_V4@            source->lensingPSF->smear->e1  = psFitsTableGetF32 (&status, table, i, "LENS_E1_SM_PSF");
+            @>PS1_V4@            source->lensingPSF->smear->e2  = psFitsTableGetF32 (&status, table, i, "LENS_E2_SM_PSF");
+            @>PS1_V4@            source->lensingPSF->shear->X11 = psFitsTableGetF32 (&status, table, i, "LENS_X11_SH_PSF");
+            @>PS1_V4@            source->lensingPSF->shear->X12 = psFitsTableGetF32 (&status, table, i, "LENS_X12_SH_PSF");
+            @>PS1_V4@            source->lensingPSF->shear->X22 = psFitsTableGetF32 (&status, table, i, "LENS_X22_SH_PSF");
+            @>PS1_V4@            source->lensingPSF->shear->e1  = psFitsTableGetF32 (&status, table, i, "LENS_E1_SH_PSF");
+            @>PS1_V4@            source->lensingPSF->shear->e2  = psFitsTableGetF32 (&status, table, i, "LENS_E2_SH_PSF");
+            @>PS1_V4@            source->lensingPSF->e1         = psFitsTableGetF32 (&status, table, i, "LENS_E1_PSF");
+            @>PS1_V4@            source->lensingPSF->e2         = psFitsTableGetF32 (&status, table, i, "LENS_E2_PSF");
+        }
+
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psFitsTableGetF32 (&status, table, i, "MOMENTS_R1");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrh         = psFitsTableGetF32 (&status, table, i, "MOMENTS_RH");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFlux    = psFitsTableGetF32 (&status, table, i, "KRON_FLUX");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFluxErr = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_ERR");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFinner  = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_INNER");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->KronFouter  = psFitsTableGetF32 (&status, table, i, "KRON_FLUX_OUTER");
+
+        @>PS1_V3@                  source->skyRadius            = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_RAD");
+        @>PS1_V3@                  source->skyFlux              = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_FLUX");
+        @>PS1_V3@                  source->skySlope             = psFitsTableGetF32 (&status, table, i, "SKY_LIMIT_SLOPE");
+
+        @PS1_DV?@                  int nPos                     = psFitsTableGetS32 (&status, table, i, "DIFF_NPOS");
+        @PS1_DV?@  if (nPos) {
+        @PS1_DV?@      source->diffStats                        = pmSourceDiffStatsAlloc();
+        @PS1_DV?@      source->diffStats->nGood                 = nPos;
+        @PS1_DV?@      source->diffStats->fRatio                = psFitsTableGetF32 (&status, table, i, "DIFF_FRATIO");
+        @PS1_DV?@      source->diffStats->nRatioBad             = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_BAD");
+        @PS1_DV?@      source->diffStats->nRatioMask            = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_MASK");
+        @PS1_DV?@      source->diffStats->nRatioAll             = psFitsTableGetF32 (&status, table, i, "DIFF_NRATIO_ALL");
+        
+        @>PS1_DV1@      source->diffStats->Rp                   = psFitsTableGetF32 (&status, table, i, "DIFF_R_P");
+        @>PS1_DV1@      source->diffStats->SNp                  = psFitsTableGetF32 (&status, table, i, "DIFF_SN_P");
+        @>PS1_DV1@      source->diffStats->Rm                   = psFitsTableGetF32 (&status, table, i, "DIFF_R_M");
+        @>PS1_DV1@      source->diffStats->SNm                  = psFitsTableGetF32 (&status, table, i, "DIFF_SN_M");
+        @PS1_DV?@  }
+
+        @ALL@                      source->mode                 = psFitsTableGetU32 (&status, table, i, "FLAGS");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->mode2                = psFitsTableGetU32 (&status, table, i, "FLAGS2");
+        @ALL@                      source->nFrames              = psFitsTableGetU16 (&status, table, i, "N_FRAMES");
+        assert (status);
+
+        sources->data[i] = source;
+    }
+    psFree (table);
+    return sources;
+}
+
+// read in a readout from the fits file
+psArray *pmSourcesRead_CMF_@CMFMODE@_Old (psFits *fits, psMetadata *header)
+{
+    // fprintf (stderr, "reading with %s\n", __func__);
+
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(header, false);
+
+    bool status;
+    psF32 *PAR, *dPAR;
+    psEllipseAxes axes;
+
+    // define PSF model type
+    int defaultModelType = pmModelClassGetType ("PS_MODEL_GAUSS");
+    int modelType = -1;
+
+    // if header does not define the model, default to a gaussian
+    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSFMODEL");
+    if (PSF_NAME != NULL) {
+        modelType = pmModelClassGetType (PSF_NAME);
+    }
+    // work around bug in psphotFullForce
+    if (modelType < 0) {
+        modelType = defaultModelType;
+    }
+    // assert (modelType > -1);
+
+    // do we expect to find lensing parameters?
+    bool haveLensOBJ = psMetadataLookupBool (&status, header, "LENS_OBJ");
+    bool haveLensPSF = psMetadataLookupBool (&status, header, "LENS_PSF");
+
+    // We get the size of the table, and allocate the array of sources first because the table
+    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
+    // the table is read blocks the free.  In fact, it's better to read the table row by row.
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
+
+    // convert the table to the pmSource entries
     for (int i = 0; i < numSources; i++) {
         psMetadata *row = psFitsReadTableRow(fits, i); // Table row
@@ -382 +1041 @@
         @ALL@     axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");
         @ALL@     axes.theta        = axes.theta * PS_RAD_DEG;
-        
-        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
-        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
-        @>PS1_V4,>PS1_SV2,>PS1_DV3@ } 
+        
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ if (model->params->n > PM_PAR_7) {
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@     PAR[PM_PAR_7] = psMetadataLookupF32 (&status, row, "PSF_CORE");
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ } 
 
         @ALL@     PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
@@ -396 +1055 @@
         @ALL@     source->psfMagErr = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
         @ALL@     source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
-        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
-        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux = psMetadataLookupF32 (&status, row, "AP_FLUX");
-        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@  source->apMagRaw  = psMetadataLookupF32 (&status, row, "AP_MAG_RAW");
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFlux    = psMetadataLookupF32 (&status, row, "AP_FLUX");
+        @>PS1_DV1,>PS1_V3,>PS1_SV1@ source->apFluxErr = psMetadataLookupF32 (&status, row, "AP_FLUX_SIG");
 
         // XXX use these to determine PAR[PM_PAR_I0] if they exist?
-        // XXX add these to PS1_SV1?
-        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
-        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
+        // XXX add these to PS1_SV1?
+        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFlux   = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
+        @>PS1_V2,PS1_SV?,PS1_DV?@ source->psfFluxErr= psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
 
         // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
@@ -424 +1083 @@
 
         @ALL@     source->pixWeightNotBad = psMetadataLookupF32 (&status, row, "PSF_QF");
-        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
+        @>PS1_V2,PS1_SV?,>PS1_DV1@ source->pixWeightNotPoor = psMetadataLookupF32 (&status, row, "PSF_QF_PERFECT");
         @ALL@     source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
         @ALL@     source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
         @ALL@     source->apRadius  = psMetadataLookupF32 (&status, row, "AP_MAG_RADIUS");
-        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
+        @>PS1_V4,>PS1_SV2,>PS1_DV3@ source->apNpixels = psMetadataLookupS32 (&status, row, "AP_NPIX");
 
         // note that some older versions used PSF_PROBABILITY: this was not well defined.
@@ -443 +1102 @@
         @ALL@     source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
 
-        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
-        // we are storing enough information so the output will be consistent with the input
+        // XXX we do not save all of the 3rd and 4th moment parameters. when we load in data,
+        // we are storing enough information so the output will be consistent with the input
         @>PS1_V2,PS1_SV?@ source->moments->Mxxx = +1.0 * psMetadataLookupF32 (&status, row, "MOMENTS_M3C");
         @>PS1_V2,PS1_SV?@ source->moments->Mxxy = 0.0;
@@ -456 +1115 @@
         @>PS1_V2,PS1_SV?@ source->moments->Myyyy = 0.0;
 
-        // Lensing parameters (on read if PS1_V5+)
-        if (haveLensOBJ) {
-          source->lensingOBJ = pmSourceLensingAlloc ();
-          source->lensingOBJ->smear = pmLensingParsAlloc();
-          source->lensingOBJ->shear = pmLensingParsAlloc();
-
-          @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
-          @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
-          @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
-          @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
-          @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
-          @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
-          @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
-          @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
-          @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
-          @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
-        }
-
-        @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM");
-        @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X");
-        @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y");
-
-        if (haveLensPSF) {
-          source->lensingPSF = pmSourceLensingAlloc ();
-          source->lensingPSF->smear = pmLensingParsAlloc();
-          source->lensingPSF->shear = pmLensingParsAlloc();
-
-          @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
-          @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
-          @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
-          @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
-          @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
-          @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
-          @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
-          @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
-          @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
-          @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
-          @>PS1_V4@ source->lensingPSF->e1         = psMetadataLookupF32 (&status, row, "LENS_E1_PSF");
-          @>PS1_V4@ source->lensingPSF->e2         = psMetadataLookupF32 (&status, row, "LENS_E2_PSF");
-        }
+        // Lensing parameters (on read if PS1_V5+)
+        if (haveLensOBJ) {
+          source->lensingOBJ = pmSourceLensingAlloc ();
+          source->lensingOBJ->smear = pmLensingParsAlloc();
+          source->lensingOBJ->shear = pmLensingParsAlloc();
+
+          @>PS1_V4@ source->lensingOBJ->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_OBJ");
+          @>PS1_V4@ source->lensingOBJ->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_OBJ");
+          @>PS1_V4@ source->lensingOBJ->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_OBJ");
+          @>PS1_V4@ source->lensingOBJ->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_OBJ");
+          @>PS1_V4@ source->lensingOBJ->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_OBJ");
+          @>PS1_V4@ source->lensingOBJ->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_OBJ");
+          @>PS1_V4@ source->lensingOBJ->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_OBJ");
+          @>PS1_V4@ source->lensingOBJ->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_OBJ");
+          @>PS1_V4@ source->lensingOBJ->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_OBJ");
+          @>PS1_V4@ source->lensingOBJ->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_OBJ");
+        }
+
+        @>PS1_V4@ source->chipNum = psMetadataLookupS16 (&status, row, "SRC_CHIP_NUM");
+        @>PS1_V4@ source->chipX = psMetadataLookupS16 (&status, row, "SRC_CHIP_X");
+        @>PS1_V4@ source->chipY = psMetadataLookupS16 (&status, row, "SRC_CHIP_Y");
+
+        if (haveLensPSF) {
+          source->lensingPSF = pmSourceLensingAlloc ();
+          source->lensingPSF->smear = pmLensingParsAlloc();
+          source->lensingPSF->shear = pmLensingParsAlloc();
+
+          @>PS1_V4@ source->lensingPSF->smear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SM_PSF");
+          @>PS1_V4@ source->lensingPSF->smear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SM_PSF");
+          @>PS1_V4@ source->lensingPSF->smear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SM_PSF");
+          @>PS1_V4@ source->lensingPSF->smear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SM_PSF");
+          @>PS1_V4@ source->lensingPSF->smear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SM_PSF");
+          @>PS1_V4@ source->lensingPSF->shear->X11 = psMetadataLookupF32 (&status, row, "LENS_X11_SH_PSF");
+          @>PS1_V4@ source->lensingPSF->shear->X12 = psMetadataLookupF32 (&status, row, "LENS_X12_SH_PSF");
+          @>PS1_V4@ source->lensingPSF->shear->X22 = psMetadataLookupF32 (&status, row, "LENS_X22_SH_PSF");
+          @>PS1_V4@ source->lensingPSF->shear->e1  = psMetadataLookupF32 (&status, row, "LENS_E1_SH_PSF");
+          @>PS1_V4@ source->lensingPSF->shear->e2  = psMetadataLookupF32 (&status, row, "LENS_E2_SH_PSF");
+          @>PS1_V4@ source->lensingPSF->e1         = psMetadataLookupF32 (&status, row, "LENS_E1_PSF");
+          @>PS1_V4@ source->lensingPSF->e2         = psMetadataLookupF32 (&status, row, "LENS_E2_PSF");
+        }
 
         @>PS1_V2,PS1_SV?,>PS1_DV1@ source->moments->Mrf         = psMetadataLookupF32 (&status, row, "MOMENTS_R1");
@@ -508 +1167 @@
         @>PS1_V3@ source->skySlope             = psMetadataLookupF32 (&status, row, "SKY_LIMIT_SLOPE");
 
-        @PS1_DV?@  int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
-        @PS1_DV?@  if (nPos) {
-        @PS1_DV?@      source->diffStats = pmSourceDiffStatsAlloc();
-        @PS1_DV?@      source->diffStats->nGood      = nPos;
-        @PS1_DV?@      source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
-        @PS1_DV?@      source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
-        @PS1_DV?@      source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
-        @PS1_DV?@      source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
-        
-        @>PS1_DV1@      source->diffStats->Rp         = psMetadataLookupF32 (&status, row, "DIFF_R_P");
-        @>PS1_DV1@      source->diffStats->SNp        = psMetadataLookupF32 (&status, row, "DIFF_SN_P");
-        @>PS1_DV1@      source->diffStats->Rm         = psMetadataLookupF32 (&status, row, "DIFF_R_M");
-        @>PS1_DV1@      source->diffStats->SNm        = psMetadataLookupF32 (&status, row, "DIFF_SN_M");
-        @PS1_DV?@  }
+        @PS1_DV?@  int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
+        @PS1_DV?@  if (nPos) {
+        @PS1_DV?@      source->diffStats = pmSourceDiffStatsAlloc();
+        @PS1_DV?@      source->diffStats->nGood      = nPos;
+        @PS1_DV?@      source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
+        @PS1_DV?@      source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
+        @PS1_DV?@      source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
+        @PS1_DV?@      source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
+        
+        @>PS1_DV1@      source->diffStats->Rp         = psMetadataLookupF32 (&status, row, "DIFF_R_P");
+        @>PS1_DV1@      source->diffStats->SNp        = psMetadataLookupF32 (&status, row, "DIFF_SN_P");
+        @>PS1_DV1@      source->diffStats->Rm         = psMetadataLookupF32 (&status, row, "DIFF_R_M");
+        @>PS1_DV1@      source->diffStats->SNm        = psMetadataLookupF32 (&status, row, "DIFF_SN_M");
+        @PS1_DV?@  }
 
         @ALL@                      source->mode       = psMetadataLookupU32 (&status, row, "FLAGS");
@@ -533 +1192 @@
 
     return sources;
+}
+
+psArray *pmSourcesRead_CMF_@CMFMODE@ (psFits *fits, psMetadata *header) {
+  // psArray *array = pmSourcesRead_CMF_@CMFMODE@_Old (fits, header);
+      psArray *array = pmSourcesRead_CMF_@CMFMODE@_New (fits, header);
+    return array;
 }
 
@@ -599 +1264 @@
     // write the radial profile apertures to header
     for (int i = 0; i < radMax->n; i++) {
-        sprintf (keyword1, "RMIN_%02d", i);
-        sprintf (keyword2, "RMAX_%02d", i);
-        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
-        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]);
+        sprintf (keyword1, "RMIN_%02d", i);
+        sprintf (keyword2, "RMAX_%02d", i);
+        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword1, PS_META_REPLACE, "min radius for SB profile", radMin->data.F32[i]);
+        psMetadataAddF32 (outhead, PS_LIST_TAIL, keyword2, PS_META_REPLACE, "max radius for SB profile", radMax->data.F32[i]);
     }
 
     // we write out all sources, regardless of quality.  the source flags tell us the state
     for (int i = 0; i < sources->n; i++) {
-        // this is the source associated with this image
+        // this is the source associated with this image
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // skip sources without measurements
@@ -639 +1304 @@
         psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
 
-        float AxialRatio = NAN;
-        float AxialTheta = NAN;
-        pmSourceExtendedPars *extpars = source->extpars;
-        if (extpars) {
-            AxialRatio = extpars->axes.minor / extpars->axes.major;
-            AxialTheta = extpars->axes.theta;
-        }
+        float AxialRatio = NAN;
+        float AxialTheta = NAN;
+        pmSourceExtendedPars *extpars = source->extpars;
+        if (extpars) {
+            AxialRatio = extpars->axes.minor / extpars->axes.major;
+            AxialTheta = extpars->axes.theta;
+        }
         psMetadataAdd (row, PS_LIST_TAIL, "F25_ARATIO",       PS_DATA_F32, "Axial Ratio of radial profile",              AxialRatio);
         psMetadataAdd (row, PS_LIST_TAIL, "F25_THETA",        PS_DATA_F32, "Angle of radial profile ellipse",                  AxialTheta);
@@ -651 +1316 @@
         // Petrosian measurements
         // XXX insert header data: petrosian ref radius, flux ratio
-        // XXX check flags to see if Pet was measured
+        // XXX check flags to see if Pet was measured
         if (doPetrosian) {
-            pmSourceExtendedPars *extpars = source->extpars;
+            pmSourceExtendedPars *extpars = source->extpars;
             if (extpars) {
-                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point 
-                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
-                // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired
-                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
-                if (repairMagErrors) {
-                  // I need to add the kron error in quadrature becasue pet_error ignores the object flux
-                  float Krf  = source->moments ? source->moments->KronFlux : NAN;
-                  float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
-                  if (isfinite (Krf) && isfinite (dKrf)) {
-                    magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf));
-                  } else {
-                    magErr = 1.0 / magErr;
-                  }
-                }
+                // XXX note that this mag is either calibrated or instrumental depending on existence of zero point 
+                float mag = (extpars->petrosianFlux > 0.0) ? -2.5*log10(extpars->petrosianFlux) + magOffset : NAN; // XXX zero point
+                // NOTE EAM 20140806 : PETRO_MAG_ERR was inverted!! this allows for it to be repaired
+                float magErr = (extpars->petrosianFlux > 0.0) ? extpars->petrosianFlux / extpars->petrosianFluxErr : NAN; // XXX zero point
+                if (repairMagErrors) {
+                  // I need to add the kron error in quadrature becasue pet_error ignores the object flux
+                  float Krf  = source->moments ? source->moments->KronFlux : NAN;
+                  float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
+                  if (isfinite (Krf) && isfinite (dKrf)) {
+                    magErr = sqrt(PS_SQR(1.0 / magErr) + PS_SQR(dKrf / Krf));
+                  } else {
+                    magErr = 1.0 / magErr;
+                  }
+                }
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude", mag);
                 psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", magErr);
@@ -705 +1370 @@
         // Flux Annuli (if we have extended source measurements, we have these.  only optionally save them)
         if (doAnnuli) {
-            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
-            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
-            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
-            psVectorInit (radSB, NAN);
-            psVectorInit (radFlux, NAN);
-            psVectorInit (radFill, NAN);
-            if (!source->extpars) goto empty_annuli;
-            if (!source->extpars->radProfile) goto empty_annuli;
-            if (!source->extpars->radProfile->binSB) goto empty_annuli;
-            psAssert (source->extpars->radProfile->binSum, "programming error");
-            psAssert (source->extpars->radProfile->binFill, "programming error");
-            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
-            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
-            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
-
-            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
-            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
-                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
-                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
-                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
-            }
-
-        empty_annuli:
-            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
-            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
-            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
-            psFree (radSB);
-            psFree (radFlux);
-            psFree (radFill);
-        }
-        if (nRow < 0) {
-            nRow = row->list->n;
-        } else {
-            psAssert (nRow == row->list->n, "inconsistent row lengths");
-        }
-        psArrayAdd (table, 100, row);
-        psFree (row);
+            psVector *radSB   = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVector *radFlux = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVector *radFill = psVectorAlloc(radMin->n, PS_TYPE_F32);
+            psVectorInit (radSB, NAN);
+            psVectorInit (radFlux, NAN);
+            psVectorInit (radFill, NAN);
+            if (!source->extpars) goto empty_annuli;
+            if (!source->extpars->radProfile) goto empty_annuli;
+            if (!source->extpars->radProfile->binSB) goto empty_annuli;
+            psAssert (source->extpars->radProfile->binSum, "programming error");
+            psAssert (source->extpars->radProfile->binFill, "programming error");
+            psAssert (source->extpars->radProfile->binSB->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (source->extpars->radProfile->binSum->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (source->extpars->radProfile->binFill->n <= radFlux->n, "inconsistent vector lengths");
+
+            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
+            for (int j = 0; j < source->extpars->radProfile->binSB->n; j++) {
+                radSB->data.F32[j]   = source->extpars->radProfile->binSB->data.F32[j];
+                radFlux->data.F32[j] = source->extpars->radProfile->binSum->data.F32[j];
+                radFill->data.F32[j] = source->extpars->radProfile->binFill->data.F32[j];
+            }
+
+        empty_annuli:
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_SB", PS_DATA_VECTOR, "mean surface brightness annuli", radSB);
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FLUX", PS_DATA_VECTOR, "flux within annuli", radFlux);
+            psMetadataAdd (row, PS_LIST_TAIL, "PROF_FILL", PS_DATA_VECTOR, "fill factor of annuli", radFill);
+            psFree (radSB);
+            psFree (radFlux);
+            psFree (radFill);
+        }
+        if (nRow < 0) {
+            nRow = row->list->n;
+        } else {
+            psAssert (nRow == row->list->n, "inconsistent row lengths");
+        }
+        psArrayAdd (table, 100, row);
+        psFree (row);
     }
     
     if (table->n == 0) {
-        if (!psFitsWriteBlank (fits, outhead, extname)) {
-            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
-            psFree(outhead);
-            psFree(table);
-            return false;
-        }
-        psFree (outhead);
-        psFree (table);
-        return true;
+        if (!psFitsWriteBlank (fits, outhead, extname)) {
+            psError(psErrorCodeLast(), false, "Unable to write empty sources file.");
+            psFree(outhead);
+            psFree(table);
+            return false;
+        }
+        psFree (outhead);
+        psFree (table);
+        return true;
     }
     
     psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
     if (!psFitsWriteTable (fits, outhead, table, extname)) {
-        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
-        psFree (outhead);
-        psFree(table);
-        return false;
+        psError(psErrorCodeLast(), false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
     }
     psFree (outhead);
@@ -885 +1550 @@
         extpars->axes.theta = psMetadataLookupF32(&status, row, "F25_THETA");
 
-        // magErr may have been saved in inverted form
+        // magErr may have been saved in inverted form
         float mag = psMetadataLookupF32(&status, row, "PETRO_MAG");
         float magErr = psMetadataLookupF32(&status, row, "PETRO_MAG_ERR");
         if (isfinite(mag)) {
-          extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
-          if (isfinite(magErr)) {
-            if (repairMagErrors) {
-              // I need to add the kron error in quadrature becasue pet_error ignores the object flux
-              float Krf  = source->moments ? source->moments->KronFlux : NAN;
-              float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
-              if (isfinite (Krf) && isfinite (dKrf)) {
-                magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf));
-              } else {
-                magErr = 1.0 / magErr;
-              }
-              extpars->petrosianFluxErr = extpars->petrosianFlux * magErr;
-            } else {
-              extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
-            }
-          }
-        }
+          extpars->petrosianFlux    = pow(10., (magOffset - mag) / 2.5);
+          if (isfinite(magErr)) {
+            if (repairMagErrors) {
+              // I need to add the kron error in quadrature becasue pet_error ignores the object flux
+              float Krf  = source->moments ? source->moments->KronFlux : NAN;
+              float dKrf = source->moments ? source->moments->KronFluxErr : NAN;
+              if (isfinite (Krf) && isfinite (dKrf)) {
+                magErr = 1.0 / sqrt(PS_SQR(magErr) - PS_SQR(dKrf / Krf));
+              } else {
+                magErr = 1.0 / magErr;
+              }
+              extpars->petrosianFluxErr = extpars->petrosianFlux * magErr;
+            } else {
+              extpars->petrosianFluxErr = extpars->petrosianFlux / magErr;
+            }
+          }
+        }
 
         extpars->petrosianRadius   = psMetadataLookupF32(&status, row, "PETRO_RADIUS");
@@ -973 +1638 @@
     int nParamMax = 0;
     for (int i = 0; i < sources->n; i++) {
-        // this is the source associated with this image
+        // this is the source associated with this image
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         if (source->modelFits == NULL) continue;
@@ -1004 +1669 @@
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // XXX if no model fits are saved, write out modelEXT?
@@ -1018 +1683 @@
 
             // pmSourceExtFitPars *extPars = source->extFitPars->data[j];
-            // assert (extPars);
-
-            // skip models which were not actually fitted
-            // XXX 
-            if (model->flags & badModel) continue;
+            // assert (extPars);
+
+            // skip models which were not actually fitted
+            // XXX 
+            if (model->flags & badModel) continue;
 
             PAR = model->params->data.F32;
@@ -1029 +1694 @@
             yPos = PAR[PM_PAR_YPOS];
 
-            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
-            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
-            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
-            // TRAIL : X,Y are fitted
-            // 
-            
-            // XXX this should be based on what happened, not on the model type
-            if (model->type == modelTypeTrail) {
-                xErr = dPAR[PM_PAR_XPOS];
-                yErr = dPAR[PM_PAR_YPOS];
-            } else {
-                // this is definitely an underestimate since it does not
-                // account for the extent of the source
-                xErr = fwhmMajor * model->magErr / 2.35;
-                yErr = fwhmMinor * model->magErr / 2.35;
-            }
-
-            @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
-            @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0;
-            @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0;
-            @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
-            @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD;
-            @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD;
-            @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD;
-            @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0;
-
-            float kronFlux = source->moments ? source->moments->KronFlux : NAN;
-            float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
+            // for the extended source models, we do not always fit the centroid in the non-linear fitting process
+            // current situation (hard-wired into psphotSourceFits.c:psphotFitPCM,
+            // SERSIC, DEV, EXP : X,Y not fitted (PCM and not PCM)
+            // TRAIL : X,Y are fitted
+            // 
+            
+            // XXX this should be based on what happened, not on the model type
+            if (model->type == modelTypeTrail) {
+                xErr = dPAR[PM_PAR_XPOS];
+                yErr = dPAR[PM_PAR_YPOS];
+            } else {
+                // this is definitely an underestimate since it does not
+                // account for the extent of the source
+                xErr = fwhmMajor * model->magErr / 2.35;
+                yErr = fwhmMinor * model->magErr / 2.35;
+            }
+
+            @>PS1_DV2,>PS1_SV3@ psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
+            @>PS1_DV2,>PS1_SV3@ float posAngle = 0.0;
+            @>PS1_DV2,>PS1_SV3@ float pltScale = 0.0;
+            @>PS1_DV2,>PS1_SV3@ pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+            @>PS1_DV2,>PS1_SV3@ double raPos = ptSky.r*PS_DEG_RAD;
+            @>PS1_DV2,>PS1_SV3@ double decPos = ptSky.d*PS_DEG_RAD;
+            @>PS1_DV2,>PS1_SV3@ posAngle *= PS_DEG_RAD;
+            @>PS1_DV2,>PS1_SV3@ pltScale *= PS_DEG_RAD*3600.0;
+
+            float kronFlux = source->moments ? source->moments->KronFlux : NAN;
+            float kronMag = isfinite(kronFlux) ? -2.5*log10(kronFlux) : NAN;
 
             row = psMetadataAlloc ();
 
-            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
-            // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
+            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
+            // This set of psMetadataAdd Entries marks the "----" "Start of the XFIT segment"
             psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
             psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
@@ -1070 +1735 @@
             @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "RA_EXT",           0, "EXT model ra coordinate",                    raPos);
             @>PS1_DV2,>PS1_SV3@ psMetadataAddF32 (row, PS_LIST_TAIL, "DEC_EXT",          0, "EXT model dec coordinate",                   decPos);
-            @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
-            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
+            @>PS1_DV2@ float instFlux = isfinite(model->mag) ? pow(10.0, -0.4*model->mag) : NAN;
+            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_FLUX",    0, "EXT fit instrumental counts",                instFlux);
 
             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
             psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
 
-            @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
-            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
-            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
-            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
-            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
-
-            // EAM : adding for PV2 outputs:
-            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
+            @>PS1_DV2@ float calMag = isfinite(magOffset) ? model->mag + magOffset : NAN;
+            @>PS1_DV2@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CAL_MAG", PS_DATA_F32, "EXT Magnitude using supplied calibration",   calMag);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_CHISQ",   PS_DATA_F32, "EXT Model Chisq",   model->chisq);
+            @>PS1_DV2,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_NDOF",    PS_DATA_S32, "EXT Model num degrees of freedom",   model->nDOF);
+            @>PS1_SV1,PS1_SV?@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_MODEL_TYPE",    PS_DATA_S32, "type for chosen EXT_MODEL",   source->modelEXT ? source->modelEXT->type : -1);
+
+            // EAM : adding for PV2 outputs:
+            @>PS1_SV1@ psMetadataAdd (row, PS_LIST_TAIL, "EXT_FLAGS", PS_DATA_S16, "model fit flags (pmModelStatus)", model->flags);
 
             @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "POSANGLE",   0, "position angle at source (degrees)",         posAngle);
-            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE",   0, "plate scale at source (arcsec/pixel)",       pltScale);
+            @>PS1_DV2@ psMetadataAddF32 (row, PS_LIST_TAIL, "PLTSCALE",   0, "plate scale at source (arcsec/pixel)",       pltScale);
 
             // psMetadataAddF32 (row, PS_LIST_TAIL, "MOMENTS_XX",       0, "second moment in x",                      extPars->Mxx);
@@ -1103 +1768 @@
 
             // XXX these should be major and minor, not 'x' and 'y'
-            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]); // this is not fitted
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            NAN); // this is not fitted, so error is NAN
-                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
-            } else {
-                if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY])) {
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",     0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",     0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",         0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
-                } else {
-                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_SXX]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SYY]);
-                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_SXY]);
-                }
-            }
+            if (model->type == pmModelClassGetType("PS_MODEL_TRAIL")) {
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", PAR[PM_PAR_LENGTH]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  PAR[PM_PAR_SIGMA]); // this is not fitted
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    PAR[PM_PAR_THETA]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_LENGTH]);
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            NAN); // this is not fitted, so error is NAN
+                psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_THETA]);
+            } else {
+                if (!isfinite(PAR[PM_PAR_SXX]) || !isfinite(PAR[PM_PAR_SYY])  || !isfinite(PAR[PM_PAR_SXY])) {
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",     0, "EXT width (SXX, isnan)", PAR[PM_PAR_SXX]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",     0, "EXT width (SYY, isnan)", PAR[PM_PAR_SYY]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",         0, "EXT angle (SXY, isnan)", PAR[PM_PAR_SXY]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR", 0, "EXT width err (SXX, isnan)", dPAR[PM_PAR_SXX]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR", 0, "EXT width err (SYY, isnan)", dPAR[PM_PAR_SYY]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",     0, "EXT angle err (SXY, isnan)", dPAR[PM_PAR_SXY]);
+                } else {
+                    psEllipseAxes axes = pmPSF_ModelToAxes (PAR, model->class->useReff);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width (major axis), length for trail", axes.major);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width (minor axis), sigma for trail",  axes.minor);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                    axes.theta);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ_ERR",0, "EXT width error (major axis)",            dPAR[PM_PAR_SXX]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN_ERR",0, "EXT width error (minor axis)",            dPAR[PM_PAR_SYY]);
+                    psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA_ERR",    0, "EXT orientation angle (error)",           dPAR[PM_PAR_SXY]);
+                }
+            }
 
             // write out the other generic parameters
@@ -1140 +1805 @@
 
                 snprintf (name, 64, "EXT_PAR_%02d", k);
-               
+               
                 if (k < model->params->n) {
                     psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->params->data.F32[k]);
@@ -1148 +1813 @@
             }
 
-            // optionally, write out the covariance matrix values
-            // XXX do I need to pad this to match the biggest covar matrix?
-            if (false && model->covar) {
-                for (int iy = 0; iy < model->covar->numCols; iy++) {
-                    for (int ix = iy; ix < model->covar->numCols; ix++) {
-                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
-                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
-
-                    }
-                }                   
-            }
+            // optionally, write out the covariance matrix values
+            // XXX do I need to pad this to match the biggest covar matrix?
+            if (false && model->covar) {
+                for (int iy = 0; iy < model->covar->numCols; iy++) {
+                    for (int ix = iy; ix < model->covar->numCols; ix++) {
+                        snprintf (name, 64, "EXT_COVAR_%02d_%02d", iy, ix);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, name, 0, "", model->covar->data.F32[iy][ix]);
+
+                    }
+                }                   
+            }
             psArrayAdd (table, 100, row);
             psFree (row);
@@ -1280 +1945 @@
             PAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07");
             // XXX add an error:
-            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
-        }
-
-        // NOTE: we no longer write out the covariance matrix
-        if (false) {
-            // read the covariance matrix
-            int nparams = model->params->n;
-            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
-            for (int y = 0; y < nparams; y++) {
-                for (int x = 0; x < nparams; x++) {
-                    char name[64];
-                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
-                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
-                }
-            }
-            model->covar = covar;
-        }
-
-        // we are only saving the values stored in dPAR[SXX,etc]
+            // dPAR[7] = psMetadataLookupF32(&status, row, "EXT_PAR_07_");
+        }
+
+        // NOTE: we no longer write out the covariance matrix
+        if (false) {
+            // read the covariance matrix
+            int nparams = model->params->n;
+            psImage *covar = psImageAlloc(nparams, nparams, PS_TYPE_F32);
+            for (int y = 0; y < nparams; y++) {
+                for (int x = 0; x < nparams; x++) {
+                    char name[64];
+                    snprintf(name, 64, "EXT_COVAR_%02d_%02d", y, x);
+                    covar->data.F32[y][x] = psMetadataLookupF32(&status, row, name);
+                }
+            }
+            model->covar = covar;
+        }
+
+        // we are only saving the values stored in dPAR[SXX,etc]
         dPAR[PM_PAR_SXX] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MAJ_ERR");
         dPAR[PM_PAR_SYY] = psMetadataLookupF32(&status, row, "EXT_WIDTH_MIN_ERR");
         dPAR[PM_PAR_SXY] = psMetadataLookupF32(&status, row, "EXT_THETA_ERR");
 
-        // other parameters that we need to read
+        // other parameters that we need to read
         PAR[PM_PAR_SKY] = psMetadataLookupF32(&status, row, "SKY_EXT");
 
@@ -1341 +2006 @@
     // perform full non-linear fits / extended source analysis?
     if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "radial apertures were not measured, skipping\n");
+        return true;
     }
 
@@ -1385 +2050 @@
     for (int i = 0; i < sources->n; i++) {
 
-        // this is the source associated with this image
+        // this is the source associated with this image
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
 
         // skip sources without radial aper measurements (or insufficient)
-        if (source->radialAper == NULL) continue;
+        if (source->radialAper == NULL) continue;
 
         // psAssert (source->radialAper->n == fwhmValues->n, "inconsistent radial aperture set");
 
-        for (int entry = 0; entry < source->radialAper->n; entry++) {
-
-            // choose the convolved EXT model, if available, otherwise the simple one
-            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
-            assert (radialAper);
-
-            if (pmSourcePositionUseMoments(source)) {
-                xPos = source->moments->Mx;
-                yPos = source->moments->My;
-            } else {
-                xPos = source->peak->xf;
-                yPos = source->peak->yf;
-            }
-
-            row = psMetadataAlloc ();
-
-            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
-            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
-            // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment"
-            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
-            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
-            if (fwhmValues) {
-                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
-            } else {
-                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
-            }
-
-            // XXX if we have raw radial apertures, write them out here
-            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
-            psVectorInit (radFlux,    NAN);
-            psVectorInit (radFluxErr, NAN);
-            psVectorInit (radFill,    NAN);
-            if (!radialAper->flux) goto write_annuli;
-            if (!radialAper->fill) goto write_annuli;
-            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
-            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
-
-            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
-            for (int j = 0; j < radialAper->flux->n; j++) {
-                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
-                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
-                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
-                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
-            }
-
-        write_annuli:
-            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
-            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
-            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
-            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
-            psFree (radFlux);
-            psFree (radFluxErr);
-            psFree (radFluxStdev);
-            psFree (radFill);
-
-            psArrayAdd (table, 100, row);
-            psFree (row);
-        }
+        for (int entry = 0; entry < source->radialAper->n; entry++) {
+
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmSourceRadialApertures *radialAper = source->radialAper->data[entry];
+            assert (radialAper);
+
+            if (pmSourcePositionUseMoments(source)) {
+                xPos = source->moments->Mx;
+                yPos = source->moments->My;
+            } else {
+                xPos = source->peak->xf;
+                yPos = source->peak->yf;
+            }
+
+            row = psMetadataAlloc ();
+
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            // the psMetadataAdd entry and the double quotes are used by grep to select the output fields for automatic documentation
+            // This set of psMetadataAdd Entries marks the "----" "Start of the XRAD segment"
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_APER",           0, "Center of aperture measurements",            xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_APER",           0, "Center of aperture measurements",            yPos);
+            if (fwhmValues) {
+                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "FWHM of matched PSF",                    fwhmValues->data.F32[entry]);
+            } else {
+                psMetadataAddF32 (row, PS_LIST_TAIL, "PSF_FWHM",         0, "image is not FWHM-matched",              NAN);
+            }
+
+            // XXX if we have raw radial apertures, write them out here
+            psVector *radFlux      = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxErr   = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFill      = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVector *radFluxStdev = psVectorAlloc(radMax->n, PS_TYPE_F32);
+            psVectorInit (radFlux,    NAN);
+            psVectorInit (radFluxErr, NAN);
+            psVectorInit (radFill,    NAN);
+            if (!radialAper->flux) goto write_annuli;
+            if (!radialAper->fill) goto write_annuli;
+            psAssert (radialAper->flux->n <= radFlux->n, "inconsistent vector lengths");
+            psAssert (radialAper->fill->n <= radFlux->n, "inconsistent vector lengths");
+
+            // copy the data from fluxVal (which is not guaranteed to be the full length) to radFlux
+            for (int j = 0; j < radialAper->flux->n; j++) {
+                radFlux->data.F32[j]      = radialAper->flux->data.F32[j];
+                radFluxErr->data.F32[j]   = radialAper->fluxErr->data.F32[j];
+                radFluxStdev->data.F32[j] = radialAper->fluxStdev->data.F32[j];
+                radFill->data.F32[j]      = radialAper->fill->data.F32[j];
+            }
+
+        write_annuli:
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX",       PS_META_REPLACE, "flux within annuli",       radFlux);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_ERR",   PS_META_REPLACE, "flux error in annuli",     radFluxErr);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FLUX_STDEV", PS_META_REPLACE, "flux standard deviation",  radFluxStdev);
+            psMetadataAddVector (row, PS_LIST_TAIL, "APER_FILL",       PS_META_REPLACE, "fill factor of annuli",    radFill);
+            psFree (radFlux);
+            psFree (radFluxErr);
+            psFree (radFluxStdev);
+            psFree (radFill);
+
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
     }
 
@@ -1581 +2246 @@
     // perform full non-linear fits / extended source analysis?
     if (!psMetadataLookupBool (&status, recipe, "GALAXY_SHAPES")) {
-        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "galaxy shapes were not measured, skipping\n");
+        return true;
     }
 
@@ -1609 +2274 @@
         pmSource *thisSource = sources->data[i];
 
-        // this is the "real" version of this source 
-        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
-
-        // if we did not fit the galaxy model, modelFits will be NULL
+        // this is the "real" version of this source 
+        pmSource *source = thisSource->parent ? thisSource->parent : thisSource;
+
+        // if we did not fit the galaxy model, modelFits will be NULL
         if (source->modelFits == NULL) continue;
 
-        // if we did not fit the galaxy model, galaxyFits will also be NULL
+        // if we did not fit the galaxy model, galaxyFits will also be NULL
         if (source->galaxyFits == NULL) continue;
 

trunk/psModules/src/objects/pmSourceIO_Ghosts.c

@@ -136 +136 @@
     pmFPAview *view = pmFPAviewAlloc (0);
 
-    //We need to check whether we are dealing with an old style ghost_model, or a new style model. Check if the mirror_rad polynomial exists
+    // We need to check whether we are dealing with an old style ghost_model, or a new style model. Check if the mirror_rad polynomial exists
     float mirCheck = 0;
     md = psMetadataLookupMetadata (&status, ghostModel, "GHOST.MIRROR.RAD"); 
@@ -153 +153 @@
     GET_1D_POLY ("GHOST.INNER.MINOR", innerMinor);
 
-
     // select the input astrometry data (also carries the refstars)
     pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
@@ -211 +210 @@
                         double theta0 = atan2(ref->FP->y,ref->FP->x);
 
-                        if(mirCheck) {
+                        // EAM: XXX we should just use the existence of mirrorRad (!= NULL) instead of carrying a different bool
+                        if (mirCheck) {
                             //TdB: first mirror the reference star positions (around the 0,0 pixel) using the radial offset coefficients and the ghost/star angle
                             double ghost_offset_rad = psPolynomial1DEval (mirrorRad, rSrc);
@@ -220 +220 @@
                             ghost->FP->x = ghost_x_fpa_mirror + psPolynomial2DEval(centerX, ghost_x_fpa_mirror, ghost_y_fpa_mirror);
                             ghost->FP->y = ghost_y_fpa_mirror + psPolynomial2DEval(centerY, ghost_x_fpa_mirror, ghost_y_fpa_mirror);
-                        }  
-                        if(!mirCheck) {
+                        } else {
                             //Use the old-style ghost position determination
                             ghost->FP->x = -ref->FP->x + psPolynomial2DEval(centerX, -ref->FP->x, -ref->FP->y);
@@ -303 +302 @@
         psFree (outerMajor);
         psFree (outerMinor);
+        psFree (mirrorRad);
         psFree (ghostModel);
         psFree (view);
@@ -327 +327 @@
     psFree (outerMajor);
     psFree (outerMinor);
+    psFree (mirrorRad);
     psFree (ghostModel);
     psFree (view);