Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/models/pmModel_SERSIC.c	(revision 28781)
@@ -126,5 +126,6 @@
         dPAR[PM_PAR_SKY]  = +1.0;
         dPAR[PM_PAR_I0]   = +f1;
-        dPAR[PM_PAR_7]    = (z < 0.01) ? -1.5*z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -1.5*z0*f2*log(z);
+        dPAR[PM_PAR_7]    = (z < 0.01) ? -z0*pow(0.01,PAR[PM_PAR_7])*log(0.01) : -z0*f2*log(z);
+	dPAR[PM_PAR_7]   *= 3.0;
 
         assert (isfinite(z1));
@@ -133,5 +134,5 @@
         dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]);
         dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]);
-        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX];
+        dPAR[PM_PAR_SXX]  = +2.0*z1*px*px/PAR[PM_PAR_SXX]; // XXX : increase drag?
         dPAR[PM_PAR_SYY]  = +2.0*z1*py*py/PAR[PM_PAR_SYY];
         dPAR[PM_PAR_SXY]  = -1.0*z1*X*Y;
@@ -224,4 +225,15 @@
     psF32     *PAR  = model->params->data.F32;
 
+    // the other parameters depend on the guess for PAR_7
+    if (!isfinite(PAR[PM_PAR_7])) {
+	PAR[PM_PAR_7]    = 0.25;
+    }    
+    float index = 0.5 / PAR[PM_PAR_7];
+
+    // the scale-length is a function of the moments and the index:
+    // Rmajor_guess = Rmajor_moments * Scale * Zero
+    float Scale = 0.70 + 0.053 * PAR[PM_PAR_7];
+    float Zero  = 1.16 - 0.615 * PAR[PM_PAR_7];
+
     psEllipseMoments emoments;
     emoments.x2 = moments->Mxx;
@@ -236,4 +248,12 @@
     if (!isfinite(axes.theta)) return false;
 
+    // set a lower limit to avoid absurd solutions..
+    float Rmajor = PS_MAX(1.0, Scale * axes.major + Zero);
+    float Rminor = Rmajor * (axes.minor / axes.major);
+
+    fprintf (stderr, "guess index: %f : %f, %f -> %f, %f\n", index, axes.major, axes.minor, Rmajor, Rminor);
+
+    axes.major = Rmajor;
+    axes.minor = Rminor;
     psEllipseShape shape = psEllipseAxesToShape (axes);
 
@@ -242,15 +262,11 @@
     if (!isfinite(shape.sxy)) return false;
 
-    // the other parameters depend on the guess for PAR_7
-    if (!isfinite(PAR[PM_PAR_7])) {
-	PAR[PM_PAR_7]    = 0.25;
-	// PAR[PM_PAR_7]    = 0.125;
-	// PAR[PM_PAR_7]    = 0.5;
-    }    
-
-    float index = 0.5 / PAR[PM_PAR_7];
     float bn = 1.9992*index - 0.3271;
     // float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
     float Io = exp(0.5*bn);
+
+    // XXX do we need this factor of sqrt(2)?
+    // float Sxx = PS_MAX(0.5, M_SQRT2*shape.sx);
+    // float Syy = PS_MAX(0.5, M_SQRT2*shape.sy);
 
     float Sxx = PS_MAX(0.5, shape.sx);
@@ -329,4 +345,6 @@
     psF64 zn = log(PAR[PM_PAR_I0] / flux);
     psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / PAR[PM_PAR_7]);
+
+    fprintf (stderr, "sersic model %f %f, n %f, radius: %f, zn: %f, f/Io: %f, major: %f\n", PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], PAR[PM_PAR_7], radius, zn, flux/PAR[PM_PAR_I0], axes.major);
 
     psAssert (isfinite(radius), "fix this code: z should not be nan for %f", PAR[PM_PAR_7]);
Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c	(revision 28781)
@@ -256,2 +256,107 @@
     return pcm;
 }
+
+// has the set of fitted terms changed?  has the fitting radius changed?
+bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model) {
+
+    bool newWindow = (source->pixels->numRows != pcm->modelFlux->numRows) || (source->pixels->numCols != pcm->modelFlux->numCols);
+
+    // re-count the number of unmasked pixels:
+    if (newWindow) {
+	for (psS32 i = 0; i < source->pixels->numRows; i++) {
+	    for (psS32 j = 0; j < source->pixels->numCols; j++) {
+		// XXX are we doing the right thing with the mask?
+		// skip masked points
+		if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
+		    continue;
+		}
+		// skip zero-variance points
+		if (source->variance->data.F32[i][j] == 0) {
+		    continue;
+		}
+		// skip nan value points
+		if (!isfinite(source->pixels->data.F32[i][j])) {
+		    continue;
+		}
+		pcm->nPix ++;
+	    }
+	}    
+    }
+
+    // if we changed the fit mode, we need to update nDOF
+    int nParams = 0;
+    // set parameter mask based on fitting mode
+    switch (fitOptions->mode) {
+      case PM_SOURCE_FIT_NORM:
+	// NORM-only model fits only source normalization (Io)
+	nParams = 1;
+	psVectorInit (pcm->constraint->paramMask, 1);
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+	break;
+      case PM_SOURCE_FIT_PSF:
+	// PSF model only fits x,y,Io
+	nParams = 3;
+	psVectorInit (pcm->constraint->paramMask, 1);
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
+	break;
+      case PM_SOURCE_FIT_EXT:
+	// EXT model fits all params (except sky)
+	nParams = model->params->n - 1;
+	psVectorInit (pcm->constraint->paramMask, 0);
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+	break;
+      case PM_SOURCE_FIT_INDEX:
+	// PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
+	psVectorInit (pcm->constraint->paramMask, 1);
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+	if (model->params->n == 7) {
+	    nParams = 1;
+	} else {
+	    nParams = 2;
+	    pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
+	}
+	break;
+      case PM_SOURCE_FIT_NO_INDEX:
+	// PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
+	psVectorInit (pcm->constraint->paramMask, 0);
+	pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+	if (model->params->n == 7) {
+	    nParams = model->params->n - 1;
+	} else {
+	    nParams = model->params->n - 2;
+	    pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+	}
+	break;
+      default:
+	psAbort("invalid fitting mode");
+    }
+
+    if (pcm->nPix <  nParams + 1) {
+        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
+        model->flags |= PM_MODEL_STATUS_BADARGS;
+	return false;
+    }
+    pcm->nDOF = pcm->nPix - nParams - 1;
+
+    // has the source pixel window changed?
+    if (newWindow) {
+
+	// adjust all supporting images:
+	pcm->modelFlux = psImageCopy (pcm->modelFlux, source->pixels, PS_TYPE_F32);
+	for (psS32 n = 0; n < pcm->dmodelsFlux->n; n++) {
+	    if ((pcm->constraint->paramMask != NULL) && (pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
+	    pcm->dmodelsFlux->data[n] = psImageCopy (pcm->dmodelsFlux->data[n], source->pixels, PS_TYPE_F32);
+	}
+
+	// adjust images for convolved model and derivative images
+	pcm->modelConvFlux = psImageCopy (pcm->modelConvFlux, source->pixels, PS_TYPE_F32);
+	for (psS32 n = 0; n < pcm->dmodelsConvFlux->n; n++) {
+	    if ((pcm->constraint->paramMask != NULL) && (pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
+	    pcm->dmodelsConvFlux->data[n] = psImageCopy (pcm->dmodelsConvFlux->data[n], source->pixels, PS_TYPE_F32);
+	}
+    }
+
+    return true;
+}
Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.h	(revision 28781)
@@ -62,4 +62,5 @@
 
 pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model, psImageMaskType maskVal, float psfSize);
+bool pmPCMupdate(pmPCMdata *pcm, pmSource *source, pmSourceFitOptions *fitOptions, pmModel *model);
 
 psImage *pmPCMdataSaveImage (pmPCMdata *pcm);
Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.c	(revision 28781)
@@ -966,5 +966,6 @@
     bool addNoise = mode & PM_MODEL_OP_NOISE;
 
-    if (source->modelFlux) {
+    // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
+    if (!addNoise && source->modelFlux) {
         // add in the pixels from the modelFlux image
         int dX = source->modelFlux->col0 - source->pixels->col0;
@@ -987,8 +988,4 @@
 
         psF32 **target = source->pixels->data.F32;
-        if (addNoise) {
-	    // when adding noise, we assume the shape and Io have been modified
-            target = source->variance->data.F32;
-        }
 
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
@@ -1005,10 +1002,8 @@
             }
         }
-	if (!addNoise) {
-	    if (add) {
-		source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-	    } else {
-		source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
-	    }
+	if (add) {
+	    source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	} else {
+	    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
 	}
         return true;
@@ -1029,4 +1024,78 @@
 	}
     }
+
+    return true;
+}
+
+// should we call pmSourceCacheModel if it does not exist?
+bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy)
+{
+    assert (mode & PM_MODEL_OP_NOISE);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+
+    if (add) {
+        psTrace ("psphot", 3, "adding noise to object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+        psTrace ("psphot", 3, "removing noise from object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmSourceNoiseOpModel (source->modelPSF, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
+
+    if (source->modelEXT) {
+	pmSourceNoiseOpModel (source->modelEXT, source, mode, FACTOR, SIZE, add, maskVal, dx, dy);
+    }
+
+    return true;
+}
+
+bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy) 
+{
+    bool status;
+    psEllipseShape oldshape;
+    psEllipseShape newshape;
+    psEllipseAxes axes;
+
+    if (add) {
+	psTrace ("psphot", 4, "adding noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    } else {
+	psTrace ("psphot", 4, "remove noise for object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    }
+
+    psF32 *PAR = model->params->data.F32;
+
+    // save original values
+    float oldI0  = PAR[PM_PAR_I0];
+    oldshape.sx  = PAR[PM_PAR_SXX];
+    oldshape.sy  = PAR[PM_PAR_SYY];
+    oldshape.sxy = PAR[PM_PAR_SXY];
+
+    // XXX can this be done more intelligently?
+    if (oldI0 == 0.0) return false;
+    if (!isfinite(oldI0)) return false;
+
+    // increase size and height of source
+    axes = psEllipseShapeToAxes (oldshape, 20.0);
+    axes.major *= SIZE;
+    axes.minor *= SIZE;
+    newshape = psEllipseAxesToShape (axes);
+    PAR[PM_PAR_I0]  = FACTOR*oldI0;
+    PAR[PM_PAR_SXX] = newshape.sx;
+    PAR[PM_PAR_SYY] = newshape.sy;
+    PAR[PM_PAR_SXY] = newshape.sxy;
+
+    psImage *target = source->variance;
+
+    if (add) {
+	status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    } else {
+	status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    }
+
+    // restore original values
+    PAR[PM_PAR_I0]  = oldI0;
+    PAR[PM_PAR_SXX] = oldshape.sx;
+    PAR[PM_PAR_SYY] = oldshape.sy;
+    PAR[PM_PAR_SXY] = oldshape.sxy;
 
     return true;
Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSource.h	(revision 28781)
@@ -236,4 +236,7 @@
 bool pmSourceSubWithOffset (pmSource *source, pmModelOpMode mode, psImageMaskType maskVal, int dx, int dy);
 
+bool pmSourceNoiseOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
+bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
+
 bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
 bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);
Index: /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c
===================================================================
--- /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c	(revision 28780)
+++ /branches/eam_branches/ipp-20100621/psModules/src/objects/pmSourceFitPCM.c	(revision 28781)
@@ -107,5 +107,8 @@
 bool pmSourceModelGuessPCM (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    pcm->modelConv->modelGuess(pcm->modelConv, source);
+    if (!pcm->modelConv->modelGuess(pcm->modelConv, source)) {
+	return false;
+    }
+    return true;
 
     // generate copy of the model
@@ -114,9 +117,14 @@
 
     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
+    // SXX,SYY in model parameters are sqrt(2) * shape.Sxx,Syy
     psEllipseShape psfShape;
     psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
     psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
     psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+
     psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
+    if (!isfinite(psfAxes.major)) return false;
+    if (!isfinite(psfAxes.minor)) return false;
+    if (!isfinite(psfAxes.theta)) return false;
 
     // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
@@ -125,13 +133,21 @@
     extShape.sxy = pcm->modelConv->params->data.F32[PM_PAR_SXY];
     extShape.sy  = pcm->modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+
     psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
+    if (!isfinite(extAxes.major)) return false;
+    if (!isfinite(extAxes.minor)) return false;
+    if (!isfinite(extAxes.theta)) return false;
 
     // decrease the initial guess ellipse by psf_minor axis:
     psEllipseAxes extAxesMod;
-    extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
-    extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
+    extAxesMod.major = sqrt (PS_MAX (0.25, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
+    extAxesMod.minor = sqrt (PS_MAX (0.25, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
     extAxesMod.theta = extAxes.theta;
 
     psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
+    if (!isfinite(extShapeMod.sx))  return false;
+    if (!isfinite(extShapeMod.sy))  return false;
+    if (!isfinite(extShapeMod.sxy)) return false;
+
     pcm->modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
     pcm->modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;