Index: trunk/psModules/src/objects/pmPCM_MinimizeChisq.c
===================================================================
--- trunk/psModules/src/objects/pmPCM_MinimizeChisq.c	(revision 34403)
+++ trunk/psModules/src/objects/pmPCM_MinimizeChisq.c	(revision 35768)
@@ -81,8 +81,9 @@
         psAbort ("programming error: no unmasked parameters to be fit\n");
     }
+    psAssert (pcm->nPar == Beta->n, "did we set the masked parameters correctly??");
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
+    psImage *alpha   = psImageAlloc(pcm->nPar, pcm->nPar, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(pcm->nPar, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -90,4 +91,5 @@
     psF32 lambda = 0.001;
     psF32 dLinear = 0.0;
+    psF32 nu = 2.0;
 
 # if (USE_FFT && PRE_CONVOLVE)
@@ -120,6 +122,11 @@
     bool done = (min->iter >= min->maxIter);
     while (!done) {
-        psTrace("psModules.objects", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psModules.objects", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        psTrace(FACILITY, 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+
+	if (min->chisqConvergence) {
+	    psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+	} else {
+	    psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*pcm->nPar, min->maxTol*pcm->nPar);
+	}
 
         // set a new guess for Alpha, Beta, Params
@@ -140,4 +147,27 @@
             p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
         }
+
+	// calculate the parameter change (rParDelta) and error radius (rParSigma)
+	//    rParDelta : radius of parameter change;
+	//    rParSigma : radius of parameter error 
+	
+	// note that (before SetABX) Alpha[i][i] is the covariance matrix and
+	// Beta is the actual parameter change for this pass
+
+	// note that Alpha & Beta only represent unmasked parameters, while params and Params have all 
+
+	// dParSigma = Alpha[i][i] : error (squared) on parameter i
+	// dParDelta = Params->data.F32[i] - params->data.F32[i]     : change on parameter i
+	float rParSigma = 0.0;
+        for (int j = 0, J = 0; j < Params->n; j++) {
+	    if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+		continue;
+	    }
+	    rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J];
+	    J++;
+	}
+	rParSigma = sqrt(rParSigma);
+	psTrace(FACILITY, 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+	// fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
 
         // calculate Chisq for new guess, update Alpha & Beta
@@ -164,27 +194,77 @@
         }
 
-        /* if (Chisq < min->value) {  */
+	// change in chisq/nDOF since last minimum
+	min->lastDelta = (min->value - Chisq) / pcm->nDOF;
+
+        // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho)
+
+	// XXX the old version of lambda changes:
+	// XXX : Madsen gives suggestion for better use of rho
+        // rho is positive if the new chisq is smaller
         if (rho >= -1e-6) {
-            min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n);
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
             beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
             params = psVectorCopy(params, Params, PS_TYPE_F32);
-            lambda *= 0.25;
 
             // save the new convolved model image
             psFree (source->modelFlux);
             source->modelFlux = pmPCMdataSaveImage(pcm);
-        } else {
-            lambda *= 10.0;
-        }
+        } 
+	switch (min->gainFactorMode) {
+	  case 0:
+	    if (rho >= -1e-6) {
+		lambda *= 0.25;
+	    } else {
+		lambda *= 10.0;
+	    }
+	    break;
+
+	  case 1:
+	    // adjust the gain ratio (lambda) based on rho
+	    if (rho < 0.25) {
+		lambda *= 2.0;
+	    } 
+	    if (rho > 0.75) {
+		lambda *= 0.333;
+	    }
+	    break;
+
+	  case 2:
+	    if (rho > 0.0) {
+		lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
+		nu = 2.0;
+	    } else {
+		lambda *= nu;
+		nu *= 2.0;
+	    }
+	    break;
+	}
         min->iter++;
 
+	// ending conditions:
+	// 1) hard limit : too many iterations
 	done = (min->iter >= min->maxIter);
 	
-	// check for convergence:
+	// 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
+	if (min->lastDelta < -1e-6) {
+	    continue;
+	}
+
+	// save this value in case we stop iterating
+	min->rParSigma = rParSigma;
+
+	// 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN)
+	// keep iterating regardless of rParSigma in this case
 	float chisqDOF = Chisq / pcm->nDOF;
-	if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+	if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
+	    continue;
+	}
+
+	// delta-chisq or rParSigma ?
+	if (min->chisqConvergence) {
 	    done |= (min->lastDelta < min->minTol);
+	} else {
+	    done |= (rParSigma < min->minTol*pcm->nPar);
 	}
     }
@@ -220,7 +300,14 @@
 
     // if the last improvement was at least as good as maxTol, accept the fit:
-    if (min->lastDelta <= min->maxTol) {
-	psTrace(FACILITY, 6, "---- end (true) ----\n");
-        return(true);
+    if (min->chisqConvergence) {
+	if (min->lastDelta <= min->maxTol) {
+	    psTrace(FACILITY, 6, "---- end (true) ----\n");
+	    return(true);
+	}
+    } else {
+	if (min->rParSigma <= min->maxTol*pcm->nPar) {
+	    psTrace(FACILITY, 6, "---- end (true) ----\n");
+	    return(true);
+	}
     }
     psTrace(FACILITY, 6, "---- end (false) ----\n");