Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/Makefile.am
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/Makefile.am	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/Makefile.am	(revision 37825)
@@ -97,6 +97,6 @@
 # generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf
 psphotFullForce_SOURCES = \
-        psphotFullForce.c             \
-	psphotFullForceArguments.c    \
+        psphotFullForce.c          \
+	psphotFullForceArguments.c \
 	psphotParseCamera.c        \
 	psphotImageLoop.c	   \
@@ -259,5 +259,7 @@
 	psphotEfficiency.c	       \
 	psphotSetNFrames.c	       \
-	psphotSourceMemory.c
+	psphotSourceMemory.c	       \
+	psphotChipParams.c		\
+	psphotGalaxyParams.c
 
 # not currently used
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphot.h
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphot.h	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphot.h	(revision 37825)
@@ -138,9 +138,13 @@
 bool            psphotExtendedSourceFits_Threaded (psThreadJob *job);
 
+bool            psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotGalaxyParams_Threaded (psThreadJob *job);
+
 bool            psphotApResid (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
 bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool            psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf, int index);
 bool            psphotMagnitudes_Threaded (psThreadJob *job);
 
@@ -168,4 +172,5 @@
 // in psphotChoosePSF.c:
 bool            psphotPSFstats (pmReadout *readout, pmPSF *psf);
+bool            psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf);
 bool            psphotMomentsStats (pmReadout *readout, psArray *sources);
 
@@ -232,6 +237,6 @@
 int             psphotSaveImage (psMetadata *header, psImage *image, char *filename);
 bool            psphotDumpConfig (pmConfig *config);
-pmReadout      *psphotSelectBackground (pmConfig *config, const pmFPAview *view);
-pmReadout      *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view);
+pmReadout      *psphotSelectBackground (pmConfig *config, const pmFPAview *view, int index);
+pmReadout      *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view, int index);
 
 // PSF / DBL / EXT evaluation functions
@@ -571,6 +576,10 @@
 bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index);
 
+bool psphotChipParams(pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotChipParams_Threaded(psThreadJob *job);
+
 const char * psphotGetFilerule(const char *baseRule);
 extern bool psphotINpsphotStack;
 
+
 #endif
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChipParams.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChipParams.c	(revision 37825)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChipParams.c	(revision 37825)
@@ -0,0 +1,293 @@
+# include "psphotInternal.h"
+
+bool psphotChipParamsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources);
+bool psphotChipParams_Threaded (psThreadJob *job);
+bool psphotChipParamsSource (pmSource *source, psArray *backmaps);
+
+typedef struct {
+    int         index;
+    short       chipNum;
+    psString    name;
+    float       C0x;
+    float       C1x;
+    float       C2x;
+    float       C0y;
+    float       C1y;
+    float       C2y;
+    int         xMin;
+    int         xMax;
+    int         yMin;
+    int         yMax;
+} psphotChipBackmap ;
+
+psArray * psphotChipParamsParseBackmaps(pmReadout *readout);
+
+bool psphotChipParams (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Chip Parameters ---");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        if (!psphotChipParamsReadout (config, recipe, view, filerule, readout, sources)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
+            return false;
+        }
+    }
+
+    return true;
+}
+
+bool psphotChipParamsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources) {
+
+    bool status = false;
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping chipParams");
+        return true;
+    }
+
+    psTimerStart ("psphot.chipparams");
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+    psArray *backmaps = psphotChipParamsParseBackmaps(readout);
+    if (!backmaps) {
+        psError (PSPHOT_ERR_UNKNOWN, false, "failed to parse chip backmaps  for %s", filerule);
+        return false;
+    }
+
+    // threaded measurement of the source magnitudes
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_CHIP_PARAMS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, backmaps); // sources
+
+// set this to 0 to run without threading
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(backmaps);
+                return false;
+            }
+# else
+	    if (!psphotChipParams_Threaded(job)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		// psFree(AnalysisRegion);
+                psFree(backmaps);
+		return false;
+	    }
+	    psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(backmaps);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(backmaps);
+
+    psLogMsg ("psphot.chip.params", PS_LOG_WARN, "measure chip paramters : %f sec for %ld objects\n", psTimerMark ("psphot.chipparams"), sources->n);
+    return true;
+}
+
+bool psphotChipParams_Threaded (psThreadJob *job) {
+
+    psArray *sources                = job->args->data[0];
+    psArray *backmaps               = job->args->data[1];
+
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+	if (!source->peak) continue; // XXX how can we have a peak-less source?
+
+        psphotChipParamsSource(source, backmaps);
+    }
+
+    return true;
+}
+
+
+bool psphotChipParamsSource (pmSource *source, psArray *backmaps) {
+
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
+
+    float warp_x = source->peak->xf;
+    float warp_y = source->peak->yf;
+
+    int chipForSource = -1;
+    for (int i=0; i < backmaps->n; i++) {
+        psphotChipBackmap *thisChip = backmaps->data[i];
+
+        float chip_x = thisChip->C0x + thisChip->C1x * warp_x + thisChip->C2x * warp_y;
+        float chip_y = thisChip->C0y + thisChip->C1y * warp_x + thisChip->C2y * warp_y;
+        if (chip_x >= thisChip->xMin && chip_x <= thisChip->xMax &&
+            chip_y >= thisChip->yMin && chip_y <= thisChip->yMax) {
+
+            chipForSource = thisChip->chipNum;
+            if (0) {
+                printf("source %8d on chip: %s %2d at %6.1f, %6.1f\n", source->chipNum, thisChip->name, 
+                    chipForSource, chip_x, chip_y);
+            }
+            source->chipNum = chipForSource;
+            source->chipX = chip_x;
+            source->chipY = chip_y;
+            // Since it is on this chip we are done for this source
+            break;
+        }
+    }
+    if (chipForSource < 0) return false;
+
+    return true;
+}
+
+
+psArray * psphotChipParamsParseBackmaps(pmReadout *readout) {
+
+    // XXX: I'm assuming the structure of the FPA here. Find a more usual way
+    // to find header
+    psMetadata *header = readout->parent->parent->parent->hdu->header;
+    if (!header) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find header for readout.");
+        return NULL;
+    }
+
+    psArray *backmaps = psArrayAllocEmpty(6);
+
+    // find the keywords that define the warp to chip transformations stored by
+    // pswarpUpdateMetadata()
+
+    for (int i = 0; ; i++) {
+        bool mdok;
+        char keyword[16];
+
+        sprintf(keyword, "SRC_%04d", i);
+        psString chipName = psMetadataLookupStr(&mdok, header, keyword);
+        if (!chipName) {
+            break;
+        }
+        psphotChipBackmap *thisChip = psAlloc(sizeof(psphotChipBackmap));
+        psArrayAdd(backmaps, 1, thisChip);
+        psFree(thisChip);
+
+        thisChip->index = i;
+        thisChip->name = chipName; // not incrementing the reference this gets freed before header
+
+        // XXX: If chipName begins with XY assume gpc1 style class_id and use number as ID
+        if (!strncmp(chipName, "XY", 2)) {
+            thisChip->chipNum = atoi(chipName+2);
+        } else {
+            // otherwise just use the index
+            thisChip->chipNum = thisChip->index;
+        }
+
+        sprintf(keyword, "SEC_%04d", i);
+        psString chipRegionStr = psMetadataLookupStr(&mdok, header, keyword);
+        if (!chipRegionStr) {
+            psError(PS_ERR_UNKNOWN, true, "failed to find required keyword %s", keyword);
+            psFree(backmaps);
+            return NULL;
+        }
+        int nvals = sscanf(chipRegionStr, "[%d:%d,%d:%d]", &thisChip->xMin, &thisChip->xMax, &thisChip->yMin, &thisChip->yMax);
+        if (nvals != 4) {
+            psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, chipRegionStr);
+            psFree(backmaps);
+            return NULL;
+        }
+        // convert from fits extent convention to zero offset array
+        #define COORD_SLACK 2
+        thisChip->xMin -= COORD_SLACK; thisChip->xMax += COORD_SLACK; 
+        thisChip->yMin -= COORD_SLACK; thisChip->yMax += COORD_SLACK;
+
+        sprintf(keyword, "MPX_%04d", i);
+        psString mpx = psMetadataLookupStr(&mdok, header, keyword);
+        if (!mpx) {
+            psError(PS_ERR_UNKNOWN, true, "failed to find %s", keyword);
+            psFree(backmaps);
+            return NULL;
+        }
+        nvals = sscanf(mpx, "[%f,%f,%f]", &thisChip->C0x, &thisChip->C1x, &thisChip->C2x);
+        if (nvals != 3) {
+            psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, mpx);
+            psFree(backmaps);
+            return NULL;
+        }
+        sprintf(keyword, "MPY_%04d", i);
+        psString mpy = psMetadataLookupStr(&mdok, header, keyword);
+        if (!mpy) {
+            psError(PS_ERR_UNKNOWN, true, "failed to find %s", keyword);
+            psFree(backmaps);
+            return NULL;
+        }
+        nvals = sscanf(mpy, "[%f,%f,%f]", &thisChip->C0y, &thisChip->C1y, &thisChip->C2y);
+        if (nvals != 3) {
+            psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, mpx);
+            psFree(backmaps);
+            return NULL;
+        }
+    }
+
+    if (!backmaps->n) {
+        psFree(backmaps);
+        return NULL;
+    }
+
+    return backmaps;
+}
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChooseAnalysisOptions.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChooseAnalysisOptions.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChooseAnalysisOptions.c	(revision 37825)
@@ -606,4 +606,10 @@
     psLogMsg ("psphot.options", PS_LOG_WARN, "choose analysis options for %ld objects: %f sec\n", objects->n, psTimerMark ("psphot.options"));
 
+    psFree(exptime);
+    psFree(zeropt);
+    psFree(chips);
+    psFree(inputFilters);
+    psFree(petroFluxLim);
+    psFree(extFitFluxLim);
     psFree(toGal);
 
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChoosePSF.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChoosePSF.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotChoosePSF.c	(revision 37825)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf);
 
 // generate a PSF model for inputs without PSF models already loaded
@@ -388,5 +389,7 @@
     psFree (models);
 
-    if (!psphotPSFstats (readout, psf)) {
+    // if (!psphotPSFstats (readout, psf)) {
+
+    if (!psphotPSFstatsSources (readout, sources, psf)) {
         psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
         psFree(options);
@@ -466,8 +469,8 @@
 
 	    if (modelPSF->params->n > 7) {
-	      psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
+		psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
 	    }
 	    if (modelPSF->params->n > 8) {
-	      psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
+		psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
 	    }
             psFree (modelPSF);
@@ -522,19 +525,172 @@
 
     if (psfExtra1->n) {
-      if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+	if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+	    goto escape;
+	}
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
+	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    if (psfExtra2->n) {
+	if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
+	    goto escape;
+	}
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
+	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
+    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
+
+    // psLogMsg ("psphot", PS_LOG_DETAIL, "PSF angle: %f, nstars: %d", axes.theta, psf->nPSFstars);
+
+    char *psfModelName = pmModelClassGetName(psf->type);
+    psMetadataAddStr(readout->analysis,  PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName);
+    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK",   PS_META_REPLACE, "Valid PSF Model?", true);
+
+    int nParams = pmModelClassParameterCount(psf->type);
+    psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR",   PS_META_REPLACE, "Number of PSF parameters", nParams);
+
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
+    psFree (stats);
+    return true;
+
+escape:
+    psFree (fwhmMajor);
+    psFree (fwhmMinor);
+    psFree (psfExtra1);
+    psFree (psfExtra2);
+    psFree (stats);
+    return false;
+}
+
+// measure average parameters of the PSF model
+bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf) {
+
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_PTR_NON_NULL(psf, false);
+
+    psImage *image = readout->image;
+    PS_ASSERT_PTR_NON_NULL(image, false);
+
+    // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc:
+    psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32);
+    psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+
+	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+
+	float xc = source->peak->xf;
+	float yc = source->peak->yf;
+
+	// create modelPSF from this model
+	pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0);
+	if (!modelPSF) {
+	    fprintf (stderr, "?");
+	    continue;
+	}
+
+	// get the model full-width at half-max
+	float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5);
+
+	// XXX make sure this is consistent with the re-definition of PM_PAR_SXX
+	shape.sx  = modelPSF->params->data.F32[PM_PAR_SXX];
+	shape.sy  = modelPSF->params->data.F32[PM_PAR_SYY];
+	shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY];
+	axes = psEllipseShapeToAxes (shape, 20.0);
+
+	float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+	if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) {
+	    fprintf (stderr, "!");
+	    psFree (modelPSF);
+	    continue;
+	}
+	psVectorAppend (fwhmMajor, FWHM_MAJOR);
+	psVectorAppend (fwhmMinor, FWHM_MINOR);
+
+	if (modelPSF->params->n > 7) {
+	    psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]);
+	}
+	if (modelPSF->params->n > 8) {
+	    psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]);
+	}
+	psFree (modelPSF);
+
+	// float fwhmtest = pmPSFtoFWHM(psf, xc, yc);
+	// fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest);
+    }
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
+    if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR");
 	goto escape;
-      }
-      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
-      psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+    }
+
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
+
+    float fwhmMaj = stats->sampleMean;  // FWHM on major axis
+
+    if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR");
+	goto escape;
+    }
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
+    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
+
+    float fwhmMin = stats->sampleMean;  // FWHM on minor axis
+
+    if (readout->parent) {
+
+	// we now have 2 possible measurements of the seeing : the PSF model based version and
+	// the Moments based value we need to define a definitive "CHIP.SEEING" value.  This
+	// value is used by the PSF-matching programs (ppSub and ppStack).  The moments-based
+	// value is probably more representative of the value needed by those tools (it is also
+	// more stable?)
+	bool status = false;
+	float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+	float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1");
+
+	pmChip *chip = readout->parent->parent; // Parent chip
+	psAssert(chip, "Cell should be attached to a chip.");
+	psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip
+	item->data.F32 = 0.5 * (fwhmMajorMoments + fwhmMinorMoments);
+
+	psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments);
+    }
+
+    if (psfExtra1->n) {
+	if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1");
+	    goto escape;
+	}
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean);
+	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     }
 
     if (psfExtra2->n) {
-      if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
-	goto escape;
-      }
-      psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
-      psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
+	if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) {
+	    psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2");
+	    goto escape;
+	}
+	psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean);
+	psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev);
     }
 
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotCleanup.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotCleanup.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotCleanup.c	(revision 37825)
@@ -25,5 +25,5 @@
     // don't display memory leaks unless trace is on. 
     // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psphot");
-    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psphot");
+    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks2 (0, NULL, stdout, false, 500), "psphot");
 #endif
     return;
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotEfficiency.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotEfficiency.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotEfficiency.c	(revision 37825)
@@ -502,5 +502,5 @@
 
     // measure the magnitudes and fluxes for the sources
-    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf)) {
+    if (!psphotMagnitudesReadout(config, recipe, view, readout, fakeSourcesAll, psf, index)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to measure magnitudes of fake sources.");
         psFree(fakeSources);
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceAnalysis.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceAnalysis.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceAnalysis.c	(revision 37825)
@@ -38,11 +38,11 @@
 
 /*** for the moment, this test code : it is not thread safe ***/
-int    Nall = 0;
-int  Nskip1 = 0;
-int  Nskip2 = 0;
-int  Nskip3 = 0;
-int  Nskip4 = 0;
-int  Nskip5 = 0;
-int  Nskip6 = 0;
+static int    Nall = 0;
+static int  Nskip1 = 0;
+static int  Nskip2 = 0;
+static int  Nskip3 = 0;
+static int  Nskip4 = 0;
+static int  Nskip5 = 0;
+static int  Nskip6 = 0;
 
 # define SKIP(VALUE) { VALUE++; continue; }
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceFits.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceFits.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotExtendedSourceFits.c	(revision 37825)
@@ -629,5 +629,5 @@
 	  pmSourceCacheModel (source, maskVal);
 
-# if (PS_TRACE_ON)
+# if (0 && PS_TRACE_ON)
 	  pmModel *model = source->modelFits->data[0];
 	  int flags = 0xffffffff;
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFitSourcesLinear.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFitSourcesLinear.c	(revision 37825)
@@ -263,4 +263,10 @@
 	    }
 	}
+        if (modelSum == 0 ||  maskedSum == 0) {
+            // no pixels at all drop the modelPSF (probably only an issue for external sources
+            // which are created from a fake peak (full force))
+            psFree(source->modelPSF);
+            continue;
+        }
 	if (modelSum < 0.5) continue; // skip sources with no model constraint (somewhat arbitrary limit)
 	if (modelSum < 0.8) {
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFullForceReadout.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFullForceReadout.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotFullForceReadout.c	(revision 37825)
@@ -106,4 +106,6 @@
     psphotKronFlux (config, view, filerule);
 
+    psphotChipParams (config, view, filerule);
+
     // measure petro fluxes
     psphotPetroFlux (config, view, filerule);
@@ -114,4 +116,10 @@
     // measure galaxy shapes
     psphotGalaxyShape (config, view, filerule);
+
+    // measure aperture photometry corrections
+    if (!psphotApResid (config, view, filerule)) {
+        psLogMsg ("psphot", 3, "failed on psphotApResid");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
 
     // calculate source magnitudes (psf mag and ap mag)
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotGalaxyParams.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotGalaxyParams.c	(revision 37825)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotGalaxyParams.c	(revision 37825)
@@ -0,0 +1,668 @@
+# include "psphotInternal.h"
+
+# define PRINT_STUFF 0
+# define USE_FOOTPRINTS 1
+
+// measure some properties on the exended objects
+
+bool psphotGalaxyParams (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Galaxy Parameters ---");
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    // add this 
+    bool doGalaxyParams = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_GALAXY_PARAMS");
+    bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    bool doAnnuli    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+    bool doExtModels = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS");
+
+    // All of these must be true to measure galaxy parameters
+    if (!doGalaxyParams || !doPetrosian || !doAnnuli || !doExtModels) {
+	psLogMsg ("psphot", PS_LOG_INFO, "skipping Galaxy parameter source measurements\n");
+	return true;
+    }
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (!psphotGalaxyParamsReadout (config, view, filerule, i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy parameter measurements for %s entry %d", filerule, i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+// # define COUNT_SKIPS
+# ifdef COUNT_SKIPS
+/*** for the moment, this test code : it is not thread safe ***/
+static int  Nskip1 = 0;
+static int  Nskip2 = 0;
+static int  Nskip3 = 0;
+static int  Nskip4 = 0;
+static int  Nskip5 = 0;
+static int  Nskip6 = 0;
+static int  Nskip7 = 0;
+static int  Nskip8 = 0;
+# define SKIP(VALUE) { VALUE++; continue; }
+# else
+# define SKIP(VALUE) {continue;}
+# endif //notdef
+
+// aperture-like measurements for extended sources
+bool psphotGalaxyParamsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    bool status;
+    int Next = 0;
+    int Npetro = 0;
+    int Nannuli = 0;
+
+    psTimerStart ("psphot.galaxy");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+	return true;
+    }
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // get the sky noise from the background analysis; if this is missing, get the user-supplied value
+    float skynoise = psMetadataLookupF32 (&status, readout->analysis, "SKY_STDEV");
+    if (!status) {
+	skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
+	psWarning ("failed to get sky noise level from background analysis; defaulting to user supplied value of %f\n", skynoise);
+    }
+
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+        fitNsigmaConv = 5.0;
+    }
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+
+    psImage *footprintImage = NULL;
+#if USE_FOOTPRINTS
+    footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
+    psImageInit(footprintImage, 0);
+    pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, false);
+#endif
+
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortByFlux);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_PARAMS");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, recipe);
+
+            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32);
+            psArrayAdd(job->args, 1, footprintImage);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
+
+
+// set this to 0 to run without threading
+# if (0)	    
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		psFree(AnalysisRegion);
+                return false;
+            } 
+# else
+	    if (!psphotGalaxyParams_Threaded(job)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		psFree(AnalysisRegion);
+		return false;
+	    }
+	    psScalar *scalar = NULL;
+	    scalar = job->args->data[8];
+	    Next += scalar->data.S32;
+	    scalar = job->args->data[9];
+	    Npetro += scalar->data.S32;
+	    scalar = job->args->data[10];
+	    Nannuli += scalar->data.S32;
+	    psFree(job);
+# endif
+	}
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+	    psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+		psScalar *scalar = NULL;
+		scalar = job->args->data[8];
+		Next += scalar->data.S32;
+		scalar = job->args->data[9];
+		Npetro += scalar->data.S32;
+		scalar = job->args->data[10];
+		Nannuli += scalar->data.S32;
+            }
+            psFree(job);
+	}
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+    psFree(footprintImage);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d petrosian\n", Npetro);
+    psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
+
+// # if (PS_TRACE_ON)
+# ifdef COUNT_SKIPS
+    fprintf (stderr, "ext analysis skipped @ 1  : %d\n", Nskip1);
+    fprintf (stderr, "ext analysis skipped @ 2  : %d\n", Nskip2);
+    fprintf (stderr, "ext analysis skipped @ 3  : %d\n", Nskip3);
+    fprintf (stderr, "ext analysis skipped @ 4  : %d\n", Nskip4);
+    fprintf (stderr, "ext analysis skipped @ 5  : %d\n", Nskip5);
+    fprintf (stderr, "ext analysis skipped @ 6  : %d\n", Nskip6);
+    fprintf (stderr, "ext analysis skipped @ 7  : %d\n", Nskip7);
+    fprintf (stderr, "ext analysis skipped @ 8  : %d\n", Nskip8);
+#endif
+
+    // psphotRadialProfileShowSkips ();
+
+    // psphotVisualShowResidualImage (readout, false);
+
+    // psphotVisualShowPetrosians (sources);
+
+    // XXX TEST: SAVE IMAGE:
+    // psphotSaveImage (NULL, readout->image, "resid.fits");
+
+    return true;
+}
+
+bool psphotGalaxyParams_Threaded (psThreadJob *job) {
+
+    // bool status;
+
+    int Next = 0;
+    int Npetro = 0;
+    int Nannuli = 0;
+    int Ngood = 0;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psRegion *region        = job->args->data[2];
+    psMetadata *recipe      = job->args->data[3];
+
+#if (!USE_FOOTPRINTS)
+    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
+#endif
+    float fitNsigmaConv     = PS_SCALAR_VALUE(job->args->data[5],F32);
+    float psfSize           = PS_SCALAR_VALUE(job->args->data[6],F32);
+    psImage *footprintImage = job->args->data[7];
+
+    bool status;
+//    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // don't try and use bad models
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+//    badModel |= PM_MODEL_STATUS_NONCONVERGE;  this causes most objects to get unmeasured
+    badModel |= PM_MODEL_STATUS_BADARGS;
+    badModel |= PM_MODEL_STATUS_OFFIMAGE;
+    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+    badModel |= PM_MODEL_PCM_FAIL_GUESS;
+
+    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
+
+    psVector *bgPixels = NULL;
+    psVector *bgPixelsFlip = NULL;
+#define dontMAKE_MARK_IMAGE 1
+#ifdef MAKE_MARK_IMAGE
+    psImage *markImage = NULL;
+#endif
+
+    #define MARK_SOURCE 0
+    #define MARK_MASK 1
+    #define MARK_BG 2
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+
+#define dontDUMP_IMAGES 1
+#ifdef DUMP_IMAGES
+    char *filebase = "/data/ippc18.0/home/bills/debug/psphot/images";
+    char filename[80];
+    psMetadata *dummyHeader = psMetadataAlloc();
+#endif
+
+    // choose the sources of interest
+    for (int i = 0; i < sources->n; i++) {
+
+	pmSource *source = sources->data[i];
+
+	// Nall ++;
+
+	bool testObject = FALSE;
+	// testObject = testObject || ((fabs(source->peak->xf - 5788) < 10) && (fabs(source->peak->yf - 5711) < 10));
+	// testObject = testObject || ((fabs(source->peak->xf -   96) < 10) && (fabs(source->peak->yf - 5550) < 10));
+	// testObject = testObject || ((fabs(source->peak->xf -  300) < 10) && (fabs(source->peak->yf - 5988) < 10));
+	if (testObject) {
+	  fprintf (stderr, "test object: %f, %f\n", source->peak->xf, source->peak->yf);
+	  testObject = TRUE;
+	}
+
+	// rules for measuring petrosian parameters for specific objects are set in
+	// psphotChooseAnalysisOptions.c
+	if (!(source->tmpFlags & PM_SOURCE_TMPF_PETRO)) SKIP (Nskip1);
+
+	// limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
+	if (source->peak->x < region->x0) SKIP (Nskip2);
+	if (source->peak->y < region->y0) SKIP (Nskip3);
+	if (source->peak->x > region->x1) SKIP (Nskip4);
+	if (source->peak->y > region->y1) SKIP (Nskip5);
+
+        // Should we skip objects that do not have an extended model?
+        if (!source->extpars) SKIP (Nskip6);
+        if (!source->modelFits) SKIP (Nskip7);
+
+        // Flux cut to reduce sources to those that match SDSS
+        // if (source->moments->KronFlux < 2e5) continue;
+
+        pmModel *sersicModel = NULL;
+        for (int j = 0; j < source->modelFits->n; j++) {
+            pmModel *model = source->modelFits->data[j];
+            if (model->type == sersicModelType) {
+                sersicModel = model;
+                break;
+            }
+        }
+        if (!sersicModel) SKIP (Nskip8);
+        if (sersicModel->flags & badModel) SKIP (Nskip8);
+	Next ++;
+
+        // if sersicModel is not modelEXT make it so
+        pmSourceMode saveMode = source->mode;
+        pmSourceType saveType = source->type;
+        bool isExtended = saveType == PM_SOURCE_TYPE_EXTENDED;
+        pmModel *extModel = source->modelEXT;
+
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+
+        if (!isExtended || extModel != sersicModel) {
+            // XXX: I'm not doing this right. For now just look at objects for which
+            source->modelEXT = psMemIncrRefCounter (sersicModel);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            source->mode |= PM_SOURCE_MODE_EXTMODEL;
+            // do we need this?
+            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+#ifdef notdef
+#endif
+        }
+
+
+        psF32 *PAR = sersicModel->params->data.F32;
+        psF32 Xo = PAR[PM_PAR_XPOS];
+        psF32 Yo = PAR[PM_PAR_YPOS];
+        psF32 Io = PAR[PM_PAR_I0];
+        // For now: set the half light radius to the major axis of the model
+        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff);
+        source->extpars->ghalfLightRadius = axes.major;
+        // and we will integrate out to 2 * half light radius
+        psF32 R = 2 * axes.major;
+        psF32 R2 = PS_SQR(R);
+
+	float radius = source->peak->xf - source->pixels->col0;
+	radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
+	radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
+	radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
+        // Make image at least 10 * half light radius
+        radius = PS_MAX (radius, 5*axes.major);
+	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius);
+
+        // RedefinePixels frees modelFlux
+        if (source->modelEXT->isPCM) {
+            pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
+        } else {
+            pmSourceCacheModel (source, maskVal);
+        }
+
+#ifdef DUMP_IMAGES
+        sprintf(filename, "%s/%04d-pix.fits", filebase, source->id);
+
+        psphotSaveImage(dummyHeader, source->pixels, filename);
+#endif
+
+
+        // Create the markImage and accumlate the Flux
+
+        // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
+        // xCM, yCM from pixel coords to pixel index here.
+        psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
+        psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+        psF32 **vPix = source->pixels->data.F32;
+        psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+
+        bgPixels = psVectorRecycle(bgPixels, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
+        bgPixels->n = 0;
+        bgPixelsFlip = psVectorRecycle(bgPixelsFlip, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
+        bgPixelsFlip->n = 0;
+
+#ifdef MAKE_MARK_IMAGE
+        // The markImage may no longer be needed
+        markImage = psImageRecycle(markImage, source->pixels->numCols, source->pixels->numRows, PS_TYPE_U8);
+
+        // assume all pixels are masked to start
+        psImageInit(markImage, MARK_MASK);
+        psU8 **vMark = markImage->data.U8;
+        #define SET_MARK(_val) vMark[row][col] = _val
+#else
+        #define SET_MARK(_val)
+#endif
+
+#if (USE_FOOTPRINTS)
+        psS32 **vFootprint = footprintImage->data.S32;
+#else
+        // XXX use recipe value instead of hard coded 1.5
+        psF32 threshold = skynoise * 3;
+#endif
+
+        psF32 sumI = 0;
+        psS32 row0 = source->pixels->row0;
+        psS32 col0 = source->pixels->col0;
+        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+            psF32 yDiff = row - yCM;
+            // y coordinate of mirror pixel
+            int yFlip = yCM - yDiff;
+            if (yFlip < 0) continue;
+            if (yFlip >= source->pixels->numRows) continue;
+            psS32 yImage = row + row0;
+            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+                // check mask and value for this pixel
+                if (vMsk && (vMsk[row][col] & maskVal)) continue;
+                psF32 pixel = vPix[row][col];
+                if (isnan(pixel)) continue;
+                // valid pixel
+
+                psS32 xImage = col + col0;
+
+                psF32 xDiff = col - xCM;
+
+                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
+                if (r2 <= R2) {
+                    sumI += pixel;
+                }
+
+#if (USE_FOOTPRINTS)
+		// this is a pixel in the object
+                if (vFootprint[yImage][xImage] == source->peak->footprint->id) {
+                    SET_MARK(MARK_SOURCE);
+                    continue;
+		}
+		// this is a pixel in another object
+                if (vFootprint[yImage][xImage] > 0.0) continue;
+#else
+		// this is a pixel in an object (this or the other)
+                if (fabs(pixel) >= threshold) {
+		  continue
+                }
+#endif
+                // we have a background pixel.
+                // Check whether it's mirror is a background pixel as well.
+
+                // x coordinate of mirror pixel
+                int xFlip = xCM - xDiff;
+                if (xFlip < 0) continue;
+                if (xFlip >= source->pixels->numCols) continue;
+                // check mask and value for mirror pixel
+                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+                psF32 mirrorPixel = vPix[yFlip][xFlip];
+                if (isnan(mirrorPixel)) continue;
+
+#if (USE_FOOTPRINTS)
+                if (vFootprint[yFlip + row0][xFlip + col0] == 0) {
+#else
+                if (fabs(mirrorPixel) < threshold) {
+#endif
+                    // This is a good background pixel with a good mirror Save the values 
+                    SET_MARK(MARK_BG);
+                    psVectorAppend(bgPixels, pixel);
+                    psVectorAppend(bgPixelsFlip, mirrorPixel);
+                }
+            }
+        }
+#ifdef DUMP_IMAGES
+        sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id);
+
+        psphotSaveImage(dummyHeader, markImage, filename);
+#endif
+
+        // subtract model to create the residual image
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
+            // Hmm something has gone wrong
+            fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id);
+            goto restoreModel;
+        }
+#ifdef DUMP_IMAGES
+        sprintf(filename, "%s/%04d-res.fits", filebase, source->id);
+
+        psphotSaveImage(dummyHeader, source->pixels, filename);
+#endif
+
+        // don't know if these pointers could change in the background subtraction so grab them again
+        vPix = source->pixels->data.F32;
+        vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
+        psF32 **vVar = source->variance->data.F32;
+        psF32 **vModel = source->modelFlux->data.F32;
+        // Accumulate Rt and Ra
+        psF32 sumRt = 0;
+        psF32 sumRa = 0;
+        int numPixels = 0;
+        psF32 sumB = 0;
+        psF32 sumModel = 0;
+        int numB = 0;
+        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+            psF32 yDiff = row - yCM;
+            if (fabs(yDiff) > R) continue;
+            // y coordinate of mirror pixel
+            int yFlip = yCM - yDiff;
+            if (yFlip < 0) continue;
+            if (yFlip >= source->pixels->numRows) continue;
+            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+                // check mask and value for this pixel
+                if (vMsk && (vMsk[row][col] & maskVal)) continue;
+                psF32 pixel = vPix[row][col];
+                if (isnan(pixel)) continue;
+
+                psF32 xDiff = col - xCM;
+                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
+                if (r2 > R2) continue;
+
+                // Measurement of bumpiness parameter excludes the pixels
+                // within 2 pixels of the center
+                if (r2 > 4) {
+                    psF32 variance = vVar[row][col];
+                    if (isfinite(variance)) {
+                        psF32 modelPixel = vModel[row][col];
+                        if (!isfinite(modelPixel)) {
+                            // XXX: Can this happen in a good model?
+                            fprintf(stderr, "non-finite modelPixel %4d (%4d, %4d)\n", source->id, col, row);
+                        }
+                        numB++;
+                        sumB += PS_SQR(pixel) - variance;
+                        sumModel += Io * modelPixel;
+                    }
+                }
+
+                // x coordinate of mirror pixel
+                int xFlip = xCM - xDiff;
+                if (xFlip < 0) continue;
+                if (xFlip >= source->pixels->numCols) continue;
+                // check mask and value for mirror pixel
+                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
+                psF32 mirrorPixel = vPix[yFlip][xFlip];
+                if (isnan(mirrorPixel)) continue;
+
+                sumRt += fabs(pixel + mirrorPixel);
+                sumRa += fabs(pixel - mirrorPixel);
+                ++numPixels;
+            }
+        }
+
+        if (numB > 0) {
+            source->extpars->gbumpy = sqrt(sumB/numB)  / (sumModel / numB);
+        }
+
+        if (bgPixels->n < numPixels) {
+            fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels);
+            goto restoreModel;
+        }
+        Ngood++;
+
+        // and calculate Bt and Ba
+        psF32 sumBt = 0;
+        psF32 sumBa = 0;
+        psF32 *vB = bgPixels->data.F32;
+        psF32 *vBFlip = bgPixelsFlip->data.F32;
+        long    length = bgPixels->n;
+        // Randomly select numPixels background pixels
+        for (int n = 0; n < numPixels; n++) {
+            double rand = psRandomUniform(rng);
+            long index = floor(rand * length);
+            if (index >= length) {
+                fprintf(stderr, "random index %ld is larger than vector length %ld\n",
+                    index, length);
+                continue;
+            }
+            sumBt += fabs(vB[index] + vBFlip[index]);
+            sumBa += fabs(vB[index] - vBFlip[index]);
+        }
+
+        // finally save the results
+        source->extpars->gRT = 0.5 * (sumRt - sumBt)/sumI;
+        source->extpars->gRA = 0.5 * (sumRa - sumBa)/sumI;
+        source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA;
+
+#if (PRINT_STUFF)
+        if (Ngood % 40 == 1) {
+	  fprintf (stderr, "   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
+        }
+        fprintf (stderr, "%5d %5.1f %5.1f : %5.2e %4.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
+		 source->id, source->peak->xf, source->peak->yf, 
+		 sumI, source->extpars->ghalfLightRadius, 
+		 radius, 
+		 source->extpars->gRT, 
+		 source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n);
+#endif
+
+        // Now put things back the way that they were.
+restoreModel:
+
+	if (testObject) {
+	  fprintf (stderr, "done with %d\n", source->id);
+	}
+
+# define USE_BEST_FIT TRUE
+	// XXX TEST: Keep the sersic model regardless
+        if (USE_BEST_FIT && ((extModel != sersicModel) || !isExtended)) {
+            // add the sersic model back in
+            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+            }
+
+            // replace the extended model
+            if (extModel != sersicModel) {
+                psMemDecrRefCounter (sersicModel);
+                source->modelEXT = extModel;
+            }
+            source->mode = saveMode;
+            source->type = saveType;
+            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        }
+    }
+
+#ifdef MAKE_MARK_IMAGE
+    psFree(markImage);
+#endif
+#ifdef DUMP_IMAGES
+    psFree(dummyHeader);
+    // psphotSaveImage raises an error
+    psErrorClear();
+#endif
+    psFree(bgPixels)
+    psFree(bgPixelsFlip)
+    psFree(rng);
+
+    psScalar *scalar = NULL;
+
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[8];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[9];
+    scalar->data.S32 = Npetro;
+
+    scalar = job->args->data[10];
+    scalar->data.S32 = Nannuli;
+
+    return true;
+}
+
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotMagnitudes.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotMagnitudes.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotMagnitudes.c	(revision 37825)
@@ -38,5 +38,5 @@
         psAssert (psf, "missing psf?");
 
-        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
+        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -46,5 +46,5 @@
 }
 
-bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
+bool psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
 
     bool status = false;
@@ -65,6 +65,6 @@
 
     // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
-    pmReadout *backModel = psphotSelectBackground (config, view);
-    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
+    pmReadout *backModel = psphotSelectBackground (config, view, index);
+    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view, index);
 
     // bit-masks to test for good/bad pixels
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotOutput.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotOutput.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotOutput.c	(revision 37825)
@@ -29,10 +29,7 @@
 }
 
-pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view) {
-
-    bool status;
-    
-
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL"));
+pmReadout *psphotSelectBackground (pmConfig *config, const pmFPAview *view, int index) {
+
+    pmFPAfile *file = pmFPAfileSelectSingle (config->files, psphotGetFilerule("PSPHOT.BACKMDL"), index);
     if (!file) return NULL;
 
@@ -41,9 +38,7 @@
 }
 
-pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view) {
-
-    bool status;
-
-    pmFPAfile *file = psMetadataLookupPtr (&status, config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"));
+pmReadout *psphotSelectBackgroundStdev (pmConfig *config, const pmFPAview *view, int index) {
+
+    pmFPAfile *file = pmFPAfileSelectSingle (config->files, psphotGetFilerule("PSPHOT.BACKMDL.STDEV"), index);
     if (!file) return NULL;
 
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotPetrosianAnalysis.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotPetrosianAnalysis.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotPetrosianAnalysis.c	(revision 37825)
@@ -12,11 +12,4 @@
     // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
     float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
-
-# if (0)
-    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
-    // XXX use this to set skynoise
-    pmReadout *backModel = psphotSelectBackground (config, view);
-    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
-# endif
 
     // S/N limit to perform full non-linear fits
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadout.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadout.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadout.c	(revision 37825)
@@ -325,4 +325,6 @@
     psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
     psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
+    // measure some parameters for galaxy science
+    psphotGalaxyParams (config, view, filerule);
     psphotRadialApertures(config, view, filerule, 0);
 
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadoutCleanup.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadoutCleanup.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotReadoutCleanup.c	(revision 37825)
@@ -57,5 +57,5 @@
     // use the psf-model to measure FWHM stats
     if (psf) {
-        if (!psphotPSFstats (readout, psf)) {
+      if (!psphotPSFstatsSources (readout, sources, psf)) {
             psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
             return false;
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetNFrames.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetNFrames.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetNFrames.c	(revision 37825)
@@ -45,8 +45,14 @@
 
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-    psAssert (detections, "missing detections?");
-    psAssert (detections->allSources, "missing sources?");
+    if (!detections) {
+        psLogMsg ("psphot", PS_LOG_INFO, "No detections image input %d", index);
+        return true;
+    }
 
     psArray *sources = detections->allSources;
+    if (!sources) {
+        psLogMsg ("psphot", PS_LOG_INFO, "No sources for input %d", index);
+        return true;
+    }
 
     int col0 = expnum->col0;
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetThreads.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetThreads.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSetThreads.c	(revision 37825)
@@ -87,4 +87,14 @@
     psFree(task);
 
+    task = psThreadTaskAlloc("PSPHOT_CHIP_PARAMS", 2);
+    task->function = &psphotChipParams_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_GALAXY_PARAMS", 11);
+    task->function = &psphotGalaxyParams_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     return true;
 }
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSourceStats.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSourceStats.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotSourceStats.c	(revision 37825)
@@ -185,18 +185,18 @@
             psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
 
-            psArrayAdd(job->args, 1, cells->data[j]); // sources
-
-            PS_ARRAY_ADD_SCALAR(job->args, INNER,   PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, MIN_SN,  PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32);
-
-            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
-            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
-
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
+            psArrayAdd(job->args, 1, cells->data[j]); // 0 : sources
+
+            PS_ARRAY_ADD_SCALAR(job->args, INNER,   PS_TYPE_F32); // 1
+            PS_ARRAY_ADD_SCALAR(job->args, MIN_SN,  PS_TYPE_F32); // 2
+            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32); // 3
+            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32); // 4
+            PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS, PS_TYPE_F32); // 5
+
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); // 6
+            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); // 7
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 8  : this is used as a return value for Nmoments
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 9  : this is used as a return value for Nfail
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // 10 : this is used as a return value for Nfaint
 
             if (!psThreadJobAddPending(job)) {
@@ -316,4 +316,5 @@
 
             // XXX: this must match the above
+	    psAbort ("this code is broken: must match psphotSourceStatsReadout above");
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, recipe);
@@ -472,11 +473,11 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[7];
+    scalar = job->args->data[8];
     scalar->data.S32 = Nmoments;
 
-    scalar = job->args->data[8];
+    scalar = job->args->data[9];
     scalar->data.S32 = Nfail;
 
-    scalar = job->args->data[9];
+    scalar = job->args->data[10];
     scalar->data.S32 = Nfaint;
 
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackParseCamera.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackParseCamera.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackParseCamera.c	(revision 37825)
@@ -262,4 +262,5 @@
     }
 
+#ifdef MAKE_CHISQ_IMAGE
     // generate an pmFPAimage for the chisqImage
     {
@@ -293,4 +294,7 @@
 	chisqVariance->save = saveChisq;
     }
+#else
+    (void) saveChisq;
+#endif   // MAKE_CHISQ_IMAGE
 
     psTrace("psphot", 1, "Done with psphotStackParseCamera...\n");
Index: branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackReadout.c
===================================================================
--- branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackReadout.c	(revision 37819)
+++ branches/eam_branches/ipp-pv3-20140717-merge/psphot/src/psphotStackReadout.c	(revision 37825)
@@ -341,4 +341,7 @@
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
     psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
+
+    // measure some parameters for galxy science
+    psphotGalaxyParams (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // create source children for the OUT filerule (for radial aperture photometry and output) 
@@ -354,4 +357,7 @@
 	return psphotReadoutCleanup (config, view, STACK_RAW);
     }
+    // These arrays are no longer needed. The inputs' sources arrays have references to the data
+    psFree (objects);
+    psFree (objectsOut);
 
     bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
@@ -410,6 +416,4 @@
     // measure aperture photometry corrections
     if (!psphotApResid (config, view, STACK_RAW)) {
-	psFree (objects);
-	psFree (objectsOut);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
 	return psphotReadoutCleanup (config, view, STACK_RAW);
@@ -432,5 +436,5 @@
 
     logMemStats("final");
-#if (1)
+#if (0)
     psphotSourceMemory(config, view, STACK_RAW);
     psphotSourceMemory(config, view, STACK_OUT);
@@ -448,7 +452,4 @@
     psphotStackRemoveChisqFromInputs(config, STACK_RAW);
 #endif
-
-    psFree (objects);
-    psFree (objectsOut);
 
     // create the exported-metadata and free local data