Index: trunk/psphot/src/psphotBlendFit.c
===================================================================
--- trunk/psphot/src/psphotBlendFit.c	(revision 20453)
+++ trunk/psphot/src/psphotBlendFit.c	(revision 21166)
@@ -2,5 +2,5 @@
 
 // XXX I don't like this name
-bool psphotBlendFit (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
+bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
 
     int Nfit = 0;
@@ -12,4 +12,124 @@
     psTimerStart ("psphot.fit.nonlinear");
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+	nThreads = 0;
+    }
+    // nThreads = 0; // XXX until testing is complete, do not thread this function
+
+    psphotInitLimitsPSF (recipe, readout);
+    psphotInitLimitsEXT (recipe);
+    psphotInitRadiusPSF (recipe, psf->type);
+
+    // starts the timer, sets up the array of fitSets 
+    psphotFitInit (nThreads);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    fprintf (stderr, "starting with %ld sources\n", sources->n);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+	psArray *cells = cellGroups->data[i];
+
+	for (int j = 0; j < cells->n; j++) {
+
+	    // allocate a job -- if threads are not defined, this just runs the job
+	    psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
+	    psArray *newSources = psArrayAllocEmpty(16);
+
+	    psArrayAdd(job->args, 1, readout);
+	    psArrayAdd(job->args, 1, recipe);
+	    psArrayAdd(job->args, 1, cells->data[j]); // sources
+	    psArrayAdd(job->args, 1, psf);
+	    psArrayAdd(job->args, 1, newSources); // return for new sources
+	    psFree (newSources);
+
+	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
+	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
+	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+	    PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+	    if (!psThreadJobAddPending(job)) {
+		psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+		psFree (job);
+		return NULL;
+	    }
+	    psFree(job);
+
+	}
+
+	// wait for the threads to finish and manage results
+	if (!psThreadPoolWait (false)) {
+	    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+	    return NULL;
+	}
+
+	// we have only supplied one type of job, so we can assume the types here
+	psThreadJob *job = NULL;
+	while ((job = psThreadJobGetDone()) != NULL) {
+	    if (job->args->n < 1) {
+		fprintf (stderr, "error with job\n");
+	    } else {
+		psScalar *scalar = NULL;
+		scalar = job->args->data[5];
+		Nfit += scalar->data.S32;
+		scalar = job->args->data[6];
+		Npsf += scalar->data.S32;
+		scalar = job->args->data[7];
+		Next += scalar->data.S32;
+		scalar = job->args->data[8];
+		Nfail += scalar->data.S32;
+
+		// add these back onto sources
+		psArray *newSources = job->args->data[4];
+		for (int j = 0; j < newSources->n; j++) {
+		    psArrayAdd (sources, 16, newSources->data[j]);
+		}
+	    }
+	    psFree(job);
+	}
+    }
+    psFree (cellGroups);
+
+    if (psTraceGetLevel("psphot") >= 6) {
+      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
+    }
+
+    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
+
+    psphotVisualShowResidualImage (readout);
+    psphotVisualShowFlags (sources);
+
+    return true;
+}
+
+bool psphotBlendFit_Threaded (psThreadJob *job) {
+
+    bool status = false;
+    int Nfit = 0;
+    int Npsf = 0;
+    int Next = 0;
+    int Nfail = 0;
+    psScalar *scalar = NULL;
+
+    pmReadout *readout  = job->args->data[0];
+    psMetadata *recipe  = job->args->data[1];
+    psArray *sources    = job->args->data[2];
+    pmPSF *psf          = job->args->data[3];
+    psArray *newSources = job->args->data[4];
+
     // bit-masks to test for good/bad pixels
     psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT");
@@ -23,15 +143,6 @@
     maskVal |= markVal;
 
-    // source analysis is done in S/N order (brightest first)
-    sources = psArraySort (sources, pmSourceSortBySN);
-
     // S/N limit to perform full non-linear fits
     float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
-
-    psphotInitLimitsPSF (recipe, readout);
-    psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, psf->type);
-
-    psphotFitInit ();
 
     // option to limit analysis to a specific region
@@ -41,6 +152,4 @@
 
     for (int i = 0; i < sources->n; i++) {
-        // if (i%100 == 0) psphotFitSummary ();
-
         pmSource *source = sources->data[i];
 
@@ -57,5 +166,5 @@
         if (source->peak->SN < FIT_SN_LIM) continue;
 
-        // XXX this should use peak?
+	// exclude sources outside optional analysis region
         if (source->peak->xf < AnalysisRegion.x0) continue;
         if (source->peak->yf < AnalysisRegion.y0) continue;
@@ -66,5 +175,6 @@
         if (source->modelPSF == NULL) continue;
 
-        // skip sources which are insignificant flux?
+        // skip sources which are insignificant flux? 
+	// XXX this is somewhat ad-hoc
         if (source->modelPSF->params->data.F32[1] < 0.1) {
             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
@@ -81,15 +191,8 @@
         Nfit ++;
 
-	// XXX TEST
-	// if ((fabs(source->peak->x - 1202) < 2) && (fabs(source->peak->y - 1065) < 2)) {
-	// psTraceSetLevel ("psLib.math.psMinimizeLMChi2", 6);
-	// }
-
-        // try fitting PSFs, then try extended sources
-        // these functions subtract the resulting fitted source (XXX and update the modelFlux?)
-	// XXX re-consider conditions under which the source has EXT fit:
-	// I should try EXT if the source size measurement says it is large
+        // try fitting PSFs or extended sources depending on source->mode
+        // these functions subtract the resulting fitted source
 	if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-	    if (psphotFitBlob (readout, source, sources, psf, maskVal, markVal)) {
+	    if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
 		source->type = PM_SOURCE_TYPE_EXTENDED;
 		psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
@@ -115,12 +218,16 @@
     }
 
-    if (psTraceGetLevel("psphot") >= 6) {
-      psphotSaveImage (NULL, readout->image,  "image.v2.fits");
-    }
-
-    psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
-
-    psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
+    // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[5];
+    scalar->data.S32 = Nfit;
+
+    scalar = job->args->data[6];
+    scalar->data.S32 = Npsf;
+
+    scalar = job->args->data[7];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[8];
+    scalar->data.S32 = Nfail;
 
     return true;