Changeset 36086 for trunk/psphot

Timestamp:

Aug 31, 2013, 6:02:26 AM (13 years ago)

Author:

eugene

Message:

merge changes from EAM branch (update tap script for galaxies; count deblended satstars for log message; scale down EXT window a bit for speed (OK?); enable interactive mode for PCM fitting; allow NONCONVERGE extended fit models to succeed; init EXT model fits with EXT_AND_SKY just in case; extra verbosity and various experiments on the index iterations)

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotBlendFit.c (modified) (3 diffs)
• src/psphotDeblendSatstars.c (modified) (3 diffs)
• src/psphotExtendedSourceFits.c (modified) (10 diffs)
• src/psphotRadiusChecks.c (modified) (1 diff)
• src/psphotSourceFits.c (modified) (6 diffs)
• src/psphotStackImageLoop.c (modified) (1 prop)
• test/tap_psphot_galaxygrid.pro (modified) (18 diffs)

trunk/psphot/src/psphotBlendFit.c

@@ -273 +273 @@
         pmSource *source = sources->data[i];
 
+# if (0)
         int TEST_ON = false;
-# if (0)
 # define TEST_X 653
 # define TEST_Y 466
@@ -342 +342 @@
         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 +353 @@
         } 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);

trunk/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;
 }

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -65 +65 @@
     psphotVisualShowImage(readout);
 
+    // psphotSaveImage (NULL, readout->image, "test.01.fits");
+
     pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     psAssert (detections, "missing detections?");
@@ -123 +125 @@
     }
 
+    // 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
@@ -134 +140 @@
     fitOptions->gainFactorMode   = gainFactorMode;
     fitOptions->chisqConvergence = chisqConvergence;
+    fitOptions->isInteractive    = isInteractive;
 
     // maskVal is used to test for rejected pixels, and must include markVal
@@ -296 +303 @@
             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
         }
 
@@ -378 +386 @@
     psphotSersicModelClassCleanup();
 
+    // psphotSaveImage (NULL, readout->image, "test.02.fits");
+
     psphotVisualShowResidualImage (readout, false);
 
@@ -521 +531 @@
                   psFree (source->modelFlux);
                   source->modelFlux = NULL;
+                  // if ((fabs(source->peak->xf - 2572) < 20) && (fabs(source->peak->yf - 5874) < 20)) {
+                  //   fprintf (stderr, "test object\n");
+                  // }
+
                   modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
                   if (!modelFit) {
@@ -544 +558 @@
                       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, 
@@ -579 +591 @@
 
             if (model->flags & (PM_MODEL_STATUS_BADARGS)) continue;
-            if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
+            // if (model->flags & (PM_MODEL_STATUS_NONCONVERGE)) continue;
             if (model->flags & (PM_MODEL_STATUS_OFFIMAGE)) continue;
 
@@ -598 +610 @@
           pmSourceCacheModel (source, maskVal);
 
+          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);
           pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
@@ -615 +634 @@
         // cache the model flux
         if (source->modelEXT->isPCM) {
+            // 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);
         } 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);
         }

trunk/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
 

trunk/psphot/src/psphotSourceFits.c

@@ -559 +559 @@
 
 # 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) {
@@ -580 +586 @@
 
     float t1, t2, t4, t5;
+    t1 = t2 = t4 = t5 = 0.0;
     if (TIMING) { psTimerStart ("psphotFitPCM"); }
+
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        options.mode = PM_SOURCE_FIT_NO_INDEX; // XXX note that there may be a conflict with psphotExtendedSourceFits.c:133
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+        // 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
+    } else {
+        options.mode = PM_SOURCE_FIT_EXT_AND_SKY;
+    }
 
     pmPCMdata *pcm = pmPCMinit (source, &options, model, maskVal, psfSize);
@@ -594 +611 @@
     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;
+            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_STATUS_BADARGS;
+            psError(PS_ERR_UNKNOWN, true, "Failed to find a index & shape");
+            psErrorClear (); // clear the polynomial error
             return model;
         }
@@ -608 +636 @@
     }
 
-    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"); }
 
@@ -629 +648 @@
         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);
+    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 +771 @@
     return true;
 }
+
+// float indexGuessInv[] = {0.5, 0.33, 0.25, 0.167, 0.125, 0.083};
+float indexGuessInv[] = {0.5, 0.4, 0.3, 0.25, 0.20, 0.15, 0.125};
+# define N_INDEX_GUESS_INV 7
+
+// float reffGuess[] = {3.0, 10.0, 20.0, 30.0, 40.0};
+float reffGuess[] = {0.5, 0.75, 1.0, 1.4, 2.0};
+# 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
+    // * 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);
+
+    if (0) {
+        psEllipseAxes guessAxes;
+        guessAxes.major = source->moments->Mrf;
+        guessAxes.minor = (momentAxes.minor / momentAxes.minor) * guessAxes.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);
+    }
+
+    // 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++) {
+
+        psEllipseAxes guessAxes;
+        guessAxes.major = reffGuess[j] * 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);
+
+        for (int i = 0; i < N_INDEX_GUESS_INV; i++) {
+            PAR[PM_PAR_7] = indexGuessInv[i];
+
+            // generated the modelFlux
+            pmPCMMakeModel (source, pcm->modelConv, 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];
+            }
+            // 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);
+
+    PAR[PM_PAR_7] = indexGuessInv[nStart];
+    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] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : 0.1;
+    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] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : 0.55;
+    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[2]) {
+      if (nStart == N_INDEX_GUESS_INV - 1) {
+        PAR[PM_PAR_7] = 0.11;
+      } else {
+        PAR[PM_PAR_7] = indexGuessInv[nStart + 1];
+      }
+    } else {
+      if (nStart == 0) {
+        PAR[PM_PAR_7] = 0.52;
+      } else {
+        PAR[PM_PAR_7] = indexGuessInv[nStart - 1];
+      }
+    }
+
+    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 (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;
+
+    // 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;
+}
+
+ 

trunk/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
 
 # options for the simulated images (using the refimage for the stars)
@@ -33 +34 @@
 $FakeOptions = $FakeOptions -exptime 30.0
 $FakeOptions = $FakeOptions -D PSF.MODEL PS_MODEL_GAUSS
+$FakeOptions = $FakeOptions -nx 3000 -ny 3000
   
 # sample alternate options:
@@ -65 +67 @@
 $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
 
 list fwhm 
@@ -75 +77 @@
 end
 
-macro go.grid.set.devexp
+# generate fake images and run psphot on them
+macro normtest.mkexp.devexp
   $FakeConfig = -camera SIMTEST
   $FakeConfig = $FakeConfig -recipe PPSIM STACKTEST.RUN
@@ -84 +87 @@
   $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
 
   # $FakeConfig = $FakeConfig -D GALAXY.MODEL PS_MODEL_GAUSS
@@ -103 +108 @@
   # $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 +121 @@
           $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 +130 @@
 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 +147 @@
   $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
+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 +304 @@
       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 +312 @@
           $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 +322 @@
 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 +344 @@
 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
 
@@ -280 +509 @@
 
   # 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.concat
+  if ($0 != 4)
+    echo "USAGE: cmf.load.concat (dat) (cmf) (type)"
     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
 
   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 +537 @@
   subset RminIn = RminIn_all if (Type == 1)
 
+  subset IndexIn = IndexIn_all if (Type == 1)
+
+  $TYPE_S = 83
+  $TYPE_D = 68
+  $TYPE_E = 69
+  if ("$3" == "SERSIC")
+    $InType = $TYPE_S
+  end
+  if ("$3" == "DEV")
+    $InType = $TYPE_D
+  end
+  if ("$3" == "EXP")
+    $InType = $TYPE_E
+  end
+
   data $2
-  read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA
+
+  break -auto off
+  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)
+  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)*4 + (MODEL_TYPE:9 == 69)
+  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 +574 @@
   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 +589 @@
   reindex rin_m = RminIn using index2
   
-  foreach field X Y M T R r 
+  if ("$3" == "SERSIC")
+    reindex Iot_m = EXT_PAR_07 using index1
+    reindex Iin_m = IndexIn using index2
+   $fields = X Y M T R r MT I
+  else
+   $fields = X Y M T R r MT
+  end
+  
+  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 grid.load.sersic.test
   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
-  end
-
+    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
-
-  # 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 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 +794 @@
 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 != 2)
+    echo "USAGE: load.normdata (file)"
     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)
+  lim n dm; clear; box; plot n dm
 end