Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 13982)
+++ /trunk/psphot/src/Makefile.am	(revision 13983)
@@ -17,38 +17,46 @@
 
 libpsphot_la_SOURCES = \
-	psphotErrorCodes.c	\
-	pmFootprint.c		\
-	psphotVersion.c		\
-	psphotModelGroupInit.c	\
-	psphotMaskReadout.c	\
-	psphotDefineFiles.c	\
-	psphotReadout.c		\
-	psphotImageMedian.c	\
-	psphotFindPeaks.c	\
-	psphotSourceStats.c	\
-	psphotRoughClass.c	\
-	psphotBasicDeblend.c	\
-	psphotChoosePSF.c	\
-	psphotGuessModels.c     \
+	psphotErrorCodes.c	 \
+	pmFootprint.c		 \
+	psphotVersion.c		 \
+	psphotModelGroupInit.c	 \
+	psphotMaskReadout.c	 \
+	psphotDefineFiles.c	 \
+	psphotReadout.c		 \
+	psphotImageMedian.c	 \
+	psphotFindPeaks.c	 \
+	psphotSourceStats.c	 \
+	psphotRoughClass.c	 \
+	psphotBasicDeblend.c	 \
+	psphotChoosePSF.c	 \
+	psphotGuessModels.c      \
 	psphotFitSourcesLinear.c \
-	psphotBlendFit.c	\
-	psphotReplaceUnfit.c	\
-	psphotApResid.c		\
-	psphotMagnitudes.c	\
-	psphotSkyReplace.c	\
-	psphotEvalPSF.c		\
-	psphotEvalFLT.c		\
-	psphotSourceFits.c	\
-	psphotRadiusChecks.c	\
-	psphotSortBySN.c	\
-	psphotOutput.c		\
-	psphotGrowthCurve.c	\
-	psphotFakeSources.c	\
-	psphotModelTest.c	\
-	psphotFitSet.c		\
+	psphotBlendFit.c	 \
+	psphotReplaceUnfit.c	 \
+	psphotApResid.c		 \
+	psphotMagnitudes.c	 \
+	psphotSkyReplace.c	 \
+	psphotEvalPSF.c		 \
+	psphotEvalFLT.c		 \
+	psphotSourceFits.c	 \
+	psphotRadiusChecks.c	 \
+	psphotSortBySN.c	 \
+	psphotOutput.c		 \
+	psphotGrowthCurve.c	 \
+	psphotFakeSources.c	 \
+	psphotModelWithPSF.c     \
+	psphotExtendedSources.c	 \
+	psphotPetrosian.c	 \
+	psphotIsophotal.c	 \
+	psphotAnnuli.c		 \
+	psphotKron.c		 \
+	psphotKernelFromPSF.c	 \
+	psphotPSFConvModel.c	 \
+	psphotModelTest.c	 \
+	psphotFitSet.c		 \
 	psphotSourceFreePixels.c \
 	psphotSummaryPlots.c     \
 	psphotMergeSources.c	 \
-	psphotLoadPSF.c	 \
+	psphotLoadPSF.c	         \
 	psphotReadoutCleanup.c	 \
 	psphotSourcePlots.c	 \
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 13982)
+++ /trunk/psphot/src/psphot.h	(revision 13983)
@@ -114,3 +114,32 @@
 bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal);
 
+bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal);
+bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal);
+psKernel *psphotKernelFromPSF (pmSource *source);
+
+bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal);
+bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal);
+bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal);
+bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal);
+
+// psf-convolved model fitting
+bool psphotModelWithPSF_LMM (
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    psMinConstraint *constraint,
+    const pmSource *source,
+    const psKernel *psf,
+    psMinimizeLMChi2Func func);
+
+psF32 psphotModelWithPSF_SetABX(
+    psImage  *alpha,
+    psVector *beta,
+    const psVector *params,
+    const psVector *paramMask,
+    const pmSource *source,
+    const psKernel *psf,
+    psMinimizeLMChi2Func func);
+
+
 #endif
Index: /trunk/psphot/src/psphotAnnuli.c
===================================================================
--- /trunk/psphot/src/psphotAnnuli.c	(revision 13983)
+++ /trunk/psphot/src/psphotAnnuli.c	(revision 13983)
@@ -0,0 +1,8 @@
+# include "psphot.h"
+
+bool psphotAnnuli (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+
+  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  return true;
+
+}
Index: /trunk/psphot/src/psphotExtendedSources.c
===================================================================
--- /trunk/psphot/src/psphotExtendedSources.c	(revision 13983)
+++ /trunk/psphot/src/psphotExtendedSources.c	(revision 13983)
@@ -0,0 +1,113 @@
+# include "psphot.h"
+
+bool psphotExtendedSources (pmReadout *readout, psArray *sources, psMetadata *recipe, psMaskType maskVal) {
+
+    bool status;
+    int Next = 0;
+    int Npsf = 0;
+
+    // S/N limit to perform full non-linear fits
+    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
+	psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
+	return true;
+    }
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // 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 doPSFConvModel = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PSF_CONVOLVED_MODEL");
+    bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
+    bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+    bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // XXX some init functions for the extended source recipe options?
+
+    // choose the sources of interest
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        // skip PSF-like and non-astronomical objects
+        if (source->type == PM_SOURCE_TYPE_STAR) continue;
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+        // limit selection to some SN limit
+        assert (source->peak); // how can a source not have a peak?
+        if (source->peak->SN < SN_LIM) continue;
+
+        // XXX this should use peak?
+        if (source->peak->x < AnalysisRegion.x0) continue;
+        if (source->peak->y < AnalysisRegion.y0) continue;
+        if (source->peak->x > AnalysisRegion.x1) continue;
+        if (source->peak->y > AnalysisRegion.y1) continue;
+
+        // if model is NULL, we don't have a starting guess
+        if (source->modelEXT == NULL) continue;
+
+        // replace object in image
+        if (source->mode & PM_SOURCE_MODE_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        }
+        Next ++;
+
+	if (doPSFConvModel && !psphotPSFConvModel (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in PSF Convolved Model fit");
+	    return false;
+	} else {
+	    pmSourceCacheModel (source, maskVal); // XXX put this in the source model function?
+	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            Npsf ++;
+	}
+
+	if (doPetrosian && !psphotPetrosian (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Petrosian analysis");
+	    return false;
+	} else {
+	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            Npsf ++;
+	}
+
+	if (doIsophotal && !psphotIsophotal (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Isophotal analysis");
+	    return false;
+	} else {
+	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            Npsf ++;
+	}
+
+	if (doAnnuli && !psphotAnnuli (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Annuli analysis");
+	    return false;
+	} else {
+	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            Npsf ++;
+	}
+
+	if (doKron && !psphotKron (source, recipe, maskVal)) {
+	    psError(PSPHOT_ERR_UNKNOWN, false, "failure in Kron analysis");
+	    return false;
+	} else {
+	    psTrace ("psphot", 5, "psf-convolved model for source at %7.1f, %7.1f", source->moments->x, source->moments->y);
+            Npsf ++;
+	}
+
+        // re-subtract the object, leave local sky
+        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        source->mode |= PM_SOURCE_MODE_SUBTRACTED;
+        source->mode |= PM_SOURCE_MODE_TEMPSUB;
+    }
+
+    psLogMsg ("psphot", PS_LOG_INFO, "extended source analysis: %f sec for %d objects\n", psTimerMark ("psphot"), Next);
+    return true;
+}
Index: /trunk/psphot/src/psphotIsophotal.c
===================================================================
--- /trunk/psphot/src/psphotIsophotal.c	(revision 13983)
+++ /trunk/psphot/src/psphotIsophotal.c	(revision 13983)
@@ -0,0 +1,8 @@
+# include "psphot.h"
+
+bool psphotIsophotal (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+
+  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  return true;
+
+}
Index: /trunk/psphot/src/psphotKernelFromPSF.c
===================================================================
--- /trunk/psphot/src/psphotKernelFromPSF.c	(revision 13983)
+++ /trunk/psphot/src/psphotKernelFromPSF.c	(revision 13983)
@@ -0,0 +1,20 @@
+# include "psphot.h"
+
+psKernel *psphotKernelFromPSF (pmSource *source) {
+
+  assert (source);
+  assert (source->psfFlux); // XXX build if needed?
+
+  int x0 = source->peak->xf - source->psfFlux->col0;
+  int y0 = source->peak->yf - source->psfFlux->row0;
+
+  // need to decide on the size: dynamically? statically?
+  psKernel *psf = psKernelAlloc (-1, +1, -1, +1);
+
+  for (int j = psf->yMin; j <= psf->yMax; j++) {
+    for (int i = psf->xMin; i <= psf->xMax; i++) {
+      psf->kernel[j][i] = source->psfFlux->data.F32[y0 + j][x0 + i];
+    }
+  }
+  return psf;
+}
Index: /trunk/psphot/src/psphotKron.c
===================================================================
--- /trunk/psphot/src/psphotKron.c	(revision 13983)
+++ /trunk/psphot/src/psphotKron.c	(revision 13983)
@@ -0,0 +1,8 @@
+# include "psphot.h"
+
+bool psphotKron (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+
+  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  return true;
+
+}
Index: /trunk/psphot/src/psphotModelWithPSF.c
===================================================================
--- /trunk/psphot/src/psphotModelWithPSF.c	(revision 13982)
+++ /trunk/psphot/src/psphotModelWithPSF.c	(revision 13983)
@@ -1,6 +1,3 @@
 # include "psphot.h"
-
-// XXX elevate the p_psMinLM_ functions to psMinLM_...
-
 
 bool psphotModelWithPSF_LMM (
@@ -69,5 +66,5 @@
 
         // set a new guess for Alpha, Beta, Params
-        if (!p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)) {
             min->iter ++;
             lambda *= 10.0;
@@ -76,5 +73,5 @@
 
         // measure linear model prediction
-        psF32 dLinear = p_psMinLM_dLinear(Beta, beta, lambda);
+        psF32 dLinear = psMinLM_dLinear(Beta, beta, lambda);
 
         // dump some useful info if trace is defined
@@ -88,5 +85,5 @@
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = p_psMinLM_SetABX(Alpha, Beta, Params, paramMask, source, func);
+        Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, source, psf, func);
         if (isnan(Chisq)) {
             min->iter ++;
@@ -112,5 +109,5 @@
         /* if (Chisq < min->value) {  */
         if (rho > 0.0) {
-            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n);
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -127,5 +124,5 @@
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
-        if (!p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0)) {
+        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0)) {
             psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
         }
@@ -172,8 +169,4 @@
     }
 
-    psF32 chisq;
-    psF32 delta;
-    psF32 weight;
-    psF32 ymodel;
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -191,5 +184,4 @@
 
     // fill in the coordinate and value entries
-    nPix = 0;
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
@@ -239,5 +231,5 @@
     psImageInit (alpha, 0.0);
     psVectorInit (beta, 0.0);
-    chisq = 0.0;
+    float chisq = 0.0;
 
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
@@ -257,9 +249,9 @@
             }
 
-	    ymodel  = modelConv->data.F32[i][j];
-	    yweight = 1.0 / source->weight->data.F32[i][j];
-	    delta = ymodel - source->pixels->data.F32[i][j];
-
-	    chisq += PS_SQR(delta) * var;
+	    float ymodel  = modelConv->data.F32[i][j];
+	    float yweight = 1.0 / source->weight->data.F32[i][j];
+	    float delta = ymodel - source->pixels->data.F32[i][j];
+
+	    chisq += PS_SQR(delta) * yweight;
 
 	    if (isnan(delta))
@@ -273,5 +265,5 @@
 	      }
 	      psImage *dmodel = dmodelsConv->data[n1];
-	      weight = dmodel->data.F32[i][j] * yweight;
+	      float weight = dmodel->data.F32[i][j] * yweight;
 	      for (psS32 n2 = 0; n2 <= n1; n2++) {
                 if ((paramMask != NULL) && (paramMask->data.U8[n2])) {
@@ -281,5 +273,5 @@
                 alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j];
 	      }
-	      beta->data.F32[jn] += weight * delta;
+	      beta->data.F32[n1] += weight * delta;
 	    }
 	}
Index: /trunk/psphot/src/psphotPSFConvModel.c
===================================================================
--- /trunk/psphot/src/psphotPSFConvModel.c	(revision 13983)
+++ /trunk/psphot/src/psphotPSFConvModel.c	(revision 13983)
@@ -0,0 +1,95 @@
+# include "psphot.h"
+
+// save as static values so they may be set externally
+static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
+static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
+// static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
+
+bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+    
+    // XXX make sure we save a cached copy of the source psf
+
+    // convert the cached cached psf model for this source to a psKernel
+    psKernel *psf = psphotKernelFromPSF (source);
+
+    // generate copy of the model
+    // XXX we could modify the parameter values or even the model 
+    // here based on the observed seeing (some lookup table...)
+    pmModel *modelConv = pmModelCopy (source->modelEXT);
+    psVector *params  = modelConv->params;
+    psVector *dparams = modelConv->dparams;
+
+    // get the model function for this model
+    pmModelFunc modelFunc = pmModelFunc_GetFunction (modelConv->type);
+    if (!modelFunc)
+	psAbort("invalid model function");
+
+    // get the limits function for this model
+    pmModelLimits checkLimits = pmModelLimits_GetFunction (modelConv->type);
+    if (!checkLimits)
+	psAbort("invalid model limits function");
+
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+    constraint->checkLimits = checkLimits;
+
+    // set parameter mask based on fitting mode
+    // we fit a model without a floating sky term
+    int nParams = params->n - 1;
+    psVectorInit (constraint->paramMask, 0);
+    constraint->paramMask->data.U8[PM_PAR_SKY] = 1;
+
+    // force the floating parameters to fall within the contraint ranges
+    for (int i = 0; i < params->n; i++) {
+	checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+	checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+    }
+
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
+
+    bool fitStatus = psphotModelWithPSF_LMM (myMin, covar, params, constraint, source, psf, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+        if (psTraceGetLevel("psModules.objects") >= 4) {
+            fprintf (stderr, "%f ", params->data.F32[i]);
+        }
+        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
+    }
+    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
+
+    // save the resulting chisq, nDOF, nIter
+    modelConv->chisq = myMin->value;
+    modelConv->nIter = myMin->iter;
+
+    // XXX I actually need to count the number of unmasked pixels here
+    modelConv->nDOF  = source->pixels->numCols*source->pixels->numRows  -  nParams;
+
+    modelConv->flags |= PM_MODEL_STATUS_FITTED;
+    if (!fitStatus) modelConv->flags |= PM_MODEL_STATUS_NONCONVERGE;
+
+    // models can go insane: reject these
+    bool onPic = true;
+    onPic &= (params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
+    onPic &= (params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0);
+    onPic &= (params->data.F32[PM_PAR_YPOS] <  source->pixels->row0 + source->pixels->numRows);
+    if (!onPic) {
+        modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
+    }
+
+    source->mode |= PM_SOURCE_MODE_FITTED;
+    source->modelConv = modelConv;
+
+    psFree(myMin);
+    psFree(covar);
+    psFree(constraint);
+
+    bool retval = (onPic && fitStatus);
+    psTrace("psModules.objects", 5, "---- %s(%d) end ----\n", __func__, retval);
+    return(retval);
+}
Index: /trunk/psphot/src/psphotPetrosian.c
===================================================================
--- /trunk/psphot/src/psphotPetrosian.c	(revision 13983)
+++ /trunk/psphot/src/psphotPetrosian.c	(revision 13983)
@@ -0,0 +1,8 @@
+# include "psphot.h"
+
+bool psphotPetrosian (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+
+  psLogMsg ("psphot", PS_LOG_INFO, "not implemented\n");
+  return true;
+
+}
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 13982)
+++ /trunk/psphot/src/psphotReadout.c	(revision 13983)
@@ -246,4 +246,6 @@
     if (dump) psphotSaveImage (NULL, readout->image,  "image.v6.fits");
 
+    psphotExtendedSources (readout, sources, recipe, maskVal);
+
 finish: