Changeset 36680 for branches/eam_branches/ps2-tc3-20130727/psphot

Timestamp:

Apr 21, 2014, 5:42:34 AM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ps2-tc3-20130727

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/Makefile.am (modified) (9 diffs)
• psphot/src/psmakecff.c (copied) (copied from trunk/psphot/src/psmakecff.c )
• psphot/src/psphot.h (modified) (9 diffs)
• psphot/src/psphotApResid.c (modified) (2 diffs)
• psphot/src/psphotArguments.c (modified) (1 diff)
• psphot/src/psphotBlendFit.c (modified) (3 diffs)
• psphot/src/psphotChooseAnalysisOptions.c (copied) (copied from trunk/psphot/src/psphotChooseAnalysisOptions.c )
• psphot/src/psphotChoosePSF.c (modified) (1 diff)
• psphot/src/psphotCleanup.c (modified) (1 diff)
• psphot/src/psphotDeblendSatstars.c (modified) (3 diffs)
• psphot/src/psphotDefineFiles.c (modified) (1 diff)
• psphot/src/psphotEllipticalContour.c (modified) (1 diff)
• psphot/src/psphotExtendedSourceAnalysis.c (modified) (8 diffs)
• psphot/src/psphotExtendedSourceFits.c (modified) (20 diffs)
• psphot/src/psphotFullForce.SourceStats.c (copied) (copied from trunk/psphot/src/psphotFullForce.SourceStats.c )
• psphot/src/psphotFullForce.c (copied) (copied from trunk/psphot/src/psphotFullForce.c )
• psphot/src/psphotFullForceArguments.c (copied) (copied from trunk/psphot/src/psphotFullForceArguments.c )
• psphot/src/psphotFullForceImageLoop.c (copied) (copied from trunk/psphot/src/psphotFullForceImageLoop.c )
• psphot/src/psphotFullForceReadout.c (copied) (copied from trunk/psphot/src/psphotFullForceReadout.c )
• psphot/src/psphotFullForceSummary.c (copied) (copied from trunk/psphot/src/psphotFullForceSummary.c )
• psphot/src/psphotFullForceSummaryReadout.c (copied) (copied from trunk/psphot/src/psphotFullForceSummaryReadout.c )
• psphot/src/psphotGalaxyShape.c (copied) (copied from trunk/psphot/src/psphotGalaxyShape.c )
• psphot/src/psphotImageLoop.c (modified) (1 diff)
• psphot/src/psphotKronFlux.c (copied) (copied from trunk/psphot/src/psphotKronFlux.c )
• psphot/src/psphotKronIterate.c (modified) (1 diff)
• psphot/src/psphotLensing.c (copied) (copied from trunk/psphot/src/psphotLensing.c )
• psphot/src/psphotLoadPSF.c (modified) (2 diffs)
• psphot/src/psphotMergeSources.c (modified) (9 diffs)
• psphot/src/psphotModelTestReadout.c (modified) (1 diff)
• psphot/src/psphotOldCode.c (copied) (copied from trunk/psphot/src/psphotOldCode.c )
• psphot/src/psphotPetroFlux.c (copied) (copied from trunk/psphot/src/psphotPetroFlux.c )
• psphot/src/psphotPetrosianRadialBins.c (modified) (3 diffs)
• psphot/src/psphotPetrosianStats.c (modified) (7 diffs)
• psphot/src/psphotRadialBins.c (modified) (2 diffs)
• psphot/src/psphotRadialProfile.c (modified) (7 diffs)
• psphot/src/psphotRadiusChecks.c (modified) (1 diff)
• psphot/src/psphotReadout.c (modified) (3 diffs)
• psphot/src/psphotRoughClass.c (modified) (1 diff)
• psphot/src/psphotSetThreads.c (modified) (1 diff)
• psphot/src/psphotSourceFits.c (modified) (7 diffs)
• psphot/src/psphotSourceSize.c (modified) (15 diffs)
• psphot/src/psphotSourceStats.c (modified) (3 diffs)
• psphot/src/psphotStackImageLoop.c (modified) (9 diffs, 1 prop)
• psphot/src/psphotStackMatchPSFsNext.c (modified) (2 diffs)
• psphot/src/psphotStackMatchPSFsUtils.c (modified) (1 diff)
• psphot/src/psphotStackReadout.c (modified) (16 diffs)
• psphot/test/tap_psphot_galaxygrid.pro (modified) (19 diffs)

branches/eam_branches/ps2-tc3-20130727/psphot/src

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary

branches/eam_branches/ps2-tc3-20130727/psphot/src/Makefile.am

@@ -25 +25 @@
 libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 
-bin_PROGRAMS = psphot psphotForced psphotMinimal psphotMakePSF psphotStack psphotModelTest
+bin_PROGRAMS = psphot psphotForced psphotFullForce psphotFullForceSummary psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff
 # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy 
 
@@ -36 +36 @@
 psphotForced_LDADD = libpsphot.la
 
+psphotFullForce_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotFullForce_LDADD = libpsphot.la
+
+psphotFullForceSummary_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psphotFullForceSummary_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psphotFullForceSummary_LDADD = libpsphot.la
+
 psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
 psphotMinimal_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
@@ -51 +59 @@
 psphotModelTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
 psphotModelTest_LDADD = libpsphot.la
+
+psmakecff_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
+psmakecff_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
+psmakecff_LDADD = libpsphot.la
+
 
 # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
@@ -81 +94 @@
         psphotMosaicChip.c         \
         psphotCleanup.c
+
+# generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf
+psphotFullForce_SOURCES = \
+        psphotFullForce.c             \
+        psphotFullForceArguments.c    \
+        psphotParseCamera.c        \
+        psphotImageLoop.c          \
+        psphotMosaicChip.c         \
+        psphotCleanup.c
+
+# combine full force results from several inputs
+psphotFullForceSummary_SOURCES = \
+        psphotFullForceSummary.c \
+        psphotFullForceSummaryReadout.c
 
 # forced photometry of specified positions given a specified psf
@@ -129 +156 @@
         psphotCleanup.c
 
+# a program that takes a cmf file and makes a cff file (input file for psphotFullForce
+psmakecff_SOURCES = \
+        psmakecff.c
+
 # psphotTest_SOURCES = \
 #         psphotTest.c
@@ -155 +186 @@
         psphotReadoutMinimal.c         \
         psphotForcedReadout.c          \
+        psphotFullForceReadout.c       \
+        psphotFullForce.SourceStats.c \
         psphotMakePSFReadout.c         \
         psphotModelBackground.c        \
@@ -175 +208 @@
         psphotMakeGrowthCurve.c        \
         psphotMagnitudes.c             \
+        psphotLensing.c        \
         psphotSetMaskBits.c            \
         psphotSkyReplace.c             \
@@ -181 +215 @@
         psphotSourceFits.c             \
         psphotRadiusChecks.c           \
+        psphotChooseAnalysisOptions.c  \
         psphotOutput.c                 \
         psphotFakeSources.c            \
@@ -198 +233 @@
         psphotRadialPlot.c             \
         psphotKronIterate.c            \
+        psphotKronFlux.c            \
+        psphotPetroFlux.c            \
+        psphotGalaxyShape.c            \
         psphotRadialProfileWings.c     \
         psphotDeblendSatstars.c        \

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphot.h

@@ -18 +18 @@
     PSPHOT_SINGLE,
     PSPHOT_FORCED,
+    PSPHOT_FULL_FORCE,
     PSPHOT_MAKE_PSF,
     PSPHOT_MODEL_TEST,
@@ -81 +82 @@
 bool            psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow);
 
+bool            psphotFullForceSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow);
+bool            psphotFullForceSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool setWindow);
+
 bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
@@ -122 +126 @@
 bool            psphotAddOrSubNoise_Threaded (psThreadJob *job);
 
+bool            psphotChooseAnalysisOptions (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotChooseAnalysisOptionsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotChooseAnalysisOptionsByObject(pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objects);
+
 bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
@@ -136 +144 @@
 bool            psphotMagnitudesReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf);
 bool            psphotMagnitudes_Threaded (psThreadJob *job);
+
+bool            psphotLensing (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool            psphotLensingReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotLensingPSFtrendsReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
 
 bool            psphotEfficiency (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -341 +353 @@
 bool psphotForcedReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
 
+pmConfig *psphotFullForceArguments(int argc, char **argv);
+bool psphotFullForceReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
+
 pmConfig *psphotMinimalArguments(int argc, char **argv);
 bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -349 +364 @@
 pmConfig *psphotModelTestArguments(int argc, char **argv);
 bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule);
+
+bool psphotFullForceSummaryReadout (pmConfig * config, const pmFPAview *view);
 
 int psphotFileruleCount(const pmConfig *config, const char *filerule);
@@ -414 +431 @@
 } psphotStackOptions;
 
+typedef struct {
+    float fRmajorMin;
+    float fRmajorMax;
+    float fRmajorDel;
+    float fRminorMin;
+    float fRminorMax;
+    float fRminorDel;
+} psphotGalaxyShapeOptions;
+
 /*** psphotStackMatchPSF prototypes ***/
 bool psphotStackMatchPSFs (pmConfig *config, const pmFPAview *view);
@@ -475 +501 @@
 bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
 bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted);
 
 bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
@@ -496 +522 @@
 bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index, int pass);
 bool psphotKronIterate_Threaded (psThreadJob *job);
+
+bool psphotKronFlux (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotKronFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources);
+bool psphotKronFlux_Threaded (psThreadJob *job);
+bool psphotKronFluxSource (pmSource *source, psImageMaskType maskVal);
+
+bool psphotPetroFlux (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotPetroFluxReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources);
+bool psphotPetroFlux_Threaded (psThreadJob *job);
+bool psphotPetroFluxSource (pmSource *source, psImageMaskType maskVal);
+
+bool psphotGalaxyShape (pmConfig *config, const pmFPAview *view, const char *filerule);
+bool psphotGalaxyShapeReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf);
+bool psphotGalaxyShape_Threaded (psThreadJob *job);
+bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize);
+bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults);
+psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc();
 
 bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule);

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotApResid.c

@@ -214 +214 @@
 
         // XXX make this user-configurable?
-        if (source->psfMagErr > 0.01) continue;
+        if (source->psfMagErr > 0.03) continue;
 
         // aperture residual for this source
@@ -267 +267 @@
 
     // set the max order (0 = constant) which the number of psf stars can support:
-    // rule of thumb: require 3 stars per 'cell' (order+1)^2
+    // we require only 3 stars for n = 0, increase stars / cell for higher order
     int MaxOrderForStars = 0;
-    if (Npsf >= 12) MaxOrderForStars = 1; // 4 cells
-    if (Npsf >= 27) MaxOrderForStars = 2; // 9 cells
-    if (Npsf >= 48) MaxOrderForStars = 3; // 16 cells
-    if (Npsf >  75) MaxOrderForStars = 4; // 25 cells
+    if (Npsf >=  16) MaxOrderForStars = 1; // 4 cells
+    if (Npsf >=  54) MaxOrderForStars = 2; // 9 cells
+    if (Npsf >= 128) MaxOrderForStars = 3; // 16 cells
+    if (Npsf >= 300) MaxOrderForStars = 4; // 25 cells
+    if (Npsf >  576) MaxOrderForStars = 5; // 36 cells
 
     pmTrend2DMode mode = PM_TREND_MAP;

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotArguments.c

@@ -217 +217 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf",      "-psflist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC",        "-src",      "-srclist");
+    // XXX allow this format? pmConfigFileSetsMD (config->arguments, &argc, argv, "FORCE",      "-force",    "-forcelist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "EXPNUM",     "-expnum",   "-expnumlist");
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotBlendFit.c

@@ -273 +273 @@
         pmSource *source = sources->data[i];
 
-        int TEST_ON = false;
-# if (0)
+# if (PS_TRACE_ON)
 # define TEST_X 653
 # define TEST_Y 466
-        if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
+# define TESTING 1
+        int TEST_ON = false;
+        if (TESTING && (fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
             fprintf (stderr, "test object\n");
             psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
@@ -342 +343 @@
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
             if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
+# if (PS_TRACE_ON)
                 if (TEST_ON) {
                     psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
                 }
+# endif
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
@@ -351 +354 @@
         } else {
             if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) {
+# if (PS_TRACE_ON)
                 if (TEST_ON) {
                     psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
                 }
+# endif
                 source->type = PM_SOURCE_TYPE_STAR;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChoosePSF.c

@@ -399 +399 @@
     char *modelName = pmModelClassGetName (psf->type);
     psLogMsg ("psphot.pspsf", PS_LOG_WARN, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf"));
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f (%d x %d model)\n", modelName, psf->ApResid, psf->dApResid, psf->trendNx, psf->trendNy);
 
     psFree (options);

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotCleanup.c

@@ -22 +22 @@
     pmVisualCleanup ();
     psLibFinalize();
+#if (PS_TRACE_ON)
+    // 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");
+#endif
     return;
 }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotDeblendSatstars.c

@@ -128 +128 @@
     // int display = psphotKapaChannel (1);
     // psphotVisualScaleImage (display, (psImage *) source->pixels->parent, NULL, "image", 1.0, 0);
+
+    int Nsatstar = 0;
 
     // examine sources in decreasing SN order
@@ -179 +181 @@
         source->mode2 |= PM_SOURCE_MODE2_SATSTAR_PROFILE; // and we have in fact subtracted the profile
 
+        Nsatstar ++;
+
         // XXX visualize, model, and subtract
         // if (!psphotVisualRadialProfileSatstar (source, maskVal)) {
@@ -192 +196 @@
     psFree (index);
 
-    psLogMsg ("psphot", PS_LOG_INFO, "deblend satstar: %f sec\n", psTimerMark ("psphot.deblend.sat"));
+    psLogMsg ("psphot", PS_LOG_INFO, "deblend %d satstars: %f sec\n", Nsatstar, psTimerMark ("psphot.deblend.sat"));
     return true;
 }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotDefineFiles.c

@@ -150 +150 @@
     }
 
+    if (psMetadataLookupPtr(NULL, config->arguments, "FORCE")) {
+        if (!pmFPAfileDefineFromArgs (&status, config, "PSPHOT.INPUT.CFF", "FORCE")) {
+            psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.INPUT.CFF");
+            return status;
+        }
+    }
+
     if (psMetadataLookupPtr(NULL, config->arguments, "SRCTEXT")) {
         // XXX cannot use pmFPAfileDefineFromArgs: this is explicitly a FITS-based I/O function

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEllipticalContour.c

@@ -64 +64 @@
         psFree (y);
         psFree (yErr);
+        source->mode2 |= PM_SOURCE_MODE2_ECONTOUR_FEW_PTS;
         return false;
     }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotExtendedSourceAnalysis.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-
-// measure the elliptical radial profile and use this to measure the petrosian parameters for the sources
+void psphotRadialProfileShowSkips ();
+
+// measure the petrosian parameters for the sources
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -15 +16 @@
     psAssert (recipe, "missing recipe?");
 
-    // perform full non-linear fits / extended source analysis?
-    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
+    bool doPetrosian = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    bool doAnnuli    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+
+    // measure petrosians?
+    if (!doPetrosian && !doAnnuli) {
         psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
         return true;
@@ -33 +37 @@
 }
 
+/*** 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;
+
+# define SKIP(VALUE) { VALUE++; continue; }
+
 // aperture-like measurements for extended sources
 bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
@@ -110 +125 @@
 
 // set this to 0 to run without threading
-# if (1)            
+# if (0)            
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
@@ -164 +179 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d annuli\n", Nannuli);
 
+# if (PS_TRACE_ON)
+    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);
+#endif
+
+    psphotRadialProfileShowSkips ();
+
     psphotVisualShowResidualImage (readout, false);
 
-    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
-    if (doPetrosian) {
-        psphotVisualShowPetrosians (sources);
-    }
+    psphotVisualShowPetrosians (sources);
 
     return true;
@@ -190 +213 @@
     float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
 
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "EXTENDED_SOURCE_SN_LIM");
-
-    // which extended source analyses should we perform?
     bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
-    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
-    bool doPetroStars   = psMetadataLookupBool (&status, recipe, "PETROSIAN_FOR_STARS");
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -207 +224 @@
         pmSource *source = sources->data[i];
 
-        // if we have checked the source validity on the basis of the object set, then 
-        // we either skip these tests below or we skip the source completely
-        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_SKIP) continue;
-        if (source->tmpFlags & PM_SOURCE_TMPF_PETRO_KEEP) goto keepSource;
-
-        // skip PSF-like and non-astronomical objects
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
-        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
-
-        // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
-
-        // optionally allow non-extended objects to get petrosians as well
-        if (!doPetroStars) {
-            if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-            if (source->type == PM_SOURCE_TYPE_STAR) continue;
-        }
-
-        // limit selection to some SN limit
-        // assert (source->peak); // how can a source not have a peak?
-        // limit selection to some SN limit
-        bool skipSource = false;
-        if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            skipSource = (source->moments->KronFlux < SN_LIM * source->moments->KronFluxErr);
-        } else {
-            skipSource = (sqrt(source->peak->detValue) < SN_LIM);
-        }
-        if (skipSource) continue;
-
-        // limit selection by analysis region (this automatically apply
-        if (source->peak->x < region->x0) continue;
-        if (source->peak->y < region->y0) continue;
-        if (source->peak->x > region->x1) continue;
-        if (source->peak->y > region->y1) continue;
-
-    keepSource:
+        Nall ++;
+
+        // 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);
 
         // replace object in image
@@ -259 +249 @@
         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
 
-        // if we request any of these measurements, we require the radial profile
-        if (doPetrosian || doAnnuli) {
-            if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
-                // all measurements below require the radial profile; skip them all
-                // re-subtract the object, leave local sky
-                psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-                continue;
-            }
-            Nannuli ++;
-            source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
-        }
+        // measure the radial profile
+        if (!psphotRadialProfile (source, recipe, skynoise, maskVal)) {
+          // re-subtract the object, leave local sky
+          psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+          SKIP (Nskip6);
+        }
+
+        Nannuli ++;
+        source->mode |= PM_SOURCE_MODE_RADIAL_FLUX;
 
         // Petrosian Mags
         if (doPetrosian) {
-            if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
-                psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-            } else {
-                psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
-                Npetro ++;
-                source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
-            }
+          if (!psphotPetrosian (source, recipe, skynoise, maskVal)) {
+            psTrace ("psphot", 5, "FAILED petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+          } else {
+            psTrace ("psphot", 5, "measured petrosian flux & radius for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+            Npetro ++;
+            source->mode |= PM_SOURCE_MODE_EXTENDED_STATS;
+          }
         }
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotExtendedSourceFits.c

@@ -111 +111 @@
     }
 
+    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolution
+    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
+        fitNsigmaConv = 5.0;
+    }
+
     bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
     if (!status) {
@@ -123 +128 @@
     }
 
+    // perform full extended source non-linear fits?
+    bool isInteractive = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_INTERACTIVE");
+    if (!status) isInteractive = false;
+
     // Define source fitting parameters for extended source fits
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->mode           = PM_SOURCE_FIT_EXT_AND_SKY;
     fitOptions->saveCovariance = true;  // XXX make this a user option?
     fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
@@ -131 +140 @@
     fitOptions->minTol         = fitMinTol;
     fitOptions->maxTol         = fitMaxTol;
+    fitOptions->nsigma         = fitNsigmaConv;
 
     fitOptions->gainFactorMode   = gainFactorMode;
     fitOptions->chisqConvergence = chisqConvergence;
+    fitOptions->isInteractive    = isInteractive;
+
+    // use poissonian errors or local-sky errors
+    fitOptions->poissonErrors = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS_POISSON");
+    if (!status) fitOptions->poissonErrors = true;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -278 +293 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, models);
+
             // Allocate a metadata iterator here because psMetadataIteratorAlloc/Free are not thread safe
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
@@ -296 +312 @@
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
 
-// set this to 0 to run without threading
-# if (1)            
-            if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree(AnalysisRegion);
-                psFree (fitOptions);
-                psFree (models);
-                psphotSersicModelClassCleanup();
-                return false;
-            } 
-# else
-            if (!psphotExtendedSourceFits_Threaded(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree(AnalysisRegion);
-                psFree (fitOptions);
-                psFree (models);
-                psphotSersicModelClassCleanup();
-                return false;
+// XXX TEST 
+            if (!isInteractive) {
+                if (!psThreadJobAddPending(job)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                    psFree(AnalysisRegion);
+                    psFree (fitOptions);
+                    psFree (models);
+                    psphotSersicModelClassCleanup();
+                    return false;
+                } 
+            } else {
+                // run without threading
+                if (!psphotExtendedSourceFits_Threaded(job)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                    psFree(AnalysisRegion);
+                    psFree (fitOptions);
+                    psFree (models);
+                    psphotSersicModelClassCleanup();
+                    return false;
+                }
+                psScalar *scalar = NULL;
+                scalar = job->args->data[9];
+                Next += scalar->data.S32;
+                scalar = job->args->data[10];
+                Nconvolve += scalar->data.S32;
+                scalar = job->args->data[11];
+                NconvolvePass += scalar->data.S32;
+                scalar = job->args->data[12];
+                Nplain += scalar->data.S32;
+                scalar = job->args->data[13];
+                NplainPass += scalar->data.S32;
+                scalar = job->args->data[14];
+                Nfaint += scalar->data.S32;
+                scalar = job->args->data[15];
+                Nfail += scalar->data.S32;
+                psFree(job->args->data[3]); // iterator allocated above
+                psFree(job);
             }
-            psScalar *scalar = NULL;
-            scalar = job->args->data[9];
-            Next += scalar->data.S32;
-            scalar = job->args->data[10];
-            Nconvolve += scalar->data.S32;
-            scalar = job->args->data[11];
-            NconvolvePass += scalar->data.S32;
-            scalar = job->args->data[12];
-            Nplain += scalar->data.S32;
-            scalar = job->args->data[13];
-            NplainPass += scalar->data.S32;
-            scalar = job->args->data[14];
-            Nfaint += scalar->data.S32;
-            scalar = job->args->data[15];
-            Nfail += scalar->data.S32;
-            psFree(job->args->data[3]); // iterator allocated above
-            psFree(job);
-# endif
         }
 
@@ -418 +435 @@
     // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
+    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
+    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;
+
     // choose the sources of interest
     for (int i = 0; i < sources->n; i++) {
@@ -423 +448 @@
         pmSource *source = sources->data[i];
 
-        // skip PSF-like and non-astronomical objects
-        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-        // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
-
-        // XXX this should use peak?
+        // rules for measuring petrosian parameters for specific objects are set in
+        // psphotChooseAnalysisOptions.c
+        if (!(source->tmpFlags & PM_SOURCE_TMPF_EXT_FIT)) continue;
+
+        // limit selection by analysis region (XXX move this into psphotChooseAnalysisOption?)
         if (source->peak->x < region->x0) continue;
         if (source->peak->y < region->y0) continue;
         if (source->peak->x > region->x1) continue;
         if (source->peak->y > region->y1) continue;
-
-
-        // XXX for a test, just do the obvious trail
-        // XXX if (source->peak->xf < 1100) continue;
-        // XXX if (source->peak->xf > 1400) continue;
-        // XXX if (source->peak->yf >  245) continue;
 
         // replace object in image
@@ -452 +467 @@
         psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
 
-        // XXX WATCH OUT HERE!!
-        // fitRadius = 30;
-
         // UPDATE : we have changed the moments calculation.  There is now an iteration within 
         // psphotKronMasked to determine moments appropriate for a larger object.  The values
@@ -464 +476 @@
             source->modelFits = psArrayAllocEmpty (models->list->n);
         }
+
+# ifdef TEST_OBJECT
+        bool testObject = false;
+        testObject |= ((fabs(source->peak->xf -  179) < 5) && (fabs(source->peak->yf - 1138) < 5));
+        if (testObject) {
+            fprintf (stderr, "test object @ %f, %f\n", source->peak->xf, source->peak->yf);
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 5);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 5);
+        }
+# endif 
 
         // loop here over the models chosen for each source (exclude by S/N)
@@ -479 +501 @@
           assert (status);
 
-          // limit selection to some SN limit
-          bool skipSource = false;
-          if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-              skipSource = (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr);
-          } else {
-              skipSource = (sqrt(source->peak->detValue) < FIT_SN_LIM);
-          }
-          if (skipSource) {
+          // limit selection to some SN limit for specific models (this value only applies if > EXTENDED_SOURCE_SN_LIM)
+          if (isfinite(FIT_SN_LIM)) {
+            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) {
               Nfaint ++;
               continue;
+            }
           }
+
+          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RUN;
 
           // check on the model type
@@ -507 +527 @@
                   psTrace ("psphot", 5, "failed to fit psf-conv model for object at %f, %f", source->moments->Mx, source->moments->My);
                   Nfail ++;
+                  source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
                   continue;
               }
@@ -512 +533 @@
                        source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
               Nconvolve ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+              if (!(modelFit->flags & badModel)) {
                   NconvolvePass ++;
                   source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -526 +547 @@
                       Nfail ++;
                       doneFits = true;
+                      source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_FAIL;
                       continue;
                   }
                   psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
                   Nplain ++;
-                  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+                  if (!(modelFit->flags & badModel)) {
                       NplainPass ++;
                       source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
@@ -544 +566 @@
                       float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
 
-                      if (false && (source->peak->xf > 1100) && 
-                          (source->peak->xf < 1400) && 
-                          (source->peak->yf < 245)) {
+                      if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
                           fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 
                                    source->pixels->col0, source->pixels->row0, 
@@ -556 +576 @@
                           fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
                           psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+                          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_RETRY;
                       }
                   }
               }
           }
-          // XXX really need to do this in a cleaner way:
-          if (!modelFit) continue;
+          psAssert (modelFit, "modelFit not set?");
 
           // test for fit quality / result
@@ -570 +590 @@
         }
 
+        // we are allowed to fit both stars and non-stars here -- if we have fitted
+        // something which we think is a star, we should use that model to subtract the
+        // object from the image.
+        if (source->type == PM_SOURCE_TYPE_STAR) {
+          // ensure the modelPSF is cached
+          pmSourceCacheModel (source, maskVal);
+          pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+          continue;
+        }
+
         // evaluate the relative quality of the models, choose one
+        // the PSF model might be the best fit : allow it to succeed
         float minChisq = NAN;
         int minModel = -1;
@@ -576 +607 @@
             pmModel *model = source->modelFits->data[i];
 
+            // skip the really bad fits
             if (!(model->flags & PM_MODEL_STATUS_FITTED)) continue;
-
-            if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue;
-            if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
-            if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue;
+            if (model->flags & badModel) continue;
+
+            // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
 
             if ((minModel < 0) || (model->chisq < minChisq)) {
@@ -598 +629 @@
           pmSourceCacheModel (source, maskVal);
 
+# if (PS_TRACE_ON)
+          pmModel *model = source->modelFits->data[0];
+          int flags = 0xffffffff;
+          if (model) {
+            flags = model->flags;
+          }
+          fprintf (stderr, "failed to fit extended source model to object %d @ %f, %f (%x)\n", source->id, source->moments->Mx, source->moments->My, flags);
+#endif
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
+          source->mode2 |= PM_SOURCE_MODE2_EXT_FITS_NONE;
           continue;
         }
@@ -615 +655 @@
         // cache the model flux
         if (source->modelEXT->isPCM) {
-            pmPCMCacheModel (source, maskVal, psfSize);
+            // fprintf (stderr, "subtract PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
+            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
         } else {
+            // fprintf (stderr, "subtract non-PCM extended source model for object %d @ %f, %f\n", source->id, source->moments->Mx, source->moments->My);
             pmSourceCacheModel (source, maskVal);
         }
+        source->modelEXT->flags |= PM_MODEL_BEST_FIT;
 
         // subtract the best fit from the object, leave local sky
@@ -625 +668 @@
         psTrace ("psphot", 4, "best ext model for %f, %f : %s chisq = %f\n", source->moments->Mx, source->moments->My, pmModelClassGetName (source->modelEXT->type), source->modelEXT->chisq);
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
+
+# ifdef TEST_OBJECT
+        if (testObject) {
+            psTraceSetLevel ("psModules.objects.pmPCM_MinimizeChisq", 0);
+            psTraceSetLevel ("psphot.psphotExtendedSourceFits_Threaded", 0);
+        }
+# endif
     }
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotImageLoop.c

@@ -132 +132 @@
                     }
                     break;
+                  case PSPHOT_FULL_FORCE:
+                    if (!psphotFullForceReadout (config, view, "PSPHOT.INPUT")) {
+                        psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                        psFree (view);
+                        return false;
+                    }
+                    break;
                   case PSPHOT_MAKE_PSF:
                     if (!psphotMakePSFReadout (config, view, "PSPHOT.INPUT")) {

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotKronIterate.c

@@ -166 +166 @@
     }
 
+# if (PS_TRACE_ON)
     fprintf (stderr, "--- starting KRON ---\n");
+#endif
 
     // We measure the Kron Radius on a smoothed copy of the readout image

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLoadPSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+// PSPHOT.PSF.LOAD vs input file -- see note at top
+bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) {
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // Generate the mask and weight images, including the user-defined analysis region of interest
+        if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
 
 // NOTE : pmPSF_IO.c functions must load the psf model onto the chip->analysis metadata because
@@ -58 +75 @@
     return true;
 }
-
-// PSPHOT.PSF.LOAD vs input file -- see note at top
-bool psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule) {
-
-    int num = psphotFileruleCount(config, filerule);
-
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
-
-        // Generate the mask and weight images, including the user-defined analysis region of interest
-        if (!psphotLoadPSFReadout (config, view, filerule, "PSPHOT.PSF.LOAD", i)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i);
-            return false;
-        }
-    }
-    return true;
-}

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotMergeSources.c

@@ -39 +39 @@
     psAssert (newSources, "missing sources?");
 
-    // XXX TEST:
-    if (detections->allSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->allSources);
-    }
-    if (detections->newSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->newSources);
-    }
-
     if (!detections->allSources) {
         detections->allSources = psArrayAllocEmpty(newSources->n);
@@ -65 +57 @@
 // Merge the externally supplied sources with the existing sources.  Mark them as having mode
 // PM_SOURCE_MODE_EXTERNAL.
+
+// XXX this function needs to be updated slightly for psphotFullForce:
+// * load the additional parameters to guide the new concepts
 
 // XXX This function needs to be updated to loop over set of input files.  At the moment, we
@@ -73 +68 @@
     bool status;
     pmDetections *extCMF = NULL;
+    pmDetections *extCFF = NULL;
     psArray *extSourcesTXT = NULL;
     int index = 0;
@@ -102 +98 @@
     {
         pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF");
-        if (!readoutCMF) goto loadTXT;
+        if (!readoutCMF) goto loadCFF;
 
         extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS");
@@ -125 +121 @@
     }
 
+loadCFF:
+    // load data from input CFF file:
+    {
+        pmReadout *readoutCFF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CFF");
+        if (!readoutCFF) goto loadTXT;
+
+        extCFF = psMetadataLookupPtr (NULL, readoutCFF->analysis, "PSPHOT.DETECTIONS");
+        if (extCFF) {
+            psF32 exptime = psMetadataLookupF32(NULL, readout->parent->concepts, "CELL.EXPOSURE");
+            for (int i = 0; i < extCFF->allSources->n; i++) {
+                pmSource *source = extCFF->allSources->data[i];
+
+                // setting this bit not only tracks the inputs, it makes pmSourceMoments
+                // keep the Mx,My values for the centroid.  
+                source->mode |= PM_SOURCE_MODE_EXTERNAL;
+
+                // source->peak->detValue,rawFlux,smoothFlux all set to input flux value which is scaled
+                // to 1 second exposure time. Scale to this image's exposure.
+                source->peak->rawFlux    *= exptime;
+                source->peak->smoothFlux *= exptime;
+                source->peak->detValue   *= exptime;
+                // source->peak->xf,yf, moments->Mx,My all set to input position
+
+                // drop the loaded source modelPSF
+                psFree (source->modelPSF);
+                source->modelPSF = NULL;
+                source->imageID = index;
+
+                psArrayAdd (detections->newSources, 100, source);
+            }
+        }
+    }
+
 loadTXT:
 
@@ -157 +186 @@
     psFree (detections);
 
-    if (!extCMF && !extSourcesTXT) {
+    if (!(extCMF || extCFF || extSourcesTXT)) {
         psLogMsg ("psphot", 3, "no external sources for this readout");
         return true;
@@ -163 +192 @@
 
     int nCMF = extCMF        ? extCMF->allSources->n        : 0;
-    int nTXT = extSourcesTXT ? extSourcesTXT->n : 0;
-
-    psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d text) merged to yield %ld total sources",
-              nCMF + nTXT, nCMF, nTXT, sources->n);
+    int nCFF = extCFF        ? extCFF->allSources->n        : 0;
+    int nTXT = extSourcesTXT ? extSourcesTXT->n             : 0;
+
+    psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d cff, %d text) merged to yield %ld total sources",
+              nCMF + nCFF + nTXT, nCMF, nCFF, nTXT, sources->n);
     return true;
 }
@@ -812 +842 @@
 // array containing the child sources.  XXX currently, this is only used by psphotStackReadout
 // (sources go on allSources so that psphotChoosePSF can be called repeatedly)
-psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) {
+psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc, bool sourcesSubtracted) {
 
     bool status;
@@ -919 +949 @@
             // child sources have not been subtracted in this image, but this flag may be raised if
             // they were subtracted in the parent's image
-            sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            // XXX NOTE : in the pre-20130914 version of psphotStack, we carried a copy of the pixels 
+            // generated before the subtraction took place (and then we smoothed to match the desired PSF).  
+            // in the new version, we copy the image after subtraction; we need to distinguish these cases
+            if (!sourcesSubtracted) {
+                sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            }
 
             // set the output detections:

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotModelTestReadout.c

@@ -125 +125 @@
 
             // get the source moments
-            status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal);
+          status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal);
             if (!status) psAbort("psSourceMoments error");
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotPetrosianRadialBins.c

@@ -36 +36 @@
     pmSourceRadialProfile *profile = source->extpars->petProfile;
 
-    float skyModelErrorSQ = PS_SQR(skynoise);
+    // float skyModelErrorSQ = PS_SQR(skynoise);
 
     int nMax = radiusMax;
@@ -133 +133 @@
         if (radius->data.F32[i] > Rmax) {
             // calculate the value for the nOut bin
-            float value, dvalue;
+            float value; // , dvalue;
             if (values->n > 0) {
                 psVectorStats (stats, values, NULL, NULL, 0);
                 value = stats->robustMedian;
-                dvalue = stats->robustStdev;
+                // dvalue = stats->robustStdev;
             } else {
                 value = NAN;
-                dvalue = NAN;
+                // dvalue = NAN;
             }
 
             binSB->data.F32[nOut] = value;
-            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
+            binSBstdev->data.F32[nOut] = skynoise / sqrt(values->n);
             binFill->data.F32[nOut] = values->n / binArea->data.F32[nOut];
 
@@ -185 +186 @@
         psFree(values);
         psFree(stats);
+        source->mode2 |= PM_SOURCE_MODE2_RADBIN_NAN_CENTER;
         return true;
     }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotPetrosianStats.c

@@ -6 +6 @@
 // generate the Petrosian radius and flux from the mean surface brightness (r_i)
 
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X);
 float InterpolateValues     (float X0, float Y0, float X1, float Y1, float X);
 float InterpolateValuesErrX (float X0, float Y0, float X1, float Y1, float X, float dX0, float dX1);
@@ -20 +21 @@
     if (!profile->binSB) {
         psLogMsg ("psphot", PS_LOG_DETAIL, "no petrosian profile, skipping source %f, %f", source->peak->xf, source->peak->yf);
+        source->mode2 |= PM_SOURCE_MODE2_PETRO_NO_PROFILE;
         return true;
     }
@@ -53 +55 @@
     float dFsum2 = 0.0;
 
-    float nSigma = 3.0;
+    float nSigma = 2.0;
     int lowestSignificantRadius = 0;
     float lowestSignificantRatio = 1.0;
@@ -117 +119 @@
                 petRadius    = InterpolateValues     (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
                 petRadiusErr = InterpolateValuesErrX (1.0, 0.0, petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, 0.0, petRatioErr->data.F32[nOut]);
-            } else {
-                petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
-                petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
+                source->mode2 |= PM_SOURCE_MODE2_PETRO_RATIO_ZEROBIN;
+            } else {
+              // petRadius    = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+              if (nOut > 1) {
+                petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-2], &refRadius->data.F32[nOut-2],   PETROSIAN_RATIO);
+              } else {
+                petRadius    = InterpolateValuesQuadratic (&petRatio->data.F32[nOut-3], &refRadius->data.F32[nOut-3], PETROSIAN_RATIO);
+              }
+# if (PS_TRACE_ON)
+              float petRadiusLinear = InterpolateValues     (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO);
+              if (fabs(petRadius - petRadiusLinear) > fabs(refRadius->data.F32[nOut] - refRadius->data.F32[nOut-1])) { 
+                fprintf (stderr, "big difference : %f vs %f\n", petRadius, petRadiusLinear); 
+              }
+#endif
+              petRadiusErr = InterpolateValuesErrX (petRatio->data.F32[nOut-1], refRadius->data.F32[nOut-1], petRatio->data.F32[nOut], refRadius->data.F32[nOut], PETROSIAN_RATIO, petRatioErr->data.F32[nOut-1], petRatioErr->data.F32[nOut]);
             }
             above = false;
@@ -148 +162 @@
             fprintf (stderr, "nan pet radius\n");
         }
+        source->mode2 |= PM_SOURCE_MODE2_PETRO_INSIG_RATIO;
     }
 
@@ -180 +195 @@
     bool found50 = false;
     bool found90 = false;
+
     // XXX use bisection to do this faster:
     for (int i = 0; !(found50 && found90) && i < refRadius->n; i++) {
@@ -235 +251 @@
 }
 
+// Lagrange's form of the interpolating polynomial...
+float InterpolateValuesQuadratic (float *Xin, float *Yin, float X) {
+
+  float dx01 = Xin[0] - Xin[1];
+  float dx02 = Xin[0] - Xin[2];
+  float dx12 = Xin[1] - Xin[2];
+
+  float dx0  = X - Xin[0];
+  float dx1  = X - Xin[1];
+  float dx2  = X - Xin[2];
+
+  float y0 = Yin[0]*dx1*dx2/(dx01*dx02);
+  float y1 = Yin[1]*dx0*dx2/(dx01*dx12); // need - sign
+  float y2 = Yin[2]*dx0*dx1/(dx02*dx12); 
+
+  float Y = y0 - y1 + y2;
+  return Y;
+}
+
 float InterpolateValues (float X0, float Y0, float X1, float Y1, float X) {
     float dydx = (Y1 - Y0) / (X1 - X0);

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadialBins.c

@@ -44 +44 @@
     psVector *radMin = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
     psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    if (!radMin || !radMin->n) {
-        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMin missing or empty)");
-        return false;
-    }
-    if (!radMax || !radMax->n) {
-        psError (PSPHOT_ERR_CONFIG, true, "error in definition of annular bins (radMax missing or empty)");
-        return false;
-    }
+    psAssert (radMin, "RADIAL.ANNULAR.BINS.LOWER is missing from recipe");
+    psAssert (radMin->n, "RADIAL.ANNULAR.BINS.LOWER is empty in recipe");
+    psAssert (radMax, "RADIAL.ANNULAR.BINS.UPPER is missing from recipe");
+    psAssert (radMax->n, "RADIAL.ANNULAR.BINS.UPPER is empty in recipe");
 
     psVector *binSB      = psVectorAllocEmpty(radMin->n, PS_TYPE_F32); // surface brightness of radial bin
@@ -161 +157 @@
         psFree(values);
         psFree(stats);
+        source->mode2 |= PM_SOURCE_MODE2_RADBIN_NAN_CENTER;
         return true;
     }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadialProfile.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+
+static int Nskip1 = 0;
+static int Nskip2 = 0;
+static int Nskip3 = 0;
+static int Nskip4 = 0;
+static int Nskip5 = 0;
+
+# define SKIP(VALUE) { VALUE++; return false; }
 
 bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal) {
@@ -23 +31 @@
     if (!psphotRadialProfilesByAngles (source, Nsec, Rmax)) {
         psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
-        return false;
+        SKIP (Nskip1);
     }
     // allocate: extpars->radFlux->radii,fluxes,theta
@@ -32 +40 @@
     if (!psphotRadiiFromProfiles (source, fluxMin, fluxMax)) {
         psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
-        return false;
+        SKIP (Nskip2);
     }
     // allocate : extpars->radFlux->isophotalRadii (use profile->radii,fluxes)
@@ -40 +48 @@
     if (!psphotEllipticalContour (source)) {
         // psLogMsg ("psphot", 3, "failed to measure elliptical contour");
-        return false;
+        SKIP (Nskip3);
     }
     // use extpars->radFlux->isophotalRadii,theta (result in extpars->axes)
@@ -48 +56 @@
     if (!psphotEllipticalProfile (source, RAW_RADIUS)) {
         psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
-        return false;
+        SKIP (Nskip4);
     }
     // allocate extpars->ellipticalFlux->radiusElliptical,fluxElliptical (use axes to scale raw pixels)
@@ -55 +63 @@
     if (!psphotRadialBins (recipe, source, Rmax, skynoise)) {
         psError (PS_ERR_UNKNOWN, false, "failed to generate radial bins");
-        return false;
+        SKIP (Nskip5);
     }
     // allocate extpars->radProfile->binSB, binSBstdv, binSum, binFill, radialBins, area (small lengths)
@@ -62 +70 @@
     return true;
 }
+
+void psphotRadialProfileShowSkips () {
+# if (PS_TRACE_ON)
+  fprintf (stderr, "radial profile skipped @ 1  : %d\n", Nskip1);
+  fprintf (stderr, "radial profile skipped @ 2  : %d\n", Nskip2);
+  fprintf (stderr, "radial profile skipped @ 3  : %d\n", Nskip3);
+  fprintf (stderr, "radial profile skipped @ 4  : %d\n", Nskip4);
+  fprintf (stderr, "radial profile skipped @ 5  : %d\n", Nskip5);
+#endif
+}

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadiusChecks.c

@@ -171 +171 @@
 # undef MIN_WINDOW
 
+// XXX EAM : 20130724 : for a test, double the window size parameters
+// # define MIN_WINDOW 10.0
+// # define SCALE1 7.0
 # define MIN_WINDOW 5.0
-# define SCALE1 5.0
+# define SCALE1 3.0
 # define PAD_WINDOW 3.0
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotReadout.c

@@ -316 +316 @@
     psphotSourceSize (config, view, filerule, false); // pass 2 (detections->allSources)
 
+    // XXX currently we are doing both the analysis of the size and the assessment of "fit ext"
+    // in source size.  this overloads the bit MODE_EXT_LIMIT to mean "fit ext" not just
+    // "bigger than a PSF"
+
+    // decide which source(s) are to be fitted with the extended source analysis code.
+    psphotChooseAnalysisOptions (config, view, filerule);
+
     psphotExtendedSourceAnalysis (config, view, filerule); // pass 1 (detections->allSources)
     psphotExtendedSourceFits (config, view, filerule); // pass 1 (detections->allSources)
@@ -333 +340 @@
     // calculate source magnitudes
     if (!psphotMagnitudes(config, view, filerule)) { // pass 1 (detections->allSources)
-      psErrorStackPrint(stderr, "Unable to do magnitudes.");
+        psErrorStackPrint(stderr, "Unable to do magnitudes.");
         psErrorClear();
     }
+
+    // calculate lensing parameters
+    if (!psphotLensing(config, view, filerule)) {
+        psErrorStackPrint(stderr, "Unable to do lensing parameters.");
+        psErrorClear();
+    }
+
     if (!psphotEfficiency(config, view, filerule)) { // pass 1
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
@@ -346 +360 @@
     // replace background in residual image
     if (!psphotSkyReplace (config, view, filerule)) { // pass 1
-      psErrorStackPrint(stderr, "Unable to replace sky");
-      psErrorClear();
+        psErrorStackPrint(stderr, "Unable to replace sky");
+        psErrorClear();
     }
 
     // drop the references to the image pixels held by each source
     if (!psphotSourceFreePixels (config, view, filerule)) { // pass 1
-      psErrorStackPrint(stderr, "Unable to free source pixels");
-      psErrorClear();
+        psErrorStackPrint(stderr, "Unable to free source pixels");
+        psErrorClear();
     }
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRoughClass.c

@@ -87 +87 @@
     if (NX > NY) {
         NXuse = ScaleForClump;
-        NYuse = (int) (ScaleForClump * (NY / NX) + 0.5);
+        NYuse = (int) (ScaleForClump * (NX / NY) + 0.5);
     } else {
         NYuse = ScaleForClump;

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSetThreads.c

@@ -47 +47 @@
     psFree(task);
 
+    task = psThreadTaskAlloc("PSPHOT_KRON_FLUX", 4);
+    task->function = &psphotKronFlux_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_PETRO_FLUX", 4);
+    task->function = &psphotPetroFlux_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_GALAXY_SHAPES", 7);
+    task->function = &psphotGalaxyShape_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     task = psThreadTaskAlloc("PSPHOT_BLEND_FIT", 10);
     task->function = &psphotBlendFit_Threaded;

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceFits.c

@@ -19 +19 @@
 static int NfitIterPCM = 0;
 static int NfitPixPCM = 0;
+
+bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
+bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize);
 
 bool psphotFitInit (int nThreads) {
@@ -559 +562 @@
 
 # define TIMING 0
+# define EXTRA_VERBOSE 0
+
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize);
+bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
+bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize);
 
 pmModel *psphotFitPCM (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
@@ -573 +582 @@
     maskVal |= markVal;
 
-    // allocate the model
+    // allocate the model (this can only fail on a config error)
     pmModel *model = pmModelAlloc(modelType);
-    if (!model) {
-        return NULL;
-    }
-
-    float t1, t2, t4, t5;
+    psAssert (model, "invalid extended model name");
+
+    float t1, t2, t3, t4, t5;
+    t1 = t2 = t3 = t4 = t5 = 0.0;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
+
+    // if we are ever (in a given psphot implementation) going to fit a parameter, we must set the options here to include
+    // that parameter (otherwise pmPCMupdate will fail to allocate the dmodelFlux image)
+    // thus, if the sersic analysis below uses an index fit, need to use this EXT_AND_SKY mode for init
+
+    options.mode = PM_SOURCE_FIT_EXT;
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_DEV")) {
+        options.mode = PM_SOURCE_FIT_SHAPE;
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_EXP")) {
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -594 +618 @@
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
         // use the source moments, etc to guess basic model parameters
-        if (!psphotSersicModelClassGuessPCM (pcm, source)) {
+        if (!psphotSersicModelGuessPCM (pcm, source, maskVal, psfSize)) {
             psFree (pcm);
-            model->flags |= PM_MODEL_STATUS_BADARGS;
+            model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
+            return model;
+        }
+        if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
+
+        // psphotFitSersicShapeAndIndex (pcm, readout, source, fitOptions, maskVal, markVal, psfSize);
+        options.mode = PM_SOURCE_FIT_NO_INDEX;
+        if (!psphotFitSersicShapeAndIndexGridAuto (pcm, readout, source, &options, maskVal, markVal, psfSize)) {
+            psFree (pcm);
+            model->flags |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
+            psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape");
+            psErrorClear (); // clear the polynomial error
             return model;
         }
@@ -603 +638 @@
         if (!pmSourceModelGuessPCM (pcm, source, maskVal, markVal)) {
             psFree (pcm);
-            model->flags |= PM_MODEL_STATUS_BADARGS;
+            model->flags |= PM_MODEL_PCM_FAIL_GUESS;
             return model;
         }
     }
 
-    if (TIMING) { t2 = psTimerMark ("psphotFitPCM"); }
-
-    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-        options.mode = PM_SOURCE_FIT_NO_INDEX;
-    } else {
-        options.mode = PM_SOURCE_FIT_EXT;
-    }
-    // update the pcm elements if we have changed the circumstance (options.mode or source->pixels)
-    pmPCMupdate(pcm, source, &options, model);
-    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
+    if (TIMING) { t3 = psTimerMark ("psphotFitPCM"); }
 
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
@@ -623 +649 @@
     NfitIterPCM += pcm->modelConv->nIter;
     NfitPixPCM += pcm->modelConv->nDOF;
+    if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
+
+    // XXX we might make this more efficient by setting NITER to be fairly small.  if we hit the iteration
+    // limit, then we could do a small grid search on the size and try again from the best fit 
+
+    if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
+    // if (pcm->modelConv->nIter == fitOptions->nIter) {
+    //  psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
+    // }
     if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
 
     if (TIMING) {
         int nPixBig = source->pixels->numCols * source->pixels->numRows;
-        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t4, t5);
-    }
-
-    // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+        fprintf (stderr, "psphotFitPCM : nIter: %2d, radius: %6.1f, npix: %5d of %5d, t1: %6.4f, t2: %6.4f, t3: %6.4f, t4: %6.4f, t5: %6.4f\n", model->nIter, model->fitRadius, model->nPix, nPixBig, t1, t2, t3, t4, t5);
+    }
+    if (EXTRA_VERBOSE && !TIMING) {
+        int nPixBig = source->pixels->numCols * source->pixels->numRows;
+        float *PAR = model->params->data.F32;
+        fprintf (stderr, "source %d : %f - %f %f - %f %f %f - %f | nIter: %2d, radius: %6.1f, npix: %5d of %5d, chisq %f\n", source->id, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1], model->nIter, model->fitRadius, model->nPix, nPixBig, model->chisqNorm);
+    }
+
     psFree (pcm);
 
@@ -748 +787 @@
     return true;
 }
+
+// 0.5 / n for (1.0, 1.25, 1.66, 2.0, 3.33, 4.0)
+// float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
+
+// 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
+float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 0.125, 0.10, 0.0833};
+float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
+# define N_INDEX_GUESS_INV 8
+
+// we are going to guess in fractions about the R1-based guess
+float reffGuess[] = {0.8, 0.9, 1.0, 1.12, 1.25};
+# define N_REFF_GUESS 5
+
+// A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotSersicModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
+
+    // we get a reasonable guess from:
+    // * Reff = Kron R1 / Q(index) -- Q comes from Graham & Driver 
+    // * Rmajor / Rminor & Theta from moments
+    // * Io from total Kron flux
+
+    // the guesses are used to fill in PAR:
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // convert the moments to Major,Minor,Theta
+    psEllipseMoments moments;
+
+    if (!isfinite(source->moments->Mrf)) return false;
+    if (!isfinite(source->moments->Mxx)) return false;
+    if (!isfinite(source->moments->Mxy)) return false;
+    if (!isfinite(source->moments->Myy)) return false;
+
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+    
+    // limit axis ratio < 20.0
+    psEllipseAxes momentAxes = psEllipseMomentsToAxes (moments, 20.0);
+
+    // set the model position
+    if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) {
+      return false;
+    }
+
+    // sky is zero (no longer fitted, but not yet deprecated)
+    PAR[PM_PAR_SKY]  = 0.0;
+
+    // for the index loop, use Io = 1.0, use fitted values to determine Io
+    PAR[PM_PAR_I0] = 1.0;
+
+    float xMin = NAN;
+    float iMin = NAN;
+    float sMin = NAN;
+    float rMin = NAN;
+
+    // loop over index and Reff, keeping the ARatio and Theta constant?
+    // loop over index guesses and find the best fit
+    for (int j = 0; j < N_REFF_GUESS; j++) {
+        for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+            PAR[PM_PAR_7] = indexGuessInv[i];
+
+            psEllipseAxes guessAxes;
+            guessAxes.major = reffGuess[j] * source->moments->Mrf / indexGuessR1q[i];
+            guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
+            guessAxes.theta = momentAxes.theta;
+
+            if (!isfinite(guessAxes.major)) return false;
+            if (!isfinite(guessAxes.minor)) return false;
+            if (!isfinite(guessAxes.theta)) return false;
+
+            // convert the major,minor,theta to shape parameters for an Reff-like model
+            pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
+            // generated the modelFlux
+            // XXX note that this does not add sky to model
+            pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
+        
+            float YY = 0.0;
+            float YM = 0.0;
+            float MM = 0.0;
+            bool usePoisson = false;
+
+            for (int iy = 0; iy < source->pixels->numRows; iy++) {
+                for (int ix = 0; ix < source->pixels->numCols; ix++) {
+                    // skip masked points
+                    if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                        continue;
+                    }
+                    // skip zero-variance points
+                    if (source->variance->data.F32[iy][ix] == 0) {
+                        continue;
+                    }
+                    // skip nan value points
+                    if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                        continue;
+                    }
+
+                    float fy = source->pixels->data.F32[iy][ix];
+                    float fm = source->modelFlux->data.F32[iy][ix];
+                    float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+                    YY += PS_SQR(fy) * wt;
+                    YM += fm * fy * wt;
+                    MM += PS_SQR(fm) * wt;
+                }
+            }
+
+            float Io = YM / MM;
+            float Chisq = YY - 2 * Io * YM + Io * Io * MM;
+            if (isnan(xMin) || (Chisq < xMin)) {
+                xMin = Chisq;
+                iMin = Io;
+                sMin = indexGuessInv[i];
+                rMin = reffGuess[j] / indexGuessR1q[i];
+            }
+            if (EXTRA_VERBOSE) {
+                fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
+                fprintf (stderr, "\n");
+            }
+        }
+    }
+
+    {
+        psEllipseAxes guessAxes;
+        guessAxes.major = rMin * source->moments->Mrf;
+        guessAxes.minor = guessAxes.major * (momentAxes.minor / momentAxes.major);
+        guessAxes.theta = momentAxes.theta;
+
+        if (!isfinite(guessAxes.major)) return false;
+        if (!isfinite(guessAxes.minor)) return false;
+        if (!isfinite(guessAxes.theta)) return false;
+
+        // convert the major,minor,theta to shape parameters for an Reff-like model
+        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+    }
+
+    PAR[PM_PAR_I0] = iMin;
+    PAR[PM_PAR_7] = sMin;
+
+    return true;
+}
+
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndex (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    for (int i = 0; i < 3; i++) {
+      // fit EXT (not PSF) model (set/unset the pixel mask)
+      options.mode = PM_SOURCE_FIT_SHAPE;
+      options.nIter = 2;
+
+      // update the pcm elements if we have changed the circumstance (here, options.mode)
+      pmPCMupdate(pcm, source, &options, model);
+      
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) {
+        float *PAR = model->params->data.F32;
+        fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+      }
+      
+      // fit EXT (not PSF) model (set/unset the pixel mask)
+      options.mode = PM_SOURCE_FIT_INDEX;
+      // options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+      options.nIter = 30;
+      
+      // update the pcm elements if we have changed the circumstance (here, options.mode)
+      pmPCMupdate(pcm, source, &options, model);
+      
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) {
+        float *PAR = model->params->data.F32;
+        fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+      }
+    }
+
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    return true;
+}
+
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGridAuto (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 7;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    // we have been provided a guess at the index (P[7]) from the list of indexGuessInv
+
+    // find the matching indexGuessInv
+    int nStart = -1;
+    for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+        if (fabs(PAR[PM_PAR_7] - indexGuessInv[i]) < 0.01) {
+            nStart = i;
+            break;
+        }
+    }
+    if (nStart == -1) {
+        fprintf (stderr, "WARNING: could not find start guess %f\n", PAR[PM_PAR_7]);
+        return false;
+    }
+
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    float Sm = NAN, Sp = NAN, So = NAN;
+    if (nStart == 0) {
+        Sm = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+        Sp = indexGuessInv[nStart + 1];
+    } else if (nStart == N_INDEX_GUESS_INV - 1) {
+        Sp = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        Sm = indexGuessInv[nStart - 1];
+    } else {
+        Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        So = indexGuessInv[nStart];
+        Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+    }
+    
+    PAR[PM_PAR_7] = Sm;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    PAR[PM_PAR_7] = So;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    PAR[PM_PAR_7] = Sp;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+        psFree (chi2);
+        psFree (Sidx);
+        psFree (poly);
+        return false;
+    }
+
+    // where is the minimum of this polynomial fit?
+    float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
+
+    // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further
+    float Smx = Sm - 0.25*(So - Sm);
+    float Spx = Sp + 0.25*(Sp - So);
+    Smin = PS_MAX(PS_MIN(Smin, Smx), Spx);
+    PAR[PM_PAR_7] = Smin;
+
+    // XXX I could set the error on PAR_7 here if I knew how to roughly convert these chisq values to true chisq values
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+ 
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGridAutoScaled (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 5;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    float parStart[8];
+    for (int i = 0; i < 8; i++) parStart[i] = PAR[i];
+
+    // we start with a guess at the index (P[7])
+
+    // get chisq for P[7], P[7]*1.1, P[7]*1.25 (or *0.75 depending on the result of 1.1)
+
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    PAR[PM_PAR_7] = parStart[7];
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    float fI = 1.1;
+    PAR[PM_PAR_7] = parStart[7]*fI;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    if (chi2->data.F32[1] < chi2->data.F32[0]) {
+      fI = 1.3;
+    } else {
+      fI = 1.0 / 1.3;
+    }
+    
+    PAR[PM_PAR_7] = parStart[7]*fI;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+    psVectorAppend (chi2, model->chisqNorm);
+
+    // can we fit the 3 pts with a parabola?
+    int nTry = 0;
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    while (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+      psErrorClear (); // clear the polynomial error
+      if (nTry > 4) {
+        psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+        psFree (chi2);
+        psFree (Sidx);
+        psFree (poly);
+        return false;
+      }
+      fI = (fI < 1.0) ? fI / 1.3 : fI * 1.3;
+      PAR[PM_PAR_7] = parStart[7]*fI;
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+      psVectorAppend (Sidx, 100*PAR[PM_PAR_7]);
+      psVectorAppend (chi2, model->chisqNorm);
+      nTry ++;
+    }
+
+    // where is the minimum of this polynomial fit?
+    float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
+
+    // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
+    Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
+    PAR[PM_PAR_7] = Smin;
+
+    // pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    // if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    
+    //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i];
+    //// 
+    //// for (float fI = 0.0; fI < 0.15; fI += 0.01) {
+    ////   PAR[PM_PAR_7] = parStart[7] - fI;
+    //// 
+    ////   // fit EXT (not PSF) model (set/unset the pixel mask)
+    ////   
+    ////   pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    ////   if (TIMING) {
+    ////        float *PAR = model->params->data.F32;
+    ////        fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]);
+    ////   }
+    //// }
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+ 
+// we have a set of guess parameters, do a small number of iterations fitting only SHAPE then only INDEX 
+bool psphotFitSersicShapeAndIndexGrid (pmPCMdata *pcm, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, int psfSize) {
+
+    pmModel *model = pcm->modelConv;
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    options.mode = PM_SOURCE_FIT_SHAPE;
+    options.nIter = 10;
+    
+    // update the pcm elements if we have changed the circumstance (here, options.mode)
+    pmPCMupdate(pcm, source, &options, model);
+    
+    psVector *chi2 = psVectorAllocEmpty (16, PS_TYPE_F32);
+    psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
+
+    float par7[] = {0.100, 0.125, 0.150, 0.175, 0.200, 0.225, 0.250};
+    for (int i = 0; i < 7; i++) {
+      PAR[PM_PAR_7] = par7[i];
+      pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+      if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+      psVectorAppend (Sidx, PAR[PM_PAR_7]);
+      psVectorAppend (chi2, model->chisqNorm);
+    }
+
+    psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 2);
+    if (!psVectorFitPolynomial1D (poly, NULL, 0, chi2, NULL, Sidx)) {
+      psError(PS_ERR_UNKNOWN, true, "Failed to find a good chisq parabola");
+      psFree (chi2);
+      psFree (Sidx);
+      psFree (poly);
+      return false;
+    }
+
+    // where is the minimum of this polynomial fit?
+    fprintf (stderr, "fit1d: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+    float Smin = -0.5 * poly->coeff[1] / poly->coeff[2];
+
+    // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
+    Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
+    PAR[PM_PAR_7] = Smin;
+    pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
+    
+    //// for (int i = 0; i < 8; i++) PAR[i] = parStart[i];
+    //// 
+    //// for (float fI = 0.0; fI < 0.15; fI += 0.01) {
+    ////   PAR[PM_PAR_7] = parStart[7] - fI;
+    //// 
+    ////   // fit EXT (not PSF) model (set/unset the pixel mask)
+    ////   
+    ////   pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
+    ////   if (TIMING) {
+    ////        float *PAR = model->params->data.F32;
+    ////        fprintf (stderr, "%d %f : %f - %f %f %f - %f\n", model->nIter, model->chisqNorm, PAR[7], PAR[4], PAR[5], PAR[6], PAR[1]);
+    ////   }
+    //// }
+
+    // return to the original fitting mode (fitOptions)
+    pmPCMupdate(pcm, source, fitOptions, model);
+
+    psFree (chi2);
+    psFree (Sidx);
+    psFree (poly);
+
+    return true;
+}
+
+// # define N_REFF_CHECK 11
+// float drefCheck[] = {-0.02, -0.04, -0.06, 0.0, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25};
+
+// we have an initial fit, check to see if the current size is besst
+bool psphotPCMfitCheckSize (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, float psfSize) {
+
+    // PAR is already at my current best guess
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // store best guess as a shape
+    psEllipseAxes centerAxes;
+    pmModelParamsToAxes (&centerAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
+
+    float xMin = NAN;
+    float iMin = NAN;
+    float rMin = NAN;
+
+    // loop over Reff, keeping the ARatio and Theta constant
+    for (int j = -4; j <= 4; j++) {
+
+        float dref = j * 0.01;
+
+        psEllipseAxes guessAxes;
+        guessAxes.major = centerAxes.major + dref;
+        guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
+        guessAxes.theta = centerAxes.theta;
+
+        if (!isfinite(guessAxes.major)) return false;
+        if (!isfinite(guessAxes.minor)) return false;
+        if (!isfinite(guessAxes.theta)) return false;
+
+        // convert the major,minor,theta to shape parameters for an Reff-like model
+        pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
+        // generated the modelFlux
+            // XXX note that this does not add sky to model
+        pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
+        
+        float YY = 0.0;
+        float YM = 0.0;
+        float MM = 0.0;
+        bool usePoisson = false;
+
+        for (int iy = 0; iy < source->pixels->numRows; iy++) {
+            for (int ix = 0; ix < source->pixels->numCols; ix++) {
+                // skip masked points
+                if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                    continue;
+                }
+                // skip zero-variance points
+                if (source->variance->data.F32[iy][ix] == 0) {
+                    continue;
+                }
+                // skip nan value points
+                if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                    continue;
+                }
+
+                float fy = source->pixels->data.F32[iy][ix];
+                float fm = source->modelFlux->data.F32[iy][ix];
+                float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+                YY += PS_SQR(fy) * wt;
+                YM += fm * fy * wt;
+                MM += PS_SQR(fm) * wt;
+            }
+        }
+
+        float Io = YM / MM;
+        float Chisq = YY - 2 * Io * YM + Io * Io * MM;
+        if (isnan(xMin) || (Chisq < xMin)) {
+            xMin = Chisq;
+            iMin = Io;
+            rMin = dref;
+        }
+        // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
+    }
+
+    psEllipseAxes guessAxes;
+    guessAxes.major = centerAxes.major + rMin;
+    guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
+    guessAxes.theta = centerAxes.theta;
+
+    if (!isfinite(guessAxes.major)) return false;
+    if (!isfinite(guessAxes.minor)) return false;
+    if (!isfinite(guessAxes.theta)) return false;
+
+    // convert the major,minor,theta to shape parameters for an Reff-like model
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+    PAR[PM_PAR_I0] = iMin;
+
+    return true;
+}
+
+// we have an initial fit, check to see if the current size is besst
+bool psphotPCMfitRetry (pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal, float psfSize) {
+
+    // PAR is already at my current best guess
+    psF32 *PAR = pcm->modelConv->params->data.F32;
+
+    // store best guess as a shape
+    psEllipseAxes centerAxes;
+    pmModelParamsToAxes (&centerAxes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true);
+
+    // retry with axes smaller by 1 pixel
+    psEllipseAxes guessAxes;
+    guessAxes.major = centerAxes.major - 0.08;
+    guessAxes.minor = guessAxes.major * centerAxes.minor / centerAxes.major;
+    guessAxes.theta = centerAxes.theta;
+
+    if (!isfinite(guessAxes.major)) return false;
+    if (!isfinite(guessAxes.minor)) return false;
+    if (!isfinite(guessAxes.theta)) return false;
+
+    // convert the major,minor,theta to shape parameters for an Reff-like model
+    pmModelAxesToParams (&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], guessAxes, true);
+
+    // generated the modelFlux
+            // XXX note that this does not add sky to model
+    pmPCMMakeModel (source, pcm->modelConv, pcm->nsigma, maskVal, psfSize);
+        
+    float YY = 0.0;
+    float YM = 0.0;
+    float MM = 0.0;
+    bool usePoisson = false;
+
+    for (int iy = 0; iy < source->pixels->numRows; iy++) {
+        for (int ix = 0; ix < source->pixels->numCols; ix++) {
+            // skip masked points
+            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+                continue;
+            }
+            // skip zero-variance points
+            if (source->variance->data.F32[iy][ix] == 0) {
+                continue;
+            }
+            // skip nan value points
+            if (!isfinite(source->pixels->data.F32[iy][ix])) {
+                continue;
+            }
+
+            float fy = source->pixels->data.F32[iy][ix];
+            float fm = source->modelFlux->data.F32[iy][ix];
+            float wt = (usePoisson) ? 1.0 / source->variance->data.F32[iy][ix] : 1.0;
+
+            YY += PS_SQR(fy) * wt;
+            YM += fm * fy * wt;
+            MM += PS_SQR(fm) * wt;
+        }
+    }
+
+    float Io = YM / MM;
+    PAR[PM_PAR_I0] = Io;
+
+    pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
+
+    return true;
+}
+
+

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceSize.c

@@ -1 +1 @@
 # include "psphotInternal.h"
-# include <gsl/gsl_sf_gamma.h>
-
+
+// this structure is only used internally to simplify the function parameters
 typedef struct {
     psImageMaskType maskVal;
@@ -13 +13 @@
     bool altDiffExt;
     float altDiffExtThresh;
-    bool extFitAll;
-    bool extFitAllReadout;
-    float extFitAllThresh;
     float soft;
     int grow;
@@ -33 +30 @@
 bool psphotSourceSelectCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal, int maxWindowCR);
-bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
 int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
 float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction);
@@ -41 +37 @@
 
 // we need to call this function after sources have been fitted to the PSF model and
-// subtracted.  To determine the CR-nature, this function examines the 9 pixels in the 3x3
-// square containing the peak and compares the observed flux to the model.  To determine
-// the EXT-nature, this function measures the amount of positive or negative total
-// deviation from the psf model at the r = FWHM/2 position
+// subtracted.  
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -121 +114 @@
     assert (status);
 
-    // XXX recipe name is not great
+    // XXX recipe name is not great (NOTE : not used!)
     options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS");
     assert (status);
 
-    // Optional extended source measurement algorithm to improve diff image trails
+    // Optional algorithm to define if a source is extended (used by DIFF analysis)
     options.altDiffExt = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE");
     assert (status);
+
     // Threshold for this alternate method
     options.altDiffExtThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.DIFF.ALTERNATE.THRESH");
     assert (status);
-    // Option to enable fitting of all objects with extended model.
-    options.extFitAll = psMetadataLookupBool(&status, recipe, "PSPHOT.EXT.FIT.ALL.SOURCES");
-    assert (status);
-    // Fitting everything is fine, but if the source density is high, we probably shouldn't.
-    options.extFitAllThresh = psMetadataLookupF32(&status, recipe, "PSPHOT.EXT.FIT.ALL.THRESH");
-    assert (status);
-    
+
     // location of a single test source
     options.xtest = psMetadataLookupS32 (&status, recipe, "PSPHOT.CRMASK.XTEST");
@@ -422 +410 @@
 
     psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %4s %4s %4s %4s %4s", "Npsf", "Next", "Nsat", "Ncr", "Nskip");
-    // Determine if this readout is above the threshold to ext fit all sources
-    options->extFitAllReadout = false;
-    if (options->extFitAll) {
-      float maskFrac = psMetadataLookupF32(&status,readout->analysis,"READOUT.MASK.FRAC");
-      if (status) {
-        maskFrac = 0.0;
-      }
-      if (sources->n * (1.0 - maskFrac) > options->extFitAllThresh) {
-        options->extFitAllReadout = false;
-      }
-      else {
-        options->extFitAllReadout = true;
-      }
-    }
-    
+
     if (!psphotSourceClassRegion (NULL, &psfClump, sources, recipe, psf, options)) {
         psLogMsg ("psphot", 4, "Failed to determine source classification for full image\n");
@@ -444 +418 @@
     return true;
     
+    // NOTE : this section is deactivated (EAM : I think we were getting poor boundary effects?)
     int nRegions = psMetadataLookupS32 (&status, readout->analysis, "PSF.CLUMP.NREGIONS");
     for (int i = 0; i < nRegions; i ++) {
@@ -582 +557 @@
         // * SAT stars should not be faint, but defects may?
 
-        // If the recipe requests we do extended source fits to everything, set
-        // the EXT_LIMIT flag
-        if (options->extFitAllReadout) {
-          psTrace("psphotSourceClassRegion.EXTALT",10,"In extFitAll: %d %d\n",options->extFitAll,options->extFitAllReadout);
-          source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
-        }
         // Defects may not always match CRs from peak curvature analysis
         // Defects may also be marked as SATSTAR -- XXX deactivate this flag?
@@ -595 +564 @@
 
         // saturated star (too many saturated pixels or peak above saturation limit).  These
-        // may also be saturated galaxies, or just large saturated regions.
+        // may also be saturated galaxies, or just large saturated regions.  They are never
+        // marked as 'extended'
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
             psTrace("psphotSourceClassRegion.SAT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g SAT\n",
@@ -604 +574 @@
         }
 
-        // any sources missing a large fraction should just be treated as PSFs
+        // any sources missing a large fraction should just be treated as PSFs. They are never
+        // marked as 'extended'
         if ((source->pixWeightNotBad < 0.9) || (source->pixWeightNotPoor < 0.9)) {
             psTrace("psphotSourceClassRegion.PSF",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
@@ -634 +605 @@
             psTrace("psphotSourceClassRegion.EXT",4,"CLASS: %g %g\t%g %g\t%g %g\t%g %g\t%g EXT\n",
                     source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
@@ -639 +611 @@
             continue;
         }
+
         // Alternate extended source limit calculation
         if (options->altDiffExt) {
@@ -652 +625 @@
                     source->peak->xf, source->peak->yf, Mminor, kMag, dMag, nSigmaMAG, options->sizeLimitCR, options->magLimitCR, options->nSigmaApResid,
                     momentRatioVeres,options->altDiffExtThresh);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
             source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
@@ -968 +942 @@
 }
 
-bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
-#ifdef CHECK_FOOTPRINTS
-    // This gets really expensive for complex images
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-        pmPeak *peak = source->peak;
-        pmFootprint *footprint = peak->footprint;
-        if (!footprint) continue;
-        for (int j = 0; j < footprint->spans->n; j++) {
-            pmSpan *sp = footprint->spans->data[j];
-            psAssert (sp, "missing span");
-        }
-    }
-#endif
-    return true;
-}
-
-/**** ------ old versions of cosmic ray masking ----- ****/
-
-bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
-
-// This attempt to mask the cosmic rays used the isophotal boundary
-bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
-
-    // replace the source flux
-    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
-
-    // flag this as a CR
-    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    pmPeak *peak = source->peak;
-    psAssert (peak, "NULL peak");
-
-    // grab the matching footprint
-    pmFootprint *footprint = peak->footprint;
-    if (!footprint) {
-      psTrace("psphot.czw",2,"Using isophot CR mask code.");
-
-        // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
-        return true;
-    }
-
-    if (!footprint->spans) {
-      psTrace("psphot.czw",2,"Using isophot CR mask code.");
-
-        // if we have no footprint, use the old code to mask by isophot
-        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
-        return true;
-    }
-    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
-    // mask all of the pixels covered by the spans of the footprint
-    for (int j = 1; j < footprint->spans->n; j++) {
-        pmSpan *span1 = footprint->spans->data[j];
-
-        int iy = span1->y;
-        int xs = span1->x0;
-        int xe = span1->x1;
-
-        for (int ix = xs; ix < xe; ix++) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
-    return true;
-}
-
 # define VERBOSE 0
 int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
@@ -1176 +1085 @@
 }
 
-bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
-
-    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    pmPeak *peak = source->peak;
-    psAssert (peak, "NULL peak");
-
-    psImage *mask   = source->maskView;
-    psImage *pixels = source->pixels;
-    psImage *variance = source->variance;
-
-    // XXX This should be a recipe variable
-# define SN_LIMIT 5.0
-
-    int xo = peak->x - pixels->col0;
-    int yo = peak->y - pixels->row0;
-
-    // mark the pixels in this row to the left, then the right
-    for (int ix = xo; ix >= 0; ix--) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
-    }
-    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
-        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
-        if (SN > SN_LIMIT) {
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
-        }
-    }
-
-    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
-    // first go up:
-    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
-    // next go down:
-    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
-        // mark the pixels in this row to the left, then the right
-        for (int ix = 0; ix < pixels->numCols; ix++) {
-            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
-            if (SN < SN_LIMIT) continue;
-
-            bool valid = false;
-            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
-            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
-            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
-
-            if (!valid) continue;
-            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
-    return true;
-}
-
-// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
-// deviation in sigmas.  This is measured on the residual image - should we ignore negative
-// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
-// no longer used by psphot.
-float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
-                         psImageMaskType maskVal, const pmModel *model, float Ro)
-{
-    psF32 *PAR = model->params->data.F32; // Model parameters
-    float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
-
-    // We treat the contour as an ellipse:
-    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
-    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
-    // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
-    // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
-    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
-
-    // min/max value of x is where B^2 - 4AC = 0; solve this for x
-    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
-    if (Q < 0.0) {
-        // ellipse is imaginary
-        return NAN;
-    }
-
-    int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
-    int nPts = 0;                       // Number of points in ellipse
-    float nSigma = 0.0;                 //
-
-    for (int x = -radius; x <= radius; x++) {
-        // Polynomial coefficients
-        // XXX Should we be using the centre of the pixel as x or x+0.5?
-        float A = PS_SQR (1.0 / syy);
-        float B = x * sxy;
-        float C = PS_SQR (x / sxx) - Ro;
-        float T = PS_SQR(B) - 4*A*C;
-        if (T < 0.0) {
-            continue;
-        }
-
-        // y position in source frame
-        float yP = (-B + sqrt (T)) / (2.0 * A);
-        float yM = (-B - sqrt (T)) / (2.0 * A);
-
-        // Get the closest pixel positions (image frame)
-        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
-        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
-        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
-
-        if (xPix < 0 || xPix >= image->numCols) {
-            continue;
-        }
-
-        if (yPixM >= 0 && yPixM < image->numRows &&
-            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
-            nSigma += dSigma;
-            nPts++;
-        }
-
-        if (yPixM == yPixP) {
-            continue;
-        }
-
-        if (yPixP >= 0 && yPixP < image->numRows &&
-            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
-            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
-            nSigma += dSigma;
-            nPts++;
-        }
-    }
-    nSigma /= nPts;
-    return nSigma;
-}
-
-// this was an old attempt to identify cosmic rays based on the peak curvature
-bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
-
-    // classify the sources based on the CR test (place this in a function?)
-    // XXX use an internal flag to mark sources which have already been measured
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // skip source if it was already measured
-        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
-            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
-            continue;
-        }
-
-        // source must have been subtracted
-        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
-            source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED;
-            psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n");
-            continue;
-        }
-
-        psF32 **resid  = source->pixels->data.F32;
-        psF32 **variance = source->variance->data.F32;
-        psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
-
-        // Integer position of peak
-        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
-        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
-
-        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
-        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
-            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
-            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
-            continue;
-        }
-
-        // Skip sources with masked pixels.  These are mostly caught as DEFECT
-        bool keep = true;
-        for (int iy = -1; (iy <= +1) && keep; iy++) {
-            for (int ix = -1; (ix <= +1) && keep; ix++) {
-                if (mask[yPeak+iy][xPeak+ix] & options->maskVal) {
-                    keep = false;
-                }
-            }
-        }
-        if (!keep) {
-            psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n");
-            continue;
-        }
-
-        // Compare the central pixel with those on either side, for the four possible lines through it.
-
-        // Soften variances (add systematic error)
-        float softening = options->soft * PS_SQR(source->peak->rawFlux); // Softening for variances
-
-        // Across the middle: y = 0
-        float cX = 2*resid[yPeak][xPeak]   - resid[yPeak+0][xPeak-1]  - resid[yPeak+0][xPeak+1];
-        float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1];
-        float nX = cX / sqrtf(dcX + softening);
-
-        // Up the centre: x = 0
-        float cY = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak+0]  - resid[yPeak+1][xPeak+0];
-        float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0];
-        float nY = cY / sqrtf(dcY + softening);
-
-        // Diagonal: x = y
-        float cL = 2*resid[yPeak][xPeak]   - resid[yPeak-1][xPeak-1]  - resid[yPeak+1][xPeak+1];
-        float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1];
-        float nL = cL / sqrtf(dcL + softening);
-
-        // Diagonal: x = - y
-        float cR = 2*resid[yPeak][xPeak]   - resid[yPeak+1][xPeak-1]  - resid[yPeak-1][xPeak+1];
-        float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1];
-        float nR = cR / sqrtf(dcR + softening);
-
-        // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2)
-        // Ndof = 4 ? (four measurements, no free parameters)
-        // XXX this value is going to be biased low because of systematic errors.
-        // we need to calibrate it somehow
-        // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq);
-
-        // not strictly accurate: overcounts the chisq contribution from the center pixel (by
-        // factor of 4); also biases a bit low if any pixels are masked
-        // XXX I am not sure I want to keep this value...
-        source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR);
-
-        float fCR = 0.0;
-        int nCR = 0;
-        if (nX > 0.0) {
-            fCR += nX;
-            nCR ++;
-        }
-        if (nY > 0.0) {
-            fCR += nY;
-            nCR ++;
-        }
-        if (nL > 0.0) {
-            fCR += nL;
-            nCR ++;
-        }
-        if (nR > 0.0) {
-            fCR += nR;
-            nCR ++;
-        }
-        source->crNsigma  = (nCR > 0)  ? fCR / nCR : 0.0;
-        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
-
-        if (!isfinite(source->crNsigma)) {
-            continue;
-        }
-
-        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
-        if (source->crNsigma > options->nSigmaCR) {
-            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-            // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
-            // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask);
-        }
-    }
-
-    // now that we have masked pixels associated with CRs, we can grow the mask
-    if (options->grow > 0) {
-        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
-        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
-        psImageConvolveSetThreads(oldThreads);
-        if (!newMask) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
-            return false;
-        }
-        psFree(readout->mask);
-        readout->mask = newMask;
-    }
-    return true;
-}
-
 float psphotSourceSizeFindThreshold (psVector *value, psVector *mask, int maskValue, float minValue, float maxValue, float delta, float guess, float fraction) {
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceStats.c

@@ -58 +58 @@
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
-
-    // XXX TEST:
-    if (detections->allSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->allSources);
-    }
-    if (detections->newSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->newSources);
-    }
 
     // determine the number of allowed threads
@@ -246 +238 @@
     psphotVisualShowMoments (sources);
 
-    // clear the mark bits
-    // psImageMaskPixels (readout->mask, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-    if (detections->allSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->allSources);
-    }
-    if (detections->newSources) {
-        psphotMaskCosmicRayFootprintCheck(detections->newSources);
-    }
-
     return true;
 }
@@ -572 +554 @@
 
             // measure basic source moments (no S/N clipping on input pixels)
+            // sources with (mode & MODE_EXTERNAL) or (mode2 & MODE2_MATCHED) use the 
+            // supplied Mx,My value for the centroid (not recalculated)
             status = pmSourceMoments (source, 4*sigma[i], sigma[i], 0.0, 0.0, maskVal);
         }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotStackImageLoop.c

@@ -22 +22 @@
 
     pmFPAview *view = pmFPAviewAlloc (0);
-    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
-    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
-    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
 
     if (!input) {
@@ -32 +30 @@
 
     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    if (useRaw && inputRaw == NULL) {
-        psLogMsg ("psphot", 1, "PSPHOT.STACK.USE.RAW set but no raw input.");
-        useRaw = false;
-    }
-
-    bool radial_apertures = psMetadataLookupBool(NULL, recipe, "RADIAL_APERTURES");
-    bool match_psfs = radial_apertures;
-    bool needConvolved = radial_apertures || !useRaw;
-    if (!needConvolved) {
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.INPUT.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.MASK.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.VARIANCE.CNV");
-        pmFPAfileActivate (config->files, false, "PSPHOT.STACK.PSF.CNV");
-    }
 
     // just load the full set of images up front except for EXPNUM which we defer
     pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.RAW");
-    pmFPAfileActivate (config->files, false, "PSPHOT.STACK.EXPNUM.CNV");
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
 
@@ -71 +53 @@
                 psMemDump("load");
 
-                if (match_psfs) {
-                    // Generate the 1st PSF-matched image set (larger target PSFs are generated by smoothing this image)
-                    if (!psphotStackMatchPSFs (config, view)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
-                } else {
-                    if (!psphotStackAllocateOutput (config, view, recipe)) {
-                        psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
-                        psFree (view);
-                        return false;
-                    }
+                if (!psphotStackAllocateOutput (config, view, recipe)) {
+                    psError(psErrorCodeLast(), false, "failure in psphotStackAllocateOutput for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
+                    psFree (view);
+                    return false;
                 }
                 psMemDump("stackmatch");
@@ -126 +99 @@
 
     // Load the appropriate EXPNUM image
-    pmFPAfileActivate (config->files, true, useRaw ? "PSPHOT.STACK.EXPNUM.RAW" : "PSPHOT.STACK.EXPNUM.CNV");
+    pmFPAfileActivate (config->files, true, "PSPHOT.STACK.EXPNUM.RAW");
 
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa EXPNUM in psphot.");
@@ -146 +119 @@
                 if (! readout->data_exists) { continue; }
 
-                if (!psphotSetNFrames (config, view, useRaw ? inputRaw->name : inputCnv->name)) ESCAPE ("failed to setNFrames.");
+                if (!psphotSetNFrames (config, view, input->name)) ESCAPE ("failed to setNFrames.");
             }
         }
@@ -167 +140 @@
 bool UpdateHeadersForReadout (pmConfig *config, pmFPAview *view) {
 
-    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-
     int num = psphotFileruleCount(config, "PSPHOT.INPUT");
 
@@ -179 +149 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = useRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 
@@ -205 +173 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 
@@ -242 +208 @@
         psAssert (output, "missing file?");
 
-        pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
-        pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
-        pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+        pmFPAfile *input = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
         psAssert (input, "missing input file");
 

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotStackMatchPSFsNext.c

@@ -43 +43 @@
     for (int i = 0; i < num; i++) {
         if (!psphotStackMatchPSFsNextReadout (config, view, filerule, i, lastSize)) {
-            psError (PSPHOT_ERR_CONFIG, false, "failed to smooth image %s (%d) to target PSF", filerule, i);
+            psLogMsg ("psphot", PS_LOG_INFO, "failed to smooth image %s (%d) to target PSF", filerule, i);
             psImageConvolveSetThreads(oldThreads);
             return false;
@@ -102 +102 @@
 
     if (targetFWHM <= currentFWHM) {
-        psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+        // psError (PSPHOT_ERR_CONFIG, true, "target FWHM cannot be smaller than current FWHM");
+        psLogMsg ("psphot", PS_LOG_INFO, "target FWHM (%f) is smaller than current FWHM (%f), not smoothing\n", targetFWHM, currentFWHM);
+        fwhmValues->data.F32[lastSize + 1] = currentFWHM;
+        pmReadoutMaskInvalid(readout, maskVal, maskSat);
         return false;
     }

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotStackMatchPSFsUtils.c

@@ -325 +325 @@
     if (deconv > deconvLimit) {
         rejectReadout = true;
-        psLogMsg("psphot", PS_LOG_WARN, "Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image for matched psf analysis%d\n", deconv, deconvLimit, index);
+        psLogMsg("psphot", PS_LOG_WARN, "Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting input %d for matched psf analysis\n", deconv, deconvLimit, index);
     }
     

branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotStackReadout.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+static bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF);
+static bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *filerule, const char *fPSF, int index);
 static bool psphotStackLoadWCS(pmConfig *config, const pmFPAview *view, const char *filerule);
 static void logMemStats(const char *heading);
 
-// we have 3 possible real filesets:
+// relevant filesets:
 # define STACK_RAW "PSPHOT.STACK.INPUT.RAW"
-# define STACK_CNV "PSPHOT.STACK.INPUT.CNV"
-# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"  /* the psf-matched image */
-
-// we have 3 files on which we operate:
-// DET (detection image)       : nominally RAW (optionally CNV?)
-// SRC (source analysis image) : nominally CNV (optionally RAW)
-// OUT (psf-matched images)    : always OUT
+# define STACK_OUT "PSPHOT.STACK.OUTPUT.IMAGE"
+
+// XXX STACK_OUT currently is a copy of STACK_RAW, but should be a pointer to is as in psphot (single)
+
+// TEST CODE, can be removed
+bool psphotDumpImages (pmConfig *config, const pmFPAview *view, const char *filerule, char *base) {
+
+    // XXX do nothing
+    return true;
+
+    int num = psphotFileruleCount(config, "PSPHOT.INPUT");
+
+    for (int i = 0; i < num; i++) {
+        // 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?");
+
+        char line[256];
+        snprintf (line, 256, "%s.%d.im.fits", base, i);
+        psphotSaveImage (NULL, readout->image, line);
+
+        snprintf (line, 256, "%s.%d.wt.fits", base, i);
+        psphotSaveImage (NULL, readout->variance, line);
+
+        snprintf (line, 256, "%s.%d.mk.fits", base, i);
+        psphotSaveImage (NULL, readout->mask, line);
+    }
+    return true;
+}
 
 bool psphotStackVisualFilerule(pmConfig *config, const pmFPAview *view, const char *filerule) {
@@ -67 +94 @@
     psAssert (breakPt, "configuration error: set BREAK_POINT");
 
-    // we have 3 relevant files: RAW (unconvolved), CNV (convolved stack), OUT (psf-matched stack)
-    // select which image (RAW or CNV) is used for analysis (RAW always used for detection)
-    bool useRaw = psMetadataLookupBool (NULL, recipe, "PSPHOT.STACK.USE.RAW");
-    char *STACK_SRC = useRaw ? STACK_RAW : STACK_CNV;
-    char *STACK_DET = STACK_RAW;
-
     // load WCS 
-    if (!psphotStackLoadWCS(config, view, STACK_SRC)) {
-        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_SRC);
+    if (!psphotStackLoadWCS(config, view, STACK_RAW)) {
+        psError (PSPHOT_ERR_CONFIG, false, "trouble loading WCS for %s", STACK_RAW);
         return false;
     }
 
     // set the photcode for each image
-    if (!psphotAddPhotcode (config, view, STACK_SRC)) {
+    if (!psphotAddPhotcode (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode");
         return false;
@@ -86 +107 @@
 
     // Generate the mask and weight images (if not supplied) and set mask bits. 
-    // This also insures that all invalid pixels are masked (this is done for STACK_CNV in psphotStackMatchPSFs)
-    if (!psphotSetMaskAndVariance (config, view, STACK_DET)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSetMaskAndVariance (config, view, STACK_OUT)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+    if (!psphotSetMaskAndVariance (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "NOTHING")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // generate a background model (median, smoothed image)
-    if (!psphotModelBackground (config, view, STACK_DET)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (!psphotSubtractBackground (config, view, STACK_DET)) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    if (strcmp(STACK_SRC, STACK_DET)) {
-#define MODEL_BACKGROUND_SRC 1
-#ifdef MODEL_BACKGROUND_SRC
-        // work around the fact that the background levels on the convolved
-        // and unconvolved stacks can be different
-        if (!psphotModelBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
-#endif
-        if (!psphotSubtractBackground (config, view, STACK_SRC)) {
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+    if (!psphotModelBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_RAW)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "BACKMDL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
 #ifdef MAKE_CHISQ_IMAGE
     // also make the chisq detection image
-    if (!psphotStackChisqImage(config, view, STACK_DET, STACK_SRC)) {
+    if (!psphotStackChisqImage(config, view, STACK_RAW, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    if (!strcasecmp (breakPt, "CHISQ")) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 #endif
-    if (!strcasecmp (breakPt, "CHISQ")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
 
     // find the detections (by peak and/or footprint) in the image.
     // This finds the detections on Chisq image as well as the individuals
-    if (!psphotFindDetections (config, view, STACK_DET, true)) { // pass 1
-        // this only happens if we had an error in psphotFindDetections
+    if (!psphotFindDetections (config, view, STACK_RAW, true)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
-    // is just a copy of the container pointer; the sources on both DET and SRC are the same
-    // memory objects
-    if (strcmp(STACK_SRC, STACK_DET)) {
-        if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
-            psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-            return psphotReadoutCleanup (config, view, STACK_SRC);
-        }
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // construct sources and measure basic stats (saved on detections->newSources)
-    if (!psphotSourceStats (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotSourceStats (config, view, STACK_RAW, true)) { // pass 1
         psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PEAKS")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-    // psphotDumpTest (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+    // psphotDumpTest (config, view, STACK_RAW);
     psMemDump("sourcestats");
     logMemStats("sourcestats");
@@ -164 +157 @@
     // classify sources based on moments, brightness
     // only run this on detections from the input images, not chisq image
-    if (!psphotRoughClass (config, view, STACK_SRC)) {
+    if (!psphotRoughClass (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-    }
-
-    // If DET and SRC are different images, subtract radial profiles for the convolved
-    // image first.  The profiles found on the convolved image will be replaced by those
-    // found for the unconvolved image.  Downstream, in psphotFindDetections, we will
-    // replace the the profiles (and re-subtract them) for the detection image, so we want
-    // to keep those versions of the profiles on the sources.
-    if (strcmp(STACK_SRC, STACK_DET)) {
-      // find and subtract radial profile models for saturated stars (XXX change name eventually)
-      if (!psphotDeblendSatstars (config, view, STACK_SRC)) {
+        return psphotReadoutCleanup (config, view, STACK_RAW);
+    }
+
+    // find and subtract radial profile models for saturated stars (XXX change name eventually)
+    if (!psphotDeblendSatstars (config, view, STACK_RAW)) {
         psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
-      }
-    }
-    // find and subtract radial profile models for saturated stars (XXX change name eventually)
-    if (!psphotDeblendSatstars (config, view, STACK_DET)) {
-        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     // only run this on detections from the input images, not chisq image
-    if (!psphotImageQuality (config, view, STACK_SRC)) { // pass 1
+    if (!psphotImageQuality (config, view, STACK_RAW)) { // pass 1
         psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "MOMENTS")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // use bright stellar objects to measure PSF
-    if (!psphotChoosePSF (config, view, STACK_SRC, true)) { // pass 1
+    if (!psphotChoosePSF (config, view, STACK_RAW, true)) { // pass 1
         psLogMsg ("psphot", 3, "failure to construct a psf model");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
     if (!strcasecmp (breakPt, "PSFMODEL")) {
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // merge the newly selected sources into the existing list
     // NOTE: merge OLD and NEW
-    psphotMergeSources (config, view, STACK_SRC);
+    psphotMergeSources (config, view, STACK_RAW);
 
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    psphotFitSourcesLinear (config, view, STACK_SRC, false, false);
-    psphotStackVisualFilerule(config, view, STACK_SRC);
+    psphotFitSourcesLinear (config, view, STACK_RAW, false, false);
+    psphotStackVisualFilerule(config, view, STACK_RAW);
 
     // measure the radial profiles to the sky
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted.  this pass starts with a circular
@@ -225 +206 @@
     // but iterates to an appropriately larger size
     logMemStats("before.kron.1");
-    psphotKronIterate(config, view, STACK_SRC, 1);
+    psphotKronIterate(config, view, STACK_RAW, 1);
     logMemStats("after.kron.1");
         
     // identify CRs and extended sources
-    psphotSourceSize (config, view, STACK_SRC, true);
+    psphotSourceSize (config, view, STACK_RAW, true);
 
     // non-linear PSF and EXT fit to brighter sources
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
-    psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotBlendFit (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // replace all sources (do NOT ignore subtraction state)
-    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 1 (detections->allSources)
 
     logMemStats("pass1");
@@ -245 +226 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
-    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 2 (detections->allSources)
+    // NOTE 2 : this function subtracts the models from the given filerule
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 2 (detections->allSources)
 
     // NOTE: possibly re-measure background model here with objects subtracted / or masked
@@ -252 +233 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
-        // if DET and SRC are different images, generate children sources for all sources in
-        // the SRC image.  This operation replaces the existing DETECTION container on DET
-        // which is currently a view to the one on SRC).  children sources go to
-        // det->allSources
-        if (strcmp(STACK_SRC, STACK_DET)) {
-            psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
-
-            //  subtract all sources from DET (this will subtract using the psf model for SRC, which
-            //  will somewhat oversubtract the sources -- this is OK
-            psphotRemoveAllSources (config, view, STACK_DET, false); // do not ignore subtraction state for sources
-        }
-
         // add noise for subtracted objects
-        psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+        psphotAddNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
         // find fainter sources
         // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
-        psphotFindDetections (config, view, STACK_DET, false); // pass 2 (detections->peaks, detections->footprints)
+        psphotFindDetections (config, view, STACK_RAW, false); // pass 2 (detections->peaks, detections->footprints)
 
         // remove noise for subtracted objects (ie, return to normal noise level)
@@ -276 +245 @@
         bool footprintsUseUnsubtracted = psMetadataLookupBool(NULL, recipe, "FOOTPRINT_USE_UNSUBTRACTED");
         if (!footprintsUseUnsubtracted) {
-            psphotSubNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
+            psphotSubNoise (config, view, STACK_RAW); // pass 1 (detections->allSources)
         }
-
-        // if DET and SRC are different images, copy the detections from DET to SRC 
-        // (this operation just ensures the metadata container has a view on SRC as well
-        if (strcmp(STACK_SRC, STACK_DET)) {
-            // replace all sources in DET
-            psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
-
-            // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
-            if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
 
         // define new sources based on only the new peaks & measure moments
         // NOTE: new sources are saved on detections->newSources
-        psphotSourceStats (config, view, STACK_SRC, false); // pass 2 (detections->newSources)
+        psphotSourceStats (config, view, STACK_RAW, false); // pass 2 (detections->newSources)
 
         // set source type
         // NOTE: apply only to detections->newSources
-        if (!psphotRoughClass (config, view, STACK_SRC)) { // pass 2 (detections->newSources)
+        if (!psphotRoughClass (config, view, STACK_RAW)) { // pass 2 (detections->newSources)
             psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
-            return psphotReadoutCleanup (config, view, STACK_SRC);
+            return psphotReadoutCleanup (config, view, STACK_RAW);
         }
 
         // replace all sources so fit below applies to all at once
         // NOTE: apply only to OLD sources (which have been subtracted)
-        psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
+        psphotReplaceAllSources (config, view, STACK_RAW, false); // pass 2
 
         // merge the newly selected sources into the existing list
         // NOTE: merge OLD and NEW
         // XXX check on free of sources...
-        psphotMergeSources (config, view, STACK_SRC); // (detections->newSources + detections->allSources -> detections->allSources)
+        psphotMergeSources (config, view, STACK_RAW); // (detections->newSources + detections->allSources -> detections->allSources)
 
         // Construct an initial model for each object, set the radius to fitRadius, set circular
         // fit mask.  NOTE: only applied to sources without guess models
-        psphotGuessModels (config, view, STACK_SRC);
+        psphotGuessModels (config, view, STACK_RAW);
     }
 
@@ -325 +281 @@
     if (splitLinearFit) {
         psLogMsg ("psphot", 3, "splitting fit of detected and matched soures\n");
-        // Fit the detected sources separately from matched that wea are about to create.
+        // Fit the detected sources separately from matched ones that we are about to create.
         // NOTE: apply to ALL sources but only include sources with postitive flux in the fit
-        psphotFitSourcesLinear (config, view, STACK_SRC, true, true); // pass 3 (detections->allSources)
+        psphotFitSourcesLinear (config, view, STACK_RAW, true, true); // pass 3 (detections->allSources)
     }
 
@@ -335 +291 @@
     // this just match the detections for the chisq image, and not bother measuring the source
     // stats in that case...?
-    objects = psphotMatchSources (config, view, STACK_SRC);
+    objects = psphotMatchSources (config, view, STACK_RAW);
     psMemDump("matchsources");
 
@@ -344 +300 @@
     // Construct an initial model for each object, set the radius to fitRadius, set circular
     // fit mask.  NOTE: only applied to sources without guess models
-    psphotGuessModels (config, view, STACK_SRC);
+    psphotGuessModels (config, view, STACK_RAW);
 
     psphotStackObjectsUnifyPosition (objects);
 
-    psphotStackObjectsSelectForAnalysis (config, view, STACK_SRC, objects);
+    // psphotStackObjectsSelectForAnalysis (config, view, STACK_RAW, objects);
 
     // final linear fit. NOTE: if splitLinearFit is true above, this pass will only fit
     // the unsubtracted (matched) sources (the sources that we fit above are subtracted)
-    psphotFitSourcesLinear (config, view, STACK_SRC, true, false); // pass 4 (detections->allSources)
+    psphotFitSourcesLinear (config, view, STACK_RAW, true, false); // pass 4 (detections->allSources)
 
     // measure the radial profiles to the sky (only measures new objects)
-    psphotRadialProfileWings (config, view, STACK_SRC);
+    psphotRadialProfileWings (config, view, STACK_RAW);
 
     // re-measure the kron mags with models subtracted
     // psphotKronMasked(config, view, STACK_SRC);
     logMemStats("before.kron.2");
-    psphotKronIterate(config, view, STACK_SRC, 2);
+    psphotKronIterate(config, view, STACK_RAW, 2);
     logMemStats("after.kron.2");
 
     // measure source size for the remaining sources
     // NOTE: applies only to NEW (unmeasured) sources
-    psphotSourceSize (config, view, STACK_SRC, false); // pass 2 (detections->allSources)
+    psphotSourceSize (config, view, STACK_RAW, false); // pass 2 (detections->allSources)
 
     psMemDump("psfstats");
@@ -371 +327 @@
     // drop matched sources without any useful measurements and set kron radii for the ones
     // we decide to keep
-    psphotFilterMatchedSources (config, view, STACK_SRC, objects);
+    psphotFilterMatchedSources (config, view, STACK_RAW, objects);
 
     // measure kron fluxes for the matched sources only
-    psphotKronIterate(config, view, STACK_SRC, 3);
+    psphotKronIterate(config, view, STACK_RAW, 3);
+
+    // decide which source(s) are to be fitted with the extended source analysis code.
+    psphotChooseAnalysisOptionsByObject (config, view, STACK_RAW, objects);
 
     // measure elliptical apertures, petrosians (objects sorted by S/N)
     // psphotExtendedSourceAnalysisByObject (config, objects, view, STACK_SRC); // pass 1 (detections->allSources)
-    psphotExtendedSourceAnalysis (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceAnalysis (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // measure non-linear extended source models (exponential, deVaucouleur, Sersic) (sources sorted by S/N)
-    psphotExtendedSourceFits (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    psphotExtendedSourceFits (config, view, STACK_RAW); // pass 1 (detections->allSources)
 
     // create source children for the OUT filerule (for radial aperture photometry and output) 
-    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects);
+    // NOTE: The new source children have image arrays pointing to the readout associated with 
+    // STACK_OUT.  in psphotStackMatchPSFsetup, we copy the current pixel values from RAW to OUT, 
+    // but keep the pointers the same so we do not break these source image references
+    // XXX NOTE : if we use the pre-20130914 psphotStackReadout code, we need to use 'false' for the
+    // sourcesSubtracted argument
+    psArray *objectsOut = psphotSourceChildrenByObject (config, view, STACK_OUT, objects, true);
     if (!objectsOut) {
         psFree(objects);
         psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
@@ -396 +360 @@
         // this forces photometry on the undetected sources from other images
 
-        // NOTE: we always do the radial apertures analysis on the convolved image since
-        // those are the ones that are psf matched and are the source of STACK_OUT's pixels
-        // XXX: Actually if PSPHOT.STACK.MATCH.PSF.SOURCE were set to RAW this wouldn't be true.
-        // but in that case we don't get past the psf matching step because there is no
-        // target psf for the RAW inputs
-
-        // If useRaw copy the sources to the convolved readout
-        if (strcmp(STACK_SRC, STACK_CNV)) {
-            if (!psphotCopySources (config, view, STACK_CNV, STACK_SRC)) {
-                psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
-                return psphotReadoutCleanup (config, view, STACK_SRC);
-            }
-        }
-        // mark any inputs that we want to skip the matched apertures for
-        psphotStackSetInputsToSkip(config, view, STACK_CNV, true);
-        psphotStackSetInputsToSkip(config, view, STACK_OUT, true);
-        psphotRadialApertures (config, view, STACK_CNV, 0); // entry 0 == unmatched
-        psMemDump("extmeas");
-
+        // set up the FWHM vector
+        psphotStackMatchPSFsetup (config, view, STACK_OUT, STACK_RAW);
+        psphotDumpImages (config, view, STACK_RAW, "raw.t0");
+        psphotDumpImages (config, view, STACK_OUT, "out.t0");
 
         int nRadialEntries = psphotStackMatchPSFsEntries(config, view, STACK_OUT);
-        for (int entry = 1; entry < nRadialEntries; entry++) {
+
+        for (int entry = 0; entry < nRadialEntries; entry++) {
             // NOTE: entry 0 is the unmatched image set
 
-            // re-measure the PSF for the smoothed image (using entries in 'allSources')
-            psphotChoosePSF (config, view, STACK_OUT, false);
-
-            // this is necessary to update the models based on the new PSF
-            psphotResetModels (config, view, STACK_OUT);
-
-            // this is necessary to get the right normalization for the new models
-            psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+            char line[256];
 
             // measure circular, radial apertures (objects sorted by S/N)
             psphotRadialApertures (config, view, STACK_OUT, entry); 
+            snprintf (line, 256, "%s.%d", "out.t1", entry);
+            psphotDumpImages (config, view, STACK_OUT, line);
 
             // replace the flux in the image so it is returned to its original state
             psphotReplaceAllSources (config, view, STACK_OUT, false);
-
-            // smooth to the next FWHM, or set 'smoothAgain' to false if no more 
-            psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
-            psMemDump("matched");
+            snprintf (line, 256, "%s.%d", "out.t2", entry);
+            psphotDumpImages (config, view, STACK_OUT, line);
+
+            if (entry < nRadialEntries - 1) {
+                // smooth to the next FWHM
+                // this function does nothing if the targetFWHM is smaller than the currentFWHM
+                psphotStackMatchPSFsNext(config, view, STACK_OUT, entry);
+                snprintf (line, 256, "%s.%d", "out.t3", entry);
+                psphotDumpImages (config, view, STACK_OUT, line);
+                psMemDump("matched");
+
+                // re-measure the PSF for the smoothed image (using entries in 'allSources')
+                if (!psphotChoosePSF (config, view, STACK_OUT, false)) {
+                    psLogMsg ("psphot", 3, "failure to construct a psf model in radial aperture loop for entry :%d", entry);
+                    return psphotReadoutCleanup (config, view, STACK_RAW);
+                }
+
+                // this is necessary to update the models based on the new PSF
+                psphotResetModels (config, view, STACK_OUT);
+
+                // this is necessary to get the right normalization for the new models
+                // and to subtract the sources
+                psphotFitSourcesLinear (config, view, STACK_OUT, false, false);
+                snprintf (line, 256, "%s.%d", "out.t4", entry);
+                psphotDumpImages (config, view, STACK_OUT, line);
+            }
         }
     }
-    psphotStackSetInputsToSkip(config, view, STACK_CNV, false);
-    psphotStackSetInputsToSkip(config, view, STACK_OUT, false);
 
     // measure aperture photometry corrections
-    if (!psphotApResid (config, view, STACK_SRC)) {
+    if (!psphotApResid (config, view, STACK_RAW)) {
         psFree (objects);
         psFree (objectsOut);
         psLogMsg ("psphot", 3, "failed on psphotApResid");
-        return psphotReadoutCleanup (config, view, STACK_SRC);
+        return psphotReadoutCleanup (config, view, STACK_RAW);
     }
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, STACK_SRC);
-
-    if (!useRaw) {
-        // psphotEfficiency wants to have the PSF of the image, but since we are measuring on
-        // the convolved images we need to generate PSFs for the DET images
-        if (!psphotChoosePSF (config, view, STACK_DET, false)) {
-            psLogMsg ("psphot", 3, "failure to construct a psf model for raw input");
-            return psphotReadoutCleanup (config, view, STACK_DET);
-        }
-    }
-    if (!psphotEfficiency(config, view, STACK_DET)) {
+    psphotMagnitudes(config, view, STACK_RAW);
+
+    // calculate lensing parameters
+    if (!psphotLensing(config, view, STACK_RAW)) {
+        psErrorStackPrint(stderr, "Unable to do lensing parameters.");
+        psErrorClear();
+    }
+
+    if (!psphotEfficiency(config, view, STACK_RAW)) {
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
         psErrorClear();
     }
-    psphotCopyEfficiency (config, view, STACK_OUT, STACK_DET);
+    psphotCopyEfficiency (config, view, STACK_OUT, STACK_RAW);
 
     logMemStats("final");
 #if (1)
-    psphotSourceMemory(config, view, STACK_SRC);
+    psphotSourceMemory(config, view, STACK_RAW);
     psphotSourceMemory(config, view, STACK_OUT);
 #endif
 
-    // replace failed sources?
-    // psphotReplaceUnfitSources (sources);
-
     // replace background in residual image
-    psphotSkyReplace (config, view, STACK_DET);
+    psphotSkyReplace (config, view, STACK_RAW);
 
     // drop the references to the image pixels held by each source
+    psphotSourceFreePixels (config, view, STACK_RAW);
     psphotSourceFreePixels (config, view, STACK_OUT);
-    psphotSourceFreePixels (config, view, STACK_SRC);
 
 #ifdef MAKE_CHISQ_IMAGE
     // remove chisq image from config->file:PSPHOT.INPUT
-    psphotStackRemoveChisqFromInputs(config, STACK_DET);
-    if (strcmp(STACK_SRC, STACK_DET)) {
-        psphotStackRemoveChisqFromInputs(config, STACK_SRC);
-    }
+    psphotStackRemoveChisqFromInputs(config, STACK_RAW);
 #endif
 
@@ -496 +453 @@
 
     // create the exported-metadata and free local data
-    return psphotReadoutCleanup (config, view, STACK_SRC);
+    return psphotReadoutCleanup (config, view, STACK_RAW);
 }
 
@@ -552 +509 @@
 }
 
-
-
 /* here is the process:
 
- * we have three(*) images:
- * RAW : unconvolved image stack 
- * CNV : input convolved image
+ * we have two image sets:
+ * RAW : unconvolved image stacks
 
  * OUT : psf-matched output image (there may be more than one of
@@ -649 +603 @@
 
    */
+
+
+// generate a vector fwhmValues where the first has the fwhm of the raw image and the
+// successive entries have the target values
+
+bool psphotStackMatchPSFsetup (pmConfig *config, const pmFPAview *view, const char *filerule, const char *filerulePSF) {
+
+    bool status;
+
+    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
+
+        if (!psphotStackMatchPSFsetupReadout (config, view, filerule, filerulePSF, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to define target PSF sizes");
+            return false;
+        }
+    }
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index);
+
+// copy the pixels from RAW to OUT (
+
+bool psphotStackMatchPSFsetupReadout (pmConfig *config, const pmFPAview *view, const char *fileruleOut, const char *fileruleRaw, int index) {
+
+    bool status;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, fileruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // find the currently selected readout
+    pmFPAfile *fileRaw = pmFPAfileSelectSingle(config->files, fileruleRaw, index); // File of interest
+    psAssert (fileRaw, "missing file?");
+
+    pmReadout *readoutRaw = pmFPAviewThisReadout(view, fileRaw->fpa);
+    psAssert (readoutRaw, "missing readout?");
+
+    readoutOut->image = psImageCopy(readoutOut->image, readoutRaw->image, PS_TYPE_F32);
+    if (readoutRaw->variance) {
+        readoutOut->variance = psImageCopy(readoutOut->variance, readoutRaw->variance, PS_TYPE_F32);
+    }
+    if (readoutRaw->mask) {
+        readoutOut->mask = psImageCopy(readoutOut->mask, readoutRaw->mask, PS_TYPE_IMAGE_MASK);
+    }
+
+    // pmChip *chipRaw = pmFPAviewThisChip(view, fileRaw->fpa); // The chip holds the PSF
+    // psAssert (chipRaw, "missing chip");
+
+    pmPSF *psf = psMetadataLookupPtr(&status, readoutRaw->analysis, "PSPHOT.PSF"); // PSF
+    if (!psf) {
+        // we should have a PSF by this point in psphot
+        psError(PSPHOT_ERR_PROG, true, "Unable to find PSF.");
+        return false;
+    }
+
+    float fwhmRaw = psphotPSFseeing (psf, readoutRaw, index);
+
+    psVector *fwhmValues = psVectorAllocEmpty(10, PS_TYPE_F32);
+    psVectorAppend(fwhmValues, fwhmRaw);
+
+    // is a single target FWHM specified, or a set of values?  set up the vector options->targetSeeing and the local 1st value
+    float targetSeeing = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+    if (!status) {
+        psVector *targetSeeing = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+        psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+        for (int i = 0; i < targetSeeing->n; i++) {
+            psVectorAppend(fwhmValues, targetSeeing->data.F32[i]);
+        }           
+    } else {
+        psVectorAppend(fwhmValues, targetSeeing);
+    }
+
+    psMetadataAddVector(readoutOut->analysis, PS_LIST_TAIL, "STACK.PSF.FWHM.VALUES", PS_META_REPLACE, "PSF sizes", fwhmValues);
+    psFree (fwhmValues);
+
+    return true;
+}
+
+float psphotPSFseeing (pmPSF *psf, pmReadout *readout, int index) {
+
+    psImage *image = readout->image;
+
+    int Nx = image->numCols;
+    int Ny = image->numRows;
+
+    float sumFWHM = 0.0;                  // FWHM for image
+    int numFWHM = 0;                      // Number of FWHM measurements
+    for (float x = 0; x < Nx; x += 0.25*Nx) {
+        for (float y = 0; y < Ny; y += 0.25*Ny) {
+            float fwhm = pmPSFtoFWHM(psf, x, y);
+            if (isfinite(fwhm)) {
+                sumFWHM += fwhm;
+                numFWHM++;
+            }
+        }
+    }
+    if (numFWHM == 0) {
+        psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", index);
+        return NAN;
+    } 
+
+    float fwhm = sumFWHM / (float) numFWHM;
+    psLogMsg ("psphotStack", PS_LOG_INFO, "Input Seeing for %d: %f\n", index, fwhm);
+
+    return fwhm;
+}

branches/eam_branches/ps2-tc3-20130727/psphot/test/tap_psphot_galaxygrid.pro

@@ -28 +28 @@
 $RefOptions = $RefOptions -Df STARS.DENSITY 10.0
 $RefOptions = $RefOptions -Df STARS.SIGMA.LIM 0.5
-
+$RefOptions = $RefOptions -nx 3000 -ny 3000
+
+if (not($?PSFMODEL))
+  $PSFMODEL = PS1_V1
+end
+
+macro reset.options
 # options for the simulated images (using the refimage for the stars)
 $FakeOptions = $BaseOptions
 $FakeOptions = $FakeOptions -exptime 30.0
-$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+# $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+$FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_$PSFMODEL
+$FakeOptions = $FakeOptions -nx 3000 -ny 3000
   
 # sample alternate options:
@@ -65 +73 @@
 $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
 $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
-$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
-$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
+$FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+$FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+end
+
+if (not($?FakeConfig)) reset.options
 
 list fwhm 
@@ -75 +86 @@
 end
 
-macro go.grid.set.devexp
+if (not($?CONVOLVE_NSIGMA)) set CONVOLVE_NSIGMA = 5.0
+
+# generate fake images and run psphot on them
+macro mkexp.devexp.single
+  if ($0 != 5)
+    echo "USAGE: mkexp.devexp.single (basename) (type) (Rmajor) (fwhm)"  
+    break
+  end
+
+  $basename = $1
+  $type = $2
+  $Rmajor = $3
+  $fwhm = $4
+
+  $Aratio = 1.0
+
   $FakeConfig = -camera SIMTEST
   $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
@@ -84 +110 @@
   $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
   $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
+  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
   $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
   $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 
   $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
-  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
-  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
   $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
   $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
+  $FakeConfig = $FakeConfig -Df CONVOLVE.NSIGMA $CONVOLVE_NSIGMA 
+  $BaseConfig = $FakeConfig
+
+  $FakeConfig = $BaseConfig
+  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
+  $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
+  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
+  $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
+          
+  mkexp $basename $fwhm $type
+end
+
+# generate fake images and run psphot on them
+macro normtest.mkexp.devexp
+  $FakeConfig = -camera SIMTEST
+  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
+  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
+  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
+  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
+  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
+  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
+  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
+  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
+  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
+  $BaseConfig = $FakeConfig
 
   # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
@@ -103 +161 @@
   # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
 
+  mkdir normtest
+
   $Nseq = 0
   foreach type EXP DEV
     foreach Rmajor 3 10 30
       foreach Aratio 0.25 0.5 1.0
-        foreach fwhm 0.8 1.0 1.5
+        foreach fwhm 1.0
+          $FakeConfig = $BaseConfig
           $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
           $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
@@ -113 +174 @@
           $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
           
-          sprint name "sample.%02d" $Nseq
+          sprint name "normtest/test.%02d" $Nseq
           mkexp $name $fwhm $type
-          fitexp $name $name.fit $type\_CONV
           $Nseq ++
         end
@@ -123 +183 @@
 end
 
-macro go.grid.set.sersic
+# generate fake images and run psphot on them
+macro grid.mkexp.devexp
+  if ($0 != 2)
+    echo "USAGE: grid.mkexp.devexp (dir)"
+    break
+  end
+
+  mkdir $1
+
   $FakeConfig = -camera SIMTEST
   $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
@@ -132 +200 @@
   $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
   $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
+  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
   $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
   $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 
   $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
-  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 120
-  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 120
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
+  $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
+  $BaseConfig = $FakeConfig
+
+  $Nseq = 0
+  foreach type EXP DEV
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          $FakeConfig = $BaseConfig
+          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
+          $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
+          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN $Aratio
+          $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX $Aratio
+          
+          sprint name "$1/sample.%02d" $Nseq
+          mkexp $name $fwhm $type
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+# generate fake images and run psphot on them
+macro grid.fitexp.devexp
+  if ($0 != 3)
+    echo "USAGE: grid.fitexp.devexp (srcdir) (outdir)"
+    break
+  end
+
+   #$Nseq = 0
+   #foreach type EXP DEV
+
+  mkdir $2
+
+  $Nseq = 0
+  foreach type EXP DEV
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          sprint srcname "$1/sample.%02d" $Nseq
+          sprint outname "$2/sample.%02d.fit" $Nseq
+          fitexp $srcname $outname $type\_CONV
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+macro grid.load.devexp
+  if ($0 != 3)
+    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
+    break
+  end
+
+  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
+  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
+  delete -q min_S Min_S
+
+  $Nseq = 0
+  foreach type EXP DEV
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          sprint name "sample.%02d" $Nseq
+          cmf.load.concat $1/$name.dat $2/$name.fit.cmf $type
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+macro grid.plots.devexp
+  if ($0 != 2)
+    echo "USAGE: grid.plot.devexp (version)"
+    break
+  end
+  # things to examine: theta, Rmajor, AR
+
+  # go.grid.check.devexp
+
+  # check on theta
+  set dT = Tot_s - Tin_s
+  set ARin = rin_s / Rin_s
+  # lim ARin dT; clear; box; plot ARin dT
+
+  subset dTx = dT if (ARin < 0.95)
+  histogram dTx NdT -8 8 0.1 -range dTi
+  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
+  label -x "angle offset (degrees, only non-circular)" -y "number count"
+  resize 700 320
+
+  # check on Rmajor
+  set dR = Rin_s - Rot_s
+  # lim Rin_s dR; clear; box; plot Rin_s dR
+  # lim Rot_s dR; clear; box; plot Rot_s dR
+  # lim Rin_s dR; clear; box; plot Rin_s dR
+  create n 0 dR[]
+  set dRf = dR / Rin_s
+  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
+  label -x sequence -y "1 - R_out| / R_in|"
+  resize 700 320
+
+  # check on A.Ratio
+  set ARot = rot_s / Rot_s
+  # lim ARin ARot; clear; box; plot ARin ARot
+  set fAR = ARot / ARin
+  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
+  label -x sequence -y "AR_out| / AR_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dM = Mot_s - Min_s
+  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM    
+  label -x sequence -y "M_out| - M_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dI = Iot_s - Iin_s
+  lim -n 4$1 n -1.0 1.0; clear; box; plot n dI
+  label -x sequence -y "I_out| - I_in|"
+  resize 700 320
+end
+
+macro grid.plot.stars
+  if ($0 != 2)
+    echo "USAGE: grid.plot.stars (version)"
+    break
+  end
+  # things to examine: theta, Rmajor, AR
+
+  # go.grid.check.devexp
+
+  # check on position
+  set dX = Xot_s - int(Xin_s) - 0.5
+  set dY = Yot_s - int(Yin_s) - 0.5
+  set ARin = rin_s / Rin_s
+  set dR = Rin_s - Rot_s
+
+  create n 0 dR[]
+  set dRf = dR / Rin_s
+  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
+  label -x sequence -y "1 - R_out| / R_in|"
+  resize 700 320
+
+  # check on A.Ratio
+  set ARot = rot_s / Rot_s
+  # lim ARin ARot; clear; box; plot ARin ARot
+  set fAR = ARot / ARin
+  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
+  label -x sequence -y "AR_out| / AR_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dM = Mot_s - Min_s
+  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM    
+  label -x sequence -y "M_out| - M_in|"
+  resize 700 320
+
+  lim -n 4$1 n -1.5 1.5; clear; box; plot n dX
+  label -x sequence -y "X_out| - X_in|"
+  resize 700 320
+
+  lim -n 5$1 n -1.5 1.5; clear; box; plot n dY
+  label -x sequence -y "Y_out| - Y_in|"
+  resize 700 320
+end
+
+macro grid.mkexp.sersic
+  if ($0 != 2)
+    echo "USAGE: grid.mkexp.devexp (dir)"
+    break
+  end
+
+  mkdir $1
+
+  $FakeConfig = -camera SIMTEST
+  $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
+  $FakeConfig = $FakeConfig -D PSASTRO:PSASTRO.CATDIR catdir.ref
+  $FakeConfig = $FakeConfig -Db STARS.FAKE F                         ; # only use stars from catdir.ref
+  $FakeConfig = $FakeConfig -Db MATCH.DENSITY F
+  $FakeConfig = $FakeConfig -Db PSF.CONVOLVE T
+  $FakeConfig = $FakeConfig -Db GALAXY.FAKE T                        ; # generate a "realistic" distribution of galaxies
+  $FakeConfig = $FakeConfig -Df GALAXY.MAG 17.0
+  $FakeConfig = $FakeConfig -Df GALAXY.GRID.MAG 14.5
+  $FakeConfig = $FakeConfig -Db GALAXY.GRID T                        ; # generate a grid of galaxies (constant mag)
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MIN 0 
+  $FakeConfig = $FakeConfig -Df GALAXY.THETA.MAX 180
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DX 300
+  $FakeConfig = $FakeConfig -Di GALAXY.GRID.DY 300
   $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_SERSIC
-
-  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.0
-  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.0
-  # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN 10.0
-  # $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX 10.0
-  # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MIN 0.25
-  # $FakeConfig = $FakeConfig -Df GALAXY.ARATIO.MAX 0.25
-  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MIN 1.66
-  # $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX 1.66
+  $BaseConfig = $FakeConfig
 
   $Nseq = 0
@@ -153 +407 @@
       foreach Aratio 0.25 0.5 1.0
         foreach fwhm 0.8 1.0 1.5
+          $FakeConfig = $BaseConfig
           $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MIN $Rmajor
           $FakeConfig = $FakeConfig -Df GALAXY.RMAJOR.MAX $Rmajor
@@ -160 +415 @@
           $FakeConfig = $FakeConfig -Df GALAXY.INDEX.MAX $index
           
-          sprint name "sersic.%02d" $Nseq
+          sprint name "$1/sersic.%02d" $Nseq
           mkexp $name $fwhm SERSIC
-          fitexp $name $name.fit SER\_CONV
+          echo "$Nseq : $index $Rmajor $Aratio $fwhm"
           $Nseq ++
         end
@@ -170 +425 @@
 end
 
-macro go.grid.check.devexp
-  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s
-  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s
+macro grid.fitexp.sersic.devexp
+  if ($0 != 3)
+    echo "USAGE: grid.fitexp.sersic.devexp (srcdir) (outdir)"
+    break
+  end
+
+  mkdir $2
 
   $Nseq = 0
-  foreach type EXP DEV
+  foreach index 1 2 3 4
     foreach Rmajor 3 10 30
       foreach Aratio 0.25 0.5 1.0
         foreach fwhm 0.8 1.0 1.5
-          sprint name "sample.%02d" $Nseq
-          ckgalaxy.load $name.dat $name.fit.cmf
+          sprint name "sersic.%02d" $Nseq
+          fitexp $1/$name $2/$name.fit EXP_CONV,DEV_CONV
           $Nseq ++
         end
@@ -188 +447 @@
 end
 
-macro go
-  mkexp test.exp 1.0 EXP
-  fitexp test.exp test.exp.fit EXP_CONV
-
-  mkexp test.ser 1.0 SERSIC
-  fitexp test.ser test.ser.fit SER_CONV
-
-  mkexp test.dev 1.0 DEV
-  fitexp test.dev test.dev.fit DEV_CONV
-
-  mkexp test.gau 1.0 GAUSS
-  fitexp test.gau test.gau.fit GAU_CONV
-
-  mkexp test.pg 1.0 PGAUSS
-  fitexp test.pg test.pg.fit PGA_CONV
-
-  mkexp test.qga 1.0 QGAUSS
-  fitexp test.qga test.qga.fit QGA_CONV
-
-  mkexp test.p1 1.0 PS1_V1
-  fitexp test.p1 test.p1.fit PS1_CONV
-end
-
-macro go.ckgalaxy
-  foreach type exp ser dev gau pg qga p1
-    ckgalaxy test.$type.dat test.$type.fit.cmf
-    wait $type
-  end
+macro grid.fitexp.sersic
+  if ($0 != 3)
+    echo "USAGE: grid.fitexp.sersic (srcdir) (outdir)"
+    break
+  end
+
+  mkdir $2
+
+  $Nseq = 0
+  foreach index 1 2 3 4
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          sprint name "sersic.%02d" $Nseq
+          fitexp $1/$name $2/$name.fit SER\_CONV
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+macro grid.load.sersic
+  if ($0 != 3)
+    echo "USAGE: grid.load.devexp (srcdir) (fitdir)"
+    break
+  end
+
+  delete -q Xin_s Yin_s Min_s Tin_s Rin_s rin_s MTin_s Iin_s
+  delete -q Xot_s Yot_s Mot_s Tot_s Rot_s rot_s MTot_s Iot_s
+
+  $Nseq = 0
+  foreach index 1 2 3 4
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          sprint name "sersic.%02d" $Nseq
+          cmf.load.concat $1/$name.dat $2/$name.fit.cmf SERSIC
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+# I want to make plots of Iin_s vs Mkron, Mxx, and similar things
+# this means I need to be able to join Chip.xfit things against Chip.psf
+# and to join Chip.xfit(DEV) to Chip.xfit(EXP)
+
+# I think I need a generic 'JOIN' function
+
+macro grid.plots.sersic
+  if ($0 != 2)
+    echo "USAGE: grid.plot.devexp (version)"
+    break
+  end
+  # things to examine: theta, Rmajor, AR
+
+  # go.grid.check.devexp
+
+  # check on theta
+  set dT = Tot_s - Tin_s
+  set ARin = rin_s / Rin_s
+  # lim ARin dT; clear; box; plot ARin dT
+
+  subset dTx = dT if (ARin < 0.95)
+  histogram dTx NdT -8 8 0.1 -range dTi
+  lim -n 0$1 dTi NdT; clear; box; plot -x 1 dTi NdT
+  label -x "angle offset (degrees, only non-circular)" -y "number count"
+  resize 700 320
+
+  # check on Rmajor
+  set dR = Rin_s - Rot_s
+  # lim Rin_s dR; clear; box; plot Rin_s dR
+  # lim Rot_s dR; clear; box; plot Rot_s dR
+  # lim Rin_s dR; clear; box; plot Rin_s dR
+  create n 0 dR[]
+  set dRf = dR / Rin_s
+  lim -n 1$1 n -0.5 0.5; clear; box; plot n dRf
+  label -x sequence -y "1 - R_out| / R_in|"
+  resize 700 320
+
+  # check on A.Ratio
+  set ARot = rot_s / Rot_s
+  # lim ARin ARot; clear; box; plot ARin ARot
+  set fAR = ARot / ARin
+  lim -n 2$1 n 0.5 1.5; clear; box; plot n fAR
+  label -x sequence -y "AR_out| / AR_in|"
+  resize 700 320
+
+  # check on magnitude
+  set dM = Mot_s - Min_s
+  lim -n 3$1 n -0.5 0.5; clear; box; plot n dM    
+  label -x sequence -y "M_out| - M_in|"
+  resize 700 320
+
+  # check on index
+  set dI = Iot_s - Iin_s
+  lim -n 4$1 n -0.5 0.5; clear; box; plot n dI    
+  label -x sequence -y "I_out| - I_in|"
+  resize 700 320
 end
 
@@ -258 +590 @@
 end
 
+if (not($?NSIGMA_CONV)) set NSIGMA_CONV = 5.0
+
 # create a realistic distribution of fake stars, GAUSS PSF
 macro fitexp
@@ -278 +612 @@
   $psphotConfig = $psphotConfig -Db PSPHOT:SAVE.RESID T
   $psphotConfig = $psphotConfig -D  PSPHOT:EXTENDED_SOURCE_MODELS_SELECTION $fitModel
+  $psphotConfig = $psphotConfig -D  PSPHOT:PSF_MODEL PS_MODEL_$PSFMODEL
+
+  $psphotConfig = $psphotConfig -Db PSPHOT:PSF.RESIDUALS F
+  $psphotConfig = $psphotConfig -Db PSPHOT:POISSON.ERRORS.PHOT.LMM F
+  $psphotConfig = $psphotConfig -Db PSPHOT:EXTENDED_SOURCE_FITS_POISSON F
+  $psphotConfig = $psphotConfig -Di PSPHOT:EXT_FIT_ITER 15
+  $psphotConfig = $psphotConfig -Df PSPHOT:PSF_FIT_RADIUS_SCALE 3.75
+  $psphotConfig = $psphotConfig -Df PSPHOT:EXT_FIT_NSIGMA_CONV $NSIGMA_CONV
 
   # ppImage / psphot on the output
+  break -auto off
   echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname
-  exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname
-end
-
-macro ckchip
-  if ($0 != 5)
-    echo "USAGE: ckchip (raw) (out) (output) (zpt_off)"
-    break
-  end
-
-  load.cmf $1 Chip.psf raw
-  load.cmf $2 Chip.psf out
-
-  # images generated with convolution will not have the right output positions
-  set X_raw = int(X_PSF_raw) + 0.5
-  set Y_raw = int(Y_PSF_raw) + 0.5
-  set M_raw = PSF_INST_MAG_raw + $4
-  set K_out = -2.5*log(KRON_FLUX_out)
-  match2d X_PSF_out Y_PSF_out X_PSF_raw Y_PSF_raw 1.5 -index1 index1 -index2 index2
-
-  local i NX NY nx ny N
-
-  device -n compare
-  resize 1000 1000
-
-  # plot trends as a function of mag
-  $NX = 2
-  $NY = 5
-  $nx = 0
-  $ny = 0
-  $N = 0
-  clear -s
-  for i 0 $pairs:n
-    section a$nx\$ny {$nx/$NX} {$ny/$NY} {1/$NX} {1/$NY}
-    show.pair $i
-    $ny ++
-    if ($ny == $NY)
-      $ny = 0
-      $nx ++
-    end
-    if ($nx == $NX)
-      png -name $3.$N.png
-      clear -s
-      $nx = 0
-      $ny = 0
-      $N ++
-    end
-  end
-
-  # plot (input - output) vs mag
-end
-
-macro ckgalaxy.load
-  if ($0 != 3)
-    echo "USAGE: ckgalaxy (dat) (cmf)"
+  exec echo ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >& $outname.log
+  exec ppImage $ppImageConfig $psphotConfig -file $basename.fits $outname >>& $outname.log
+  break -auto on
+end
+
+macro cmf.load.reset
+  $fields = X Y M T R r MT I
+  foreach field $fields
+    foreach set in ot
+      delete -q $field\$set\_s
+    end
+  end
+
+  delete -q min_S Min_S
+end
+
+macro cmf.load.concat
+  if ($0 != 4)
+    echo "USAGE: cmf.load.concat (dat) (cmf) (inType)"
     break
   end
 
   data $1
-  read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9
-
+  read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
+
+  $TYPE_S = 83
+  $TYPE_D = 68
+  $TYPE_E = 69
+  if ("$3" == "SERSIC")
+    $InType = $TYPE_S
+  end
+  if ("$3" == "DEV")
+    $InType = $TYPE_D
+    set IndexIn_all = 0.125 + zero(Xin_all)
+  end
+  if ("$3" == "EXP")
+    $InType = $TYPE_E
+    set IndexIn_all = 0.5 + zero(Xin_all)
+  end
+
+  # select only the galaxies 
   subset Xin = Xin_all if (Type == 1)
   subset Yin = Yin_all if (Type == 1)
   subset Min = Min_all if (Type == 1)
+  subset Fin = Fin_all if (Type == 1)
+  set min = -2.5*log(Fin)
 
   subset Tin_rad = ThetaIn_all if (Type == 1)
@@ -351 +677 @@
   subset RminIn = RminIn_all if (Type == 1)
 
+  subset IndexIn = IndexIn_all if (Type == 1)
+
   data $2
-  read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
+
+  break -auto off
+  output -err /dev/null
+  read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07
+  $reread = not($STATUS)
+  output -err stderr
+  break -auto on
+  if ($reread)
+    read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE 
+    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*0.125 + (MODEL_TYPE:9 == 69)*0.5
+  end
+
   set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
   set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
@@ -364 +703 @@
   reindex Yin_m = Yin using index2
 
+  set MTin_m = $InType + zero(Xin_m)
+  reindex MTot_m = MODEL_TYPE:9 using index1
+
   reindex Mot_m = EXT_INST_MAG using index1
   reindex Tot_m = EXT_THETA using index1
@@ -376 +718 @@
   reindex rin_m = RminIn using index2
   
-  foreach field X Y M T R r 
+  reindex Iot_m = EXT_PAR_07 using index1
+  reindex Iin_m = IndexIn using index2
+  $fields = X Y M T R r MT I
+  
+  foreach field $fields
     foreach set in ot
       concat $field\$set\_m $field\$set\_s
     end
   end
-end
-
-macro stats.pair
+
+  concat min min_S
+  concat Min Min_S
+end
+
+macro cmf.load.stars.concat
   if ($0 != 3)
-    echo "USAGE: stats.pair (N) (output)"
-    break
-  end
-
-  list word -split $spairs:$1
-  if ($word:n != 8)
-    echo "invalid pair $1"
-    break
-  end
-
-  $Nr = $word:3
-
-  reindex v1 = $word:0 using index1
-  reindex v2 = $word:1 using index2
-  reindex rv = $word:2 using index$Nr
-
-  set delta = v1 - v2
-  subset d1 = delta if ($word:4 < rv) && (rv < $word:5) && (abs(delta) < $word:7)
-  subset d2 = delta if ($word:5 < rv) && (rv < $word:6) && (abs(delta) < $word:7)
-
-  vstats -q d1 -sigma-clip 3.0
-  $M1 = $MEAN
-  $S1 = $SIGMA
-  vstats -q d2 -sigma-clip 3.0
-  $M2 = $MEAN
-  $S2 = $SIGMA
-
-  output $2
-  fprintf "%-18s  %7.4f %7.4f  %7.4f %7.4f" $word:0  $M1 $S1  $M2 $S2
-  output stdout
-end
-
-macro show.pair
-  if ($0 != 2)
-    echo "USAGE: show.pair (N)"
-    break
-  end
-
-  list word -split $pairs:$1
-  if ($word:n != 7)
-    echo "invalid pair $1"
-    break
-  end
-
-  $Nr = $word:3
-
-  reindex v1 = $word:0 using index1
-  if ("$word:6" == "V") 
-    reindex v2 = $word:1 using index2
-  end
-  if ("$word:6" == "S") 
-    set v2 = $word:1 + zero(index1)
-  end
-  reindex rv = $word:2 using index$Nr
-
-  set delta = v1 - v2
-  if (("$word:4" == "def") || ("$word:5" == "def"))
-    lim rv delta; box; plot rv delta
-  else
-    lim rv $word:4 $word:5; box; plot rv delta
-  end
-  label -y '$word:0' -x '$word:2'
-end
-
-# This list is used to compare a pair of vectors (sans error) or a
-# vector and an expected (constant) value.  The last field defines a
-# vector or constant for the comparison.  It is assumed that the
-# vector sets have been loaded and matched with match2d to generate
-# index vectors 'index1' and index2'.  The macro 'show.pair' generates
-# a plot of the range vector vs (v1 - v2).  The indices for v1, v2 are
-# index1 and 2 respectively.  The index for the range vector is defined
-# by the integer following that vector.  the y-limits of the plot are
-# given by the last two numbers
-list pairs
-  X_PSF_out             X_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
-  Y_PSF_out             Y_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
-  X_PSF_out             X_PSF_raw             PSF_INST_MAG_raw 2 -1.01 1.01 V
-  Y_PSF_out             Y_PSF_raw             PSF_INST_MAG_raw 2 -1.01 1.01 V
-  X_PSF_SIG_out         X_PSF_SIG_raw         PSF_INST_MAG_raw 2 -1.01 1.01 V
-  Y_PSF_SIG_out         Y_PSF_SIG_raw         PSF_INST_MAG_raw 2 -1.01 1.01 V
-  #PSF_INST_MAG_out      PSF_INST_MAG_raw      PSF_INST_MAG_raw 2 -1.01 1.01 V
-  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
-  PSF_INST_MAG_SIG_out  PSF_INST_MAG_SIG_raw  PSF_INST_MAG_raw 2 -1.01 1.01 V
-  #PSF_INST_FLUX_out     PSF_INST_FLUX_raw     PSF_INST_MAG_raw 2   def  def V
-  #PSF_INST_FLUX_SIG_out PSF_INST_FLUX_SIG_raw PSF_INST_MAG_raw 2   def  def V
-  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
-  AP_MAG_RAW_out        M_raw                 PSF_INST_MAG_raw 2 -1.01 1.01 V
-  AP_MAG_RADIUS_out     0.0                   PSF_INST_MAG_raw 2 -0.01 20.1 S
-  SKY_out               0.0                   PSF_INST_MAG_raw 2   def  def S
-  SKY_SIGMA_out         0.0                   PSF_INST_MAG_raw 2   def  def S
-  PSF_CHISQ_out         1.0                   PSF_INST_MAG_raw 2   def  def S
-  CR_NSIGMA_out         0.0                   PSF_INST_MAG_raw 2   def  def S
-  EXT_NSIGMA_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
-  PSF_MAJOR_out         0.0                   PSF_INST_MAG_raw 2 -0.01 5.01 S
-  PSF_MINOR_out         0.0                   PSF_INST_MAG_raw 2 -0.01 5.01 S
-  PSF_THETA_out         0.0                   PSF_INST_MAG_raw 2 -1.61 1.61 S
-  PSF_QF_out            0.0                   PSF_INST_MAG_raw 2 -0.10 1.10 S
-  PSF_QF_PERFECT_out    0.0                   PSF_INST_MAG_raw 2 -0.10 1.10 S
-  PSF_NDOF_out          0.0                   PSF_INST_MAG_raw 2  def  def  S
-  PSF_NPIX_out          0.0                   PSF_INST_MAG_raw 2  def  def  S
-  MOMENTS_XX_out        0.0                   PSF_INST_MAG_raw 2 -0.01 3.01 S
-  MOMENTS_XY_out        0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
-  MOMENTS_YY_out        0.0                   PSF_INST_MAG_raw 2 -0.01 3.01 S
-  MOMENTS_M3C_out       0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
-  MOMENTS_M3S_out       0.0                   PSF_INST_MAG_raw 2 -3.01 3.01 S
-  MOMENTS_M4C_out       0.0                   PSF_INST_MAG_raw 2 -2.01 2.01 S
-  MOMENTS_M4S_out       0.0                   PSF_INST_MAG_raw 2 -2.01 2.01 S
-  MOMENTS_R1_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
-  MOMENTS_RH_out        0.0                   PSF_INST_MAG_raw 2 -5.01 5.01 S
-  K_out                 M_raw                 PSF_INST_MAG_raw 2  def  def  V
-  KRON_FLUX_ERR_out     0.0                   PSF_INST_MAG_raw 2  def  def  S
-  KRON_FLUX_INNER_out   0.0                   PSF_INST_MAG_raw 2  def  def  S
-  KRON_FLUX_OUTER_out   0.0                   PSF_INST_MAG_raw 2  def  def  S
-# CAL_PSF_MAG          CAL_PSF_MAG          none Mraw 2 -1.01 1.01 V
-# CAL_PSF_MAG_SIG      CAL_PSF_MAG_SIG      none Mraw 2 -1.01 1.01 V
-# RA_PSF               RA_PSF               none Mraw 2 -1.01 1.01 V
-# DEC_PSF              DEC_PSF              none Mraw 2 -1.01 1.01 V
-# PEAK_FLUX_AS_MAG     PEAK_FLUX_AS_MAG     none Mraw 2 -1.01 1.01 V
-# FLAGS                FLAGS                0.0  Mraw 2 -1.01 1.01 S
-# FLAGS2               FLAGS2               0.0  Mraw 2 -1.01 1.01 S
-end
-
-macro load.cmf
-  if ($0 != 4)
-   echo "load.cmf (filename) (ext) (label)"
-   break
+    echo "USAGE: cmf.load.concat (dat) (cmf)"
+    break
   end
 
   data $1
-
-  # create the list of fields to load
-  delete -q myFields
-  for i 0 $fields:n
-    if ($?myFields) 
-      $myFields = $myFields $fields:$i
-    else
-      $myFields = $fields:$i
-    end
-  end
-
-  read -fits $2 $myFields
-
-  # rename the loaded vectors appending the supplied lable
-  for i 0 $fields:n
-    set $fields:$i\_$3 = $fields:$i
-    delete $fields:$i
-  end
-end
-
-# this list defines the fields to be loaded from file
-list fields
-  X_PSF              
-  Y_PSF              
-  X_PSF_SIG          
-  Y_PSF_SIG          
-  PSF_INST_MAG       
-  PSF_INST_MAG_SIG   
-  PSF_INST_FLUX      
-  PSF_INST_FLUX_SIG  
-  AP_MAG             
-  AP_MAG_RAW         
-  AP_MAG_RADIUS      
-  SKY                
-  SKY_SIGMA          
-  PSF_CHISQ          
-  CR_NSIGMA          
-  EXT_NSIGMA         
-  PSF_MAJOR          
-  PSF_MINOR          
-  PSF_THETA          
-  PSF_QF             
-  PSF_QF_PERFECT     
-  PSF_NDOF           
-  PSF_NPIX           
-  MOMENTS_XX         
-  MOMENTS_XY         
-  MOMENTS_YY         
-  MOMENTS_M3C        
-  MOMENTS_M3S        
-  MOMENTS_M4C        
-  MOMENTS_M4S        
-  MOMENTS_R1         
-  MOMENTS_RH         
-  KRON_FLUX          
-  KRON_FLUX_ERR      
-  KRON_FLUX_INNER    
-  KRON_FLUX_OUTER    
-#   CAL_PSF_MAG      
-#   CAL_PSF_MAG_SIG  
-#   RA_PSF           
-#   DEC_PSF          
-#   PEAK_FLUX_AS_MAG 
-#  FLAGS            
-#  FLAGS2           
-end
-
-# use these cmf entries to measure average stats of the given pairs
-list spairs
-  X_PSF_out             X_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
-  Y_PSF_out             Y_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
-  PSF_INST_MAG_out      M_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
-  AP_MAG_out            M_raw                 PSF_INST_MAG_raw 2 -15 -10 -8.5  0.1
-end
-
-macro show.dpair
-  if ($0 != 2)
-    echo "USAGE: show.dpair (N)"
-    break
-  end
-
-  list word -split $dpairs:$1
-  if ($word:n != 7)
-    echo "invalid dpair $1"
-    break
-  end
-
-  $Nr = $word:4
-
-  reindex v1 = $word:0 using index1
-  reindex dv = $word:1 using index1
-  reindex v2 = $word:2 using index2
-  reindex rv = $word:3 using index$Nr
-
-  set delta = (v1 - v2) / dv
-  if (("$word:5" == "def") || ("$word:6" == "def"))
-    lim rv delta; box; plot rv delta
-  else
-    lim rv $word:5 $word:6; box; plot rv delta
-  end
-  label -y '$word:0' -x '$word:3'
-end
-
-# this list is used to compare a pair of vectors with an error it is
-# assumed that the vector sets have been loaded and matched with
-# match2d to generate index vectors 'index1' and index2'.  the macro
-# show.dpair generates a plot of the range vector vs (v1 - v2) / dv1.
-# The indices for v1, dv1, and v2 are index1,1, and 2 respectively.  The
-# index for the range vector is defined by the integer following that
-# vector.  The y-limits of the plot are given by the last two numbers
-# (use 'def') for the full default range of the delta vector
-list dpairs
-  # v1              dv                v2   range 
-  X_PSF             X_PSF_SIG         Xraw Mraw  2 -10.0 10.0
-  Y_PSF             Y_PSF_SIG         Yraw Mraw  2 -10.0 10.0
-  PSF_INST_MAG      PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
-  PSF_INST_FLUX     PSF_INST_FLUX_SIG Fraw Mraw  2 -10.0 10.0
-  AP_MAG            PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
-  AP_MAG_RAW        PSF_INST_MAG_SIG  Mraw Mraw  2 -10.0 10.0
+  read Xin_all 1 Yin_all 2 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
+
+  subset Xin     = Xin_all     if (Type == 0)
+  subset Yin     = Yin_all     if (Type == 0)
+  subset Min     = Min_all     if (Type == 0)
+  subset RmajIn  = RmajIn_all  if (Type == 0)
+  subset RminIn  = RminIn_all  if (Type == 0)
+  subset IndexIn = IndexIn_all if (Type == 0)
+  subset Tin_rad = ThetaIn_all if (Type == 0)
+  set Tin = Tin_rad * 180 / 3.14159265
+
+  data $2
+
+  break -auto off
+  read -fits Chip.psf X_PSF Y_PSF PSF_INST_MAG PSF_MAJOR PSF_MINOR PSF_THETA
+  set PSF_THETA_ALT = PSF_THETA * (PSF_THETA >= 0.0) + (PSF_THETA + 3.14159265) * (PSF_THETA < 0.0)
+  set PSF_THETA = PSF_THETA_ALT * 180 / 3.14159265
+  
+  match2d X_PSF Y_PSF Xin Yin 1.0 -index1 index1 -index2 index2
+
+  reindex Xot_m = X_PSF using index1
+  reindex Yot_m = Y_PSF using index1
+
+  reindex Xin_m = Xin using index2
+  reindex Yin_m = Yin using index2
+
+  reindex Mot_m = PSF_INST_MAG using index1
+  reindex Tot_m = PSF_THETA using index1
+
+  reindex Min_m = Min using index2
+  reindex Tin_m = Tin using index2
+
+  reindex Rot_m = PSF_MAJOR using index1
+  reindex rot_m = PSF_MINOR using index1
+
+  reindex Rin_m = RmajIn using index2
+  reindex rin_m = RminIn using index2
+  
+  $fields = X Y M T R r
+  foreach field $fields
+    foreach set in ot
+      concat $field\$set\_m $field\$set\_s
+    end
+  end
+end
+
+macro grid.load.sersic.test
+  if ($0 != 3)
+    echo "USAGE: grid.load.devexp.test (srcdir) (fitdir)"
+    break
+  end
+
+  $fields_bt = X Y M T R r I
+  $fields_ot = Pmag Kmag Amag IDx MTot
+
+  foreach field $fields_bt
+    foreach set in ot
+      delete -q $field\$set\_dev_s
+      delete -q $field\$set\_exp_s
+    end
+  end
+  foreach field $fields_ot
+    delete -q $field\_dev_s
+    delete -q $field\_exp_s
+  end
+
+  $Nseq = 0
+  foreach index 1 2 3 4
+    foreach Rmajor 3 10 30
+      foreach Aratio 0.25 0.5 1.0
+        foreach fwhm 0.8 1.0 1.5
+          sprint name "sersic.%02d" $Nseq
+          cmf.load.sersic.test $1/$name.dat $2/$name.fit.cmf
+          $Nseq ++
+        end
+      end
+    end
+  end
+end
+
+# I have run DEV and EXP against input models of type SERSIC
+macro cmf.load.sersic.test
+  if ($0 != 3)
+    echo "USAGE: cmf.load.sersic.test (dat) (cmf)"
+    break
+  end
+
+  # input parameters
+  data $1
+  read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
+
+  # galaxies only
+  subset Xin = Xin_all if (Type == 1)
+  subset Yin = Yin_all if (Type == 1)
+  subset Min = Min_all if (Type == 1)
+  subset Fin = Fin_all if (Type == 1)
+  set min = -2.5*log(Fin)
+
+  subset Tin_rad = ThetaIn_all if (Type == 1)
+  set Tin = Tin_rad * 180 / 3.14159265
+
+  subset RmajIn = RmajIn_all if (Type == 1)
+  subset RminIn = RminIn_all if (Type == 1)
+
+  subset IndexIn = IndexIn_all if (Type == 1)
+
+  $TYPE_S = 83
+  $TYPE_D = 68
+  $TYPE_E = 69
+
+  data $2
+
+  # load measured values from xfit
+  read -fits Chip.xfit IPP_IDET X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE PSF_INST_MAG AP_MAG KRON_MAG
+  set EXT_PAR_07 = (MODEL_TYPE:9 == $TYPE_D)*4 + (MODEL_TYPE:9 == $TYPE_E)
+
+  set EXT_THETA_ALT = EXT_THETA * (EXT_THETA >= 0.0) + (EXT_THETA + 3.14159265) * (EXT_THETA < 0.0)
+  set EXT_THETA = EXT_THETA_ALT * 180 / 3.14159265
+  set IPP_IDET_EXT = IPP_IDET
+  
+  match2d X_EXT Y_EXT Xin Yin 1.0 -index1 index1 -index2 index2
+
+  reindex Xot_m = X_EXT using index1
+  reindex Yot_m = Y_EXT using index1
+
+  reindex Xin_m = Xin using index2
+  reindex Yin_m = Yin using index2
+
+  reindex Mot_m = EXT_INST_MAG using index1
+  reindex Tot_m = EXT_THETA using index1
+
+  reindex Min_m = Min using index2
+  reindex Tin_m = Tin using index2
+
+  reindex Rot_m = EXT_WIDTH_MAJ using index1
+  reindex rot_m = EXT_WIDTH_MIN using index1
+
+  reindex Rin_m = RmajIn using index2
+  reindex rin_m = RminIn using index2
+  
+  reindex Pmag_m = PSF_INST_MAG using index1
+  reindex Kmag_m = KRON_MAG using index1
+  reindex Amag_m = AP_MAG using index1
+
+  reindex IDx_m = IPP_IDET_EXT using index1
+
+  reindex MTot_m = MODEL_TYPE:9 using index1
+
+  reindex Iot_m = EXT_PAR_07 using index1
+  reindex Iin_m = IndexIn using index2
+  
+  # load moments and other kron values from Chip.psf
+  read -fits Chip.psf IPP_IDET X_PSF Y_PSF MOMENTS_XX MOMENTS_XY MOMENTS_YY KRON_FLUX_INNER MOMENTS_R1 MOMENTS_RH
+  set IPP_IDET_PSF = IPP_IDET
+
+  join -outer IPP_IDET_PSF IDx_m 
+  reindex Xp = X_PSF using index2
+  reindex Yp = Y_PSF using index2
+
+  reindex Mxx_m = MOMENTS_XX using index2
+  reindex Mxy_m = MOMENTS_XY using index2
+  reindex Myy_m = MOMENTS_YY using index2
+  reindex Mr1_m = MOMENTS_R1 using index2
+  reindex Mrh_m = MOMENTS_RH using index2
+  reindex Kfi_m = KRON_FLUX_INNER using index2
+  set Kmi_m = -2.5*log(Kfi_m)
+
+  $fields_bt = X Y M T R r I
+  $fields_ot = Pmag Kmag Amag IDx MTot Mxx Mxy Myy Mr1 Mrh Kmi
+
+  foreach field $fields_ot
+    subset $field\_exp_m = $field\_m where (MTot_m == $TYPE_E)
+    subset $field\_dev_m = $field\_m where (MTot_m == $TYPE_D)
+  end
+  foreach field $fields_bt
+    foreach set in ot
+      subset $field\$set\_exp_m = $field\$set\_m where (MTot_m == $TYPE_E)
+      subset $field\$set\_dev_m = $field\$set\_m where (MTot_m == $TYPE_D)
+    end
+  end
+
+  join IDx_exp_m IDx_dev_m
+  foreach field $fields_ot
+    reindex $field\_exp_mr = $field\_exp_m using index1
+    reindex $field\_dev_mr = $field\_dev_m using index2
+  end
+  foreach field $fields_bt
+    foreach set in ot
+      reindex $field\$set\_exp_mr = $field\$set\_exp_m using index1
+      reindex $field\$set\_dev_mr = $field\$set\_dev_m using index2
+    end
+  end
+
+  # concat
+  foreach field $fields_bt
+    foreach set in ot
+      concat $field\$set\_dev_mr $field\$set\_dev_s
+      concat $field\$set\_exp_mr $field\$set\_exp_s
+    end
+  end
+  foreach field $fields_ot
+    concat $field\_dev_mr $field\_dev_s
+    concat $field\_exp_mr $field\_exp_s
+  end
+
+  concat min min_S
+  concat Min Min_S
 end
 
@@ -659 +972 @@
 end
 
-macro plot.angles
+macro sersic.integral
   if ($0 != 2)
-    echo "USAGE: plot.angles (file.dat)"
+    echo "sersic.integral (index)"
+    break
+  end
+
+  $index = $1
+  create r 0.0 200.0 0.001
+  set q = r^(1/$index)
+  set f = exp(-q)
+  set rf = f*r
+  integrate r rf r[0] r[-1]
+end
+
+macro sersic.integral.rmax
+  if ($0 != 4)
+    echo "sersic.integral.rmax (index) (kappa) (rmax)"
+    echo "kappa is a guess for kappa"
+    break
+  end
+
+  local index kappa
+
+  $index = $1
+  $kappa = $2
+
+  create r 0.0 $3 0.003
+  set q = r^(1/$index)
+  set f = exp(-$kappa*q)
+  set rf = f*r
+  integrate r rf r[0] r[-1]
+end
+
+# S(r) = exp(-kappa*(r/Reff)^(1/index))
+# integrate S(r) r dr [ignores 2pi and change-of-variable factors (Rmaj*Rmin)]
+macro sersic.integral.reff.rmax
+  if ($0 != 5)
+    echo "sersic.integral.rmax (index) (kappa) (reff) (rmax)"
+    echo "kappa is a guess for kappa"
+    break
+  end
+
+  local index kappa Reff
+
+  $index = $1
+  $kappa = $2
+  $Reff = $3
+
+  create r 0.0 $4 0.01
+  set q = (r/$Reff)^(1/$index)
+  set f = exp(-$kappa*q)
+  set rf = f*r
+  integrate r rf r[0] r[-1]
+end
+
+# integrate to r = reff (rho = 1), applying kappa
+macro sersic.integral.reff
+  if ($0 != 3)
+    echo "sersic.integral.reff (index) (kappa)"
+    break
+  end
+
+  local index kappa
+
+  $index = $1
+  $kappa = $2
+  create r 0.0 100.0 0.001
+  set q = r^(1/$index)
+  set f = exp(-$kappa*q)
+  set rf = f*r
+  integrate r rf r[0] r[-1]
+  $S0 = $sum
+  integrate r rf r[0] 1.0
+  $S1 = $sum
+  echo $S0 $S1 {$S1 / $S0}
+end
+
+macro sersic.integral.find.reff
+  if ($0 != 3)
+    echo "sersic.integral.find.reff (index) (kappa)"
+    echo "kappa is a guess for kappa"
+    break
+  end
+
+  local kappa index Rmax i
+  $index = $1
+  $kappa = $1
+
+  $Rmin = 0; $Vmin = 0
+  $Rmax = 4000
+  sersic.integral.rmax $index $kappa $Rmax
+  $Vmax = $sum
+
+  $Rtry = 0.5*($Rmin + $Rmax)
+
+  while (abs($Rmax - $Rmin) > 0.001)
+    sersic.integral.rmax $index $kappa $Rtry
+    $Vtry = $sum
+
+    $Rold = $Rtry
+    if ($Vtry > 0.5*$Vmax)
+      $Rtry = 0.5*($Rmin + $Rtry)
+      $Rmax = $Rold
+    else
+      $Rtry = 0.5*($Rmax + $Rtry)
+      $Rmin = $Rold
+    end
+    # fprintf "%5.2f %5.2f %5.2f | %5.2f | %6.3f %6.3f | %6.3f" $Rmin $Rold $Rmax $Rtry {$Vtry / $Vmax} $Vmax $Vtry
+  end
+  # echo $Rtry
+  # echo {$kappa * $Rtry ^ (1.0/$index)}
+  $KappaReal = $kappa * $Rtry ^ (1.0/$index)
+end
+
+macro sersic.integral.index
+
+  local i kappa
+
+  delete -q sersic_sum sersic_norm 
+  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
+  # create idx 0.5 5.1 0.1
+  vlist idx 1 2 4
+  for i 0 idx[]
+    $kappa =  -0.275552 + 1.972625*idx[$i] + 0.003487*idx[$i]*idx[$i]
+    sersic.integral.rmax idx[$i] $kappa 300
+    concat $sum sersic_sum
+
+    $bn = 1.9992*idx[$i] - 0.3271;
+    $Io = exp($bn);
+            
+    # the integral of a Sersic (supposedly) has an analytical form as follows:
+    $logGamma = lgamma(2.0*idx[$i]);
+    $bnFactor = $bn^(2.0*idx[$i]);
+    $norm = idx[$i] * $Io * exp($logGamma) / $bnFactor;
+    concat $norm sersic_norm
+
+    echo idx[$i] $kappa $sum $logGamma $bn $bnFactor $norm
+  end
+  set lsersic_sum = ln(sersic_sum)
+  break
+  
+  $order = 2
+  delete -q sersic_fit fit_idx
+
+  # fit to ranges:
+  delete lsersic_fit fit_idx
+  vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
+  for i 0 {bound[] - 1}
+    subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1]) 
+    subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1]) 
+    fit idxs sums $order
+    applyfit idxs fits
+    concat fits lsersic_fit
+    concat idxs fit_idx
+  end
+  set sersic_fit = exp(lsersic_fit)
+
+  lim -n 1 idx sersic_sum; clear; box; plot idx sersic_sum; plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit
+  lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
+end
+
+macro sersic.kappa.index
+
+  local i
+
+  delete -q sersic_kappa
+  # vlist idx 0.5 0.75 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0
+  create idx 0.5 5.1 0.1
+  for i 0 idx[]
+    sersic.integral.find.reff idx[$i] 10
+    concat $KappaReal sersic_kappa
+  end
+  
+  # $order = 2
+  # delete -q sersic_kfit kfit_idx
+  # 
+  # # fit to ranges:
+  # vlist bound 0.0 1.0 2.0 3.0 4.0 5.5
+  # for i 0 {bound[] - 1}
+  #   subset idxs = idx if (idx >= bound[$i]) && (idx < bound[$i+1]) 
+  #   subset sums = lsersic_sum if (idx >= bound[$i]) && (idx < bound[$i+1]) 
+  #   fit idxs sums $order
+  #   applyfit idxs fits
+  #   concat fits lsersic_fit
+  #   concat idxs fit_idx
+  # end
+  # set sersic_fit = exp(lsersic_fit)
+
+  lim -n 1 idx sersic_kappa; clear; box; plot idx sersic_kappa; # plot -c blue -pt 7 idx sersic_norm; plot -x 0 -c red fit_idx sersic_fit
+  # lim -n 2 idx lsersic_sum; clear; box; plot idx lsersic_sum; plot -x 0 -c red fit_idx lsersic_fit
+
+  # kappa(index) :
+  # y = -0.275552 x^0 1.972625 x^1 0.003487 x^2 
+end
+
+macro sersic.central.pixel
+  if ($0 != 3)
+    echo "USAGE: sersic.central.pixel (index) (Reff)"
+    break
+  end
+
+  local index Reff
+
+  $index = $1
+  $Reff = $2
+
+  $kappa = -0.275552 + 1.972625*$index + 0.003487 * $index^2
+  # echo "kappa: $kappa"
+
+  sersic.norm $index
+  sersic.integral.reff.rmax $index $kappa $Reff {10*$Reff}; set sumFull = $sum; echo $sum {$Reff^2*$myNorm} {$sum / ($Reff^2*$myNorm)}
+  # sersic.integral.reff.rmax $index $kappa $Reff $Reff     ; # echo $sum {$sum / $sumFull}
+  # sersic.integral.reff.rmax $index $kappa $Reff 0.1       ; # echo $sum {$sum / $sumFull}
+  # sersic.integral.reff.rmax $index $kappa $Reff 0.2       ; # echo $sum {$sum / $sumFull}
+  # sersic.integral.reff.rmax $index $kappa $Reff 0.5       ; # echo $sum {$sum / $sumFull}
+  # sersic.integral.reff.rmax $index $kappa $Reff 1.0       ; # echo $sum {$sum / $sumFull}
+  # sersic.integral.reff.rmax $index $kappa $Reff 2.0       ; # echo $sum {$sum / $sumFull}
+  sersic.integral.reff.rmax $index $kappa $Reff 0.564     ; echo $index $Reff {$sum / $sumFull}
+
+end
+
+macro sersic.norm
+  if ($0 != 2)
+    echo "USAGE: sersic.norm (index)"
+    break
+  end
+
+  local index
+  $index = $1
+
+  if (($index >= 0.0) && ($index < 1.0))
+      $norm = 0.201545  - 0.950965 * $index - 0.315248 * $index^2
+      echo $norm {exp($norm)}
+      $myNorm = exp($norm)
+      return
+  end
+
+  if (($index >= 1.0) && ($index < 2.0)) 
+      $norm = 0.402084  - 1.357775 * $index - 0.105102 * $index^2
+      echo $norm {exp($norm)}
+      $myNorm = exp($norm)
+      return
+  end
+
+  if (($index >= 2.0) && ($index < 3.0))
+      $norm = 0.619093 - 1.591674 * $index - 0.041576 * $index^2
+      echo $norm {exp($norm)}
+      $myNorm = exp($norm)
+      return
+  end
+
+  if (($index >= 3.0) && ($index < 4.0))
+      $norm = 0.770263 - 1.696421 * $index - 0.023363 * $index^2
+      echo $norm {exp($norm)}
+      $myNorm = exp($norm)
+      return
+  end
+
+  if (($index >= 4.0) && ($index < 5.5)) 
+      $norm = 0.885891 - 1.755684 * $index - 0.015753 * $index^2
+      echo $norm {exp($norm)}
+      $myNorm = exp($norm)
+      return
+  end
+end
+
+macro load.model.im
+  if ($0 != 2)
+    echo "USAGE: load.model.im (N)"
+    break
+  end
+
+  rd obj$1 obj.$1.fits
+  rd cnv$1 cnv.$1.fits
+  rd var$1 var.$1.fits
+  rd msk$1 msk.$1.fits
+  for i 1 7
+    rd dpar$i.$1 dpar.$i.$1.fits
+  end
+  set dC$1 = (msk$1 == 0) * (obj$1 - cnv$1)^2 / var$1
+  tv dC$1 -0.01 3.0
+end
+
+macro load.normdata
+  if ($0 != 3)
+    echo "USAGE: load.normdata (file) [clear/noclear]"
     break
   end
 
   data $1
-  read x 1 y 2 type 4 trad 9 
-  set t = trad * 180 / 3.14159265
-  subset xs = x if (type == 1)
-  subset ys = y if (type == 1)
-  subset ts = t if (type == 1)
-
-  set cs = dcos(ts)
-  set sn = dsin(ts)
-
-  delete xp yp
-  erase red
-  for i 0 xs[]
-    point red LINE xs[$i] ys[$i] {10*cs[$i]} {10*sn[$i]}
-  end
-  plot xp yp -pt 100
-end
-
-macro check.flux
-  if ($0 != 2)
-    echo "USAGE: plot.angles (file.dat)"
-    break
-  end
-
-  data $1
-  read f 3 type 4
-  subset F = f if (type == 1)
-  set M = -2.5*log(F)
-  vstat M
+  read f 3 t 4 m 5
+  set M = -2.5*log(f)
+  subset Mg = M if (t == 1)
+  subset mg = m if (t == 1)
+  set dm = mg - Mg
+  set n = ramp(dm)
+  if ("$2" == "clear")
+    lim n dm; clear; box; 
+  end
+  plot n dm
 end
 
@@ -700 +1278 @@
   exit 0
 end
+
+# note that t1 = original confi
+# t2 = no residuals
+# t3 = const weight PSF fit
+# t4 = const weight galaxy fits
+
+macro run.radius.loop
+  local radius Nrun 
+
+  if (1)
+    $Nrun = 0
+    foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+      mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
+      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4b EXP_CONV
+      $Nrun ++
+    end
+  end
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end
+
+macro check.radius.loop
+  if ($0 != 2)
+    echo "USAGE: check.radius.loop (ext)"
+    break
+  end
+
+  # cmf.load.reset
+
+  local radius Nrun 
+
+  $Nrun = 0
+  foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
+    $Nrun ++
+  end
+
+  grid.plots.devexp a
+end
+
+macro check.fit
+  if ($0 != 4)
+    echo "USAGE: check.fit (dat) (cmf) (type)"
+    break
+  end
+
+  cmf.load.concat $1 $2 $3
+  set dM_m = Mot_m - Min_m
+  vstat -q dM_m
+  $mag_off = $MEDIAN
+  set dR_m = Rot_m - Rin_m
+  vstat -q dR_m
+  $rad_off = $MEDIAN
+  echo $mag_off $rad_off
+end
+
+macro make.circles
+
+  delete Xo Yo Ro
+
+  for ix 0 16
+    for iy $ix 16
+      $r2 = $ix^2 + $iy^2
+      if ($r2 > 15^2) continue
+      concat $ix Xo
+      concat $iy Yo
+      concat $r2 Ro
+    end
+  end
+
+  sort Ro Xo Yo
+  $Rold = -1
+  $N = -1
+  for i 0 Ro[]
+    if (Ro[$i] != $Rold)
+      $N ++
+      $Rold = Ro[$i]
+    end
+    if ((Xo[$i] == 0) && (Yo[$i] == 0))
+      fprintf "// center is 0,0"
+      continue
+    end
+    if (Xo[$i] == 0)
+      fprintf "ADD_AXIS (%3d, %2d)     // r^2 = %3d" $N Yo[$i] Ro[$i]
+      continue
+    end
+    if (Xo[$i] == Yo[$i])
+      fprintf "ADD_DIAG (%3d, %2d)     // r^2 = %3d" $N Xo[$i] Ro[$i]
+      continue
+    end
+
+    fprintf "ADD_RAND (%3d, %2d, %2d) // r^2 = %3d" $N Xo[$i] Yo[$i] Ro[$i]
+ end
+end
+
+# EXP : tests.20131120/test.nsig
+# DEV : tests.20131120/test.dev
+
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.nsig   EXP EXP_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.dev    DEV DEV_CONV
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.serexp EXP SER_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.serdev DEV SER_CONV
+
+macro run.convolve.loop
+  if ($0 != 5)
+    echo "run.convolve.loop (inName) (fitName) (inType) (fitType)"
+    break
+  end
+
+  local radius
+  $radius = 3.0
+  cmf.load.reset
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    $CONVOLVE_NSIGMA = $Cin
+    sprintf nameIn %s.%02d $1 $Cin
+    mkexp.devexp.single $nameIn $3 $radius 1.0
+    foreach Cot 3 5 7 9 11
+      $NSIGMA_CONV = $Cot
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      fitexp $nameIn $nameOt $4
+    end
+  end
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    sprintf nameIn %s.%02d $1 $Cin
+    foreach Cot 3 5 7 9 11
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      check.fit $nameIn.dat $nameOt.cmf 
+    end
+  end
+  grid.plots.devexp a
+end
+
+macro check.convolve.loop
+  if ($0 != 4)
+    echo "run.convolve.loop (inName) (fitName) (inType)"
+    break
+  end
+
+  local radius
+  $radius = 3.0
+  cmf.load.reset
+
+  # foreach Cin 3 5 7 9 11
+  foreach Cin 11
+    sprintf nameIn %s.%02d $1 $Cin
+    foreach Cot 3 5 7 9 11
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      echo -no-return $Cin $Cot " "
+      check.fit $nameIn.dat $nameOt.cmf $3
+    end
+  end
+  grid.plots.devexp a
+end

@@ -24 +24 @@
 psphotModelTest
 psphotMinimal
+psphotFullForce
+psmakecff
+psphotFullForceSummary