Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_PS1_V1.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 35639)
@@ -1,9 +1,10 @@
 /******************************************************************************
- * this file defines the PS1_V1 source shape model.  Note that these model functions are loaded
- * by pmModelClass.c using 'include', and thus need no 'include' statements of their own.  The
- * models use a psVector to represent the set of parameters, with the sequence used to specify
- * the meaning of the parameter.  The meaning of the parameters may thus vary depending on the
- * specifics of the model.  All models which are used as a PSF representations share a few
- * parameters, for which # define names are listed in pmModel.h:
+ * this file defines the PS1_V1 source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 35639)
@@ -1,9 +1,10 @@
 /******************************************************************************
- * this file defines the QGAUSS source shape model (XXX need a better name!).  Note that these
- * model functions are loaded by pmModelClass.c using 'include', and thus need no 'include'
- * statements of their own.  The models use a psVector to represent the set of parameters, with
- * the sequence used to specify the meaning of the parameter.  The meaning of the parameters
- * may thus vary depending on the specifics of the model.  All models which are used a PSF
- * representations share a few parameters, for which # define names are listed in pmModel.h:
+ * this file defines the QGAUSS source shape model.  Note that these model functions are
+ * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of
+ * their own.  The models use a psVector to represent the set of parameters, with the
+ * sequence used to specify the meaning of the parameter.  The meaning of the parameters
+ * may thus vary depending on the specifics of the model.  All models which are used as a
+ * PSF representations share a few parameters, for which # define names are listed in
+ * pmModel.h:
 
    power-law with fitted linear term
@@ -14,6 +15,6 @@
    * PM_PAR_XPOS 2  - X center of object
    * PM_PAR_YPOS 3  - Y center of object
-   * PM_PAR_SXX 4   - X^2 term of elliptical contour (sqrt(2) / SigmaX)
-   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (sqrt(2) / SigmaY)
+   * PM_PAR_SXX 4   - X^2 term of elliptical contour (SigmaX / sqrt(2))
+   * PM_PAR_SYY 5   - Y^2 term of elliptical contour (SigmaY / sqrt(2))
    * PM_PAR_SXY 6   - X*Y term of elliptical contour
    * PM_PAR_7   7   - amplitude of the linear component (k)
@@ -138,5 +139,4 @@
 # define AR_MAX 20.0
 # define AR_RATIO 0.99
-
 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta)
 {
@@ -203,5 +203,5 @@
           return true;
       }
-    default:
+      default:
         psAbort("invalid choice for limits");
     }
@@ -293,6 +293,4 @@
     if (flux >= PAR[PM_PAR_I0]) return 1.0;
 
-    // if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
-
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
     shape.sy  = PAR[PM_PAR_SYY] / M_SQRT2;
@@ -307,5 +305,4 @@
         return ( sigma * sqrt (2.0 * z) );
     }
-
     psF64 limit = flux / PAR[PM_PAR_I0];
 
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_TRAIL.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_TRAIL.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/models/pmModel_TRAIL.c	(revision 35639)
@@ -369,7 +369,10 @@
     if (!isfinite(axes.theta)) return false;
 
-    float size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
-    //    if (size > psfAxes.major) { size -= psfAxes.major; }
-    //else { size = psfAxes.major; }
+    float size = NAN;
+    if (!isfinite(source->moments->Mrf)) {
+      size = axes.major;
+    } else {
+      size = (axes.major > sqrt(source->moments->Mrf)) ? axes.major : sqrt(source->moments->Mrf);
+    }
 
     float theta, peak;
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtry.h
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtry.h	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtry.h	(revision 35639)
@@ -99,8 +99,12 @@
 /** fit EXT models to all possible psf sources */
 bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmPSFtryFitEXT_Threaded (psThreadJob *job);
 
 bool pmPSFtryMakePSF (bool *pGoodFit, pmPSFtry *psfTry);
 
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal);
+bool pmPSFtryFitPSF_Threaded (psThreadJob *job);
+
+bool pmPSFThreads (void);
 
 /** pmPSFtryMetric()
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitEXT.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitEXT.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitEXT.c	(revision 35639)
@@ -46,9 +46,26 @@
 #include "pmSourceVisual.h"
 
+bool pmPSFThreads (void) {
+
+    psThreadTask *task = NULL;
+
+    task = psThreadTaskAlloc("PSF_TRY_FIT_EXT", 6);
+    task->function = &pmPSFtryFitEXT_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSF_TRY_FIT_PSF", 6);
+    task->function = &pmPSFtryFitPSF_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    return true;
+}
+
+static int Next = 0;
+
 // Fit an EXT model to all candidates PSF sources.
 // Note: this is independent of the modeled 2D variations in the PSF.
 bool pmPSFtryFitEXT (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
-
-    bool status;
 
     psTimerStart ("psf.fit");
@@ -60,52 +77,74 @@
     maskVal |= markVal;
 
-    int Next = 0;
+    Next = 0;
     for (int i = 0; i < psfTry->sources->n; i++) {
 
         pmSource *source = psfTry->sources->data[i];
-        if (!source->moments) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-        if (!source->moments->nPixels) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
-        // If mask object does not exist, mark the source as bad.
-        // We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
-        // psphot to exit with a fault. 
-        if (source->maskObj == NULL) {
-            psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
 
-        source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
-        if (source->modelEXT == NULL) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
+	if (!source->moments) {
+	    psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+	    psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no moments\n", i, source->peak->x, source->peak->y);
+	    continue;
+	}
+	if (!source->moments->nPixels) {
+	    psTrace ("psModules.objects", 4, "masking %d (%d,%d) : no pixels\n", i, source->peak->x, source->peak->y);
+	    psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+	    continue;
+	}
+	// If mask object does not exist, mark the source as bad.
+	// We cannot proceed with it because psImageMaskPixels leaves an uncleared error code last which causes
+	// psphot to exit with a fault. 
+	if (source->maskObj == NULL) {
+	    psTrace ("psModules.objects", 4, "source %d (%d,%d) : null maskObj\n", i, source->peak->x, source->peak->y);
+	    psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+	    continue;
+	}
 
-        // set object mask to define valid pixels
-	// XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+	source->modelEXT = pmSourceModelGuess (source, options->type, maskVal, markVal);
+	if (source->modelEXT == NULL) {
+	    psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y);
+	    psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+	    continue;
+	}
 
-        // fit model as EXT, not PSF
-        status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
+	// do some actual work on this source
+	psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_EXT");
+	psArrayAdd(job->args, 1, source);
+	psArrayAdd(job->args, 1, psfTry);
+	psArrayAdd(job->args, 1, options);
+	
+	PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
 
-        // clear object mask to define valid pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+	PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+	PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
 
-        // exclude the poor fits
-        if (!status) {
-            psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
-            continue;
-        }
-        Next ++;
+# if (1)
+	if (!psThreadJobAddPending(job)) {
+	    psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+	    return false;
+	}
+# else
+	if (!pmPSFtryFitEXT_Threaded(job)) {
+	    psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+	    return false;
+	}
+	psFree(job);
+# endif
     }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to model psf");
+	return false;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+	// we have no returned data from this operation
+	if (job->args->n < 1) fprintf (stderr, "error with job\n");
+	psFree(job);
+    }
+
     psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, psfTry->sources->n);
     psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, psfTry->sources->n);
@@ -118,2 +157,34 @@
     return true;
 }
+
+bool pmPSFtryFitEXT_Threaded (psThreadJob *job) {
+
+    pmSource *source =      job->args->data[0];
+    pmPSFtry *psfTry =      job->args->data[1];
+    pmPSFOptions *options = job->args->data[2];
+
+    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+
+    // set object mask to define valid pixels
+    // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)?
+    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+
+    // fit model as EXT, not PSF
+    bool status = pmSourceFitModel (source, source->modelEXT, options->fitOptions, maskVal);
+
+    // clear object mask to define valid pixels
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+    // exclude the poor fits
+    if (!status) {
+	psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y);
+	psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL;
+	return true;
+    }
+    Next ++;
+    
+    return true;
+}
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitPSF.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitPSF.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmPSFtryFitPSF.c	(revision 35639)
@@ -43,10 +43,10 @@
 #include "pmSourceVisual.h"
 
+static int Npsf = 0;
+
 // stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
 // be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    bool status;
-
     psTimerStart ("psf.fit");
 
@@ -57,13 +57,5 @@
     maskVal |= markVal;
 
-    // DEBUG code: save the PSF model fit data in detail
-# ifdef DEBUG
-    char filename[64];
-    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
-    FILE *f = fopen (filename, "w");
-    psAssert (f, "failed open");
-# endif
-
-    int Npsf = 0;
+    Npsf = 0;
     for (int i = 0; i < psfTry->sources->n; i++) {
 
@@ -77,56 +69,68 @@
         }
 
-        // set shape for this model based on PSF
+	// set shape for this model based on PSF
 	psFree (source->modelPSF);
-        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
-        if (source->modelPSF == NULL) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-	// PSF fit and aperture mags use different radii
-        source->modelPSF->fitRadius = options->fitRadius;
-        source->apRadius = options->apRadius;
-
-        // set object mask to define valid pixels for PSF model fit
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
-
-        // fit the PSF model to the source
-        status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
-
-        // skip poor fits
-        if (!status) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // set object mask to define valid pixels for APERTURE magnitude
-	if (options->fitRadius != options->apRadius) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-	    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+	source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
+	if (source->modelPSF == NULL) {
+	    psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
+	    psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
+	    return false;
 	}
 
-	// This function calculates the psf and aperture magnitudes
-        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
-        if (!status || isnan(source->apMag)) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // clear object mask to define valid pixels
-	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-
-        psfTry->fitMag->data.F32[i] = source->psfMag;
-        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
-        psfTry->metricErr->data.F32[i] = source->psfMagErr;
-
-	// XXX this did not work: modifies shape of psf too much
-        // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
-
+	// do some actual work on this source
+	psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_PSF");
+	psArrayAdd(job->args, 1, source);
+	psArrayAdd(job->args, 1, psfTry);
+	psArrayAdd(job->args, 1, options);
+	
+	PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
+
+	PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+	PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+
+# if (1)
+	if (!psThreadJobAddPending(job)) {
+	    psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+	    return false;
+	}
+# else
+	if (!pmPSFtryFitPSF_Threaded(job)) {
+	    psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+	    return false;
+	}
+	psFree(job);
+# endif
+    }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+	psError(PS_ERR_UNKNOWN, false, "failure to model psf");
+	return false;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+	// we have no returned data from this operation
+	if (job->args->n < 1) fprintf (stderr, "error with job\n");
+	psFree(job);
+    }
+    psfTry->psf->nPSFstars = Npsf;
+
+    // DEBUG code: save the PSF model fit data in detail
 # ifdef DEBUG
+
+    char filename[64];
+    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
+    FILE *f = fopen (filename, "w");
+    psAssert (f, "failed open");
+
+    for (int i = 0; i < psfTry->sources->n; i++) {
+
+	// skip masked sources
+        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
+
+        pmSource *source = psfTry->sources->data[i];
+
 	fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
 		 source->peak->xf, source->peak->yf, 
@@ -136,12 +140,5 @@
 		 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
 		 source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
-# endif
-
-        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
-        Npsf ++;
-    }
-    psfTry->psf->nPSFstars = Npsf;
-
-# ifdef DEBUG
+    }
     fclose (f);
 # endif
@@ -159,2 +156,56 @@
     return true;
 }
+
+bool pmPSFtryFitPSF_Threaded (psThreadJob *job) {
+
+    pmSource *source =      job->args->data[0];
+    pmPSFtry *psfTry =      job->args->data[1];
+    pmPSFOptions *options = job->args->data[2];
+
+    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+
+    // PSF fit and aperture mags use different radii
+    source->modelPSF->fitRadius = options->fitRadius;
+    source->apRadius            = options->apRadius;
+
+    // set object mask to define valid pixels for PSF model fit
+    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+
+    // fit the PSF model to the source
+    bool status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
+
+    // skip poor fits
+    if (!status) {
+	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+	psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
+	psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
+	return true;
+    }
+
+    // set object mask to define valid pixels for APERTURE magnitude
+    if (options->fitRadius != options->apRadius) {
+	psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+	psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+    }
+
+    // This function calculates the psf and aperture magnitudes
+    status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+    if (!status || isnan(source->apMag)) {
+	psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
+	psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
+	return true;
+    }
+
+    psfTry->fitMag->data.F32[i] = source->psfMag;
+    psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
+    psfTry->metricErr->data.F32[i] = source->psfMagErr;
+
+    psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
+    Npsf ++;
+    return true;
+}
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.c	(revision 35639)
@@ -62,4 +62,8 @@
     opt->poissonErrors = true;
     opt->saveCovariance = false;
+
+    // we default to the old algorithm
+    opt->gainFactorMode = 0;
+    opt->chisqConvergence = true;
 
     return opt;
@@ -241,4 +245,6 @@
 
     psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
+    myMin->gainFactorMode = options->gainFactorMode;
+    myMin->chisqConvergence = options->chisqConvergence;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -274,5 +280,10 @@
     model->flags |= PM_MODEL_STATUS_FITTED;
     if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
-    if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+
+    if (myMin->chisqConvergence) {
+      if (myMin->lastDelta > myMin->minTol) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    } else {
+      if (myMin->rParSigma > myMin->minTol*nParams) model->flags |= PM_MODEL_STATUS_WEAK_FIT;
+    }
 
     // get the Gauss-Newton distance for fixed model parameters
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.h
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.h	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitModel.h	(revision 35639)
@@ -35,4 +35,6 @@
     bool poissonErrors;			///< use poisson errors for fits?
     bool saveCovariance;
+    int gainFactorMode;
+    bool chisqConvergence; 
 } pmSourceFitOptions;
 
Index: /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitSet.c
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitSet.c	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/objects/pmSourceFitSet.c	(revision 35639)
@@ -568,4 +568,6 @@
 
     psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
+    myMin->gainFactorMode = options->gainFactorMode;
+    myMin->chisqConvergence = options->chisqConvergence;
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
Index: /branches/eam_branches/ipp-20130509/psModules/src/psmodules.h
===================================================================
--- /branches/eam_branches/ipp-20130509/psModules/src/psmodules.h	(revision 35638)
+++ /branches/eam_branches/ipp-20130509/psModules/src/psmodules.h	(revision 35639)
@@ -158,4 +158,5 @@
 #include <pmReadoutFake.h>
 #include <pmPSFEnvelope.h>
+#include <pmThreadTools.h>
 
 #endif