Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 20081)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 20082)
@@ -43,7 +43,4 @@
     // array to store candidate PSF stars
     int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS");
-    assert (status);
-
-    float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
     assert (status);
 
@@ -258,65 +255,4 @@
     pmPSFtry *try = models->data[bestN];
 
-    // measure and save the median value of dparams[PM_PAR_SXX],dparams[PM_PAR_SYY]
-    // these are used by psphotEvalPSF to identify extended sources
-    // XXX is this still needed?
-    psVector *Sx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
-    psVector *Sy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
-    psVector *dSx = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
-    psVector *dSy = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
-    psVector *dSN = psVectorAllocEmpty (try->sources->n, PS_TYPE_F32);
-    for (int i = 0; i < try->sources->n; i++) {
-        // masked for: bad model fit, outlier in parameters
-        if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
-            continue;
-
-        pmSource *source = try->sources->data[i];
-        Sx->data.F32[Sx->n] = source->modelPSF->params->data.F32[PM_PAR_SXX];
-        Sy->data.F32[Sy->n] = source->modelPSF->params->data.F32[PM_PAR_SYY];
-        dSN->data.F32[dSN->n] = source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0];
-        dSx->data.F32[dSx->n] = source->modelPSF->dparams->data.F32[PM_PAR_SXX];
-        dSy->data.F32[dSy->n] = source->modelPSF->dparams->data.F32[PM_PAR_SYY];
-        Sx->n ++;
-        Sy->n ++;
-        dSN->n ++;
-        dSx->n ++;
-        dSy->n ++;
-    }
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-
-    psVectorStats (stats, dSx, NULL, NULL, 0);
-    float dSxo = stats->sampleMean;
-    psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "dSxo: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
-
-    psVectorStats (stats, dSy, NULL, NULL, 0);
-    float dSyo = stats->sampleMean;
-    psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "dSyo: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
-
-    for (int i = 0; i < dSN->n; i++) {
-        dSx->data.F32[i] = (dSx->data.F32[i] - dSxo) / PS_MAX (PSF_MIN_DS, Sx->data.F32[i]*dSN->data.F32[i]);
-        dSy->data.F32[i] = (dSy->data.F32[i] - dSyo) / PS_MAX (PSF_MIN_DS, Sy->data.F32[i]*dSN->data.F32[i]);
-    }
-
-    psVectorStats (stats, dSx, NULL, NULL, 0);
-    float nSx = stats->sampleStdev;
-    psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "ndSx: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSx->n);
-    psVectorStats (stats, dSy, NULL, NULL, 0);
-    float nSy = stats->sampleStdev;
-    psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "ndSy: mean: %f, median: %f, stdev: %f, npts: %ld", stats->sampleMean, stats->sampleMedian, stats->sampleStdev, dSy->n);
-
-    psFree (Sx);
-    psFree (Sy);
-    psFree (dSx);
-    psFree (dSy);
-    psFree (dSN);
-    psFree (stats);
-
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_MEAN", PS_META_REPLACE, "mean offset of dSXX", dSxo);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_MEAN", PS_META_REPLACE, "mean offset of dSYY", dSyo);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSX_STDV", PS_META_REPLACE, "stdev of mean-corrected dSXX", nSx);
-    psMetadataAddF32 (recipe, PS_LIST_TAIL, "DSY_STDV", PS_META_REPLACE, "stdev of mean-corrected dSYY", nSy);
-
-    // psphotCountPSFStars (sources);
-
     // unset the PSFSTAR flag for stars not used for PSF model
     for (int i = 0; i < try->sources->n; i++) {
@@ -326,6 +262,4 @@
         }
     }
-
-    // psphotCountPSFStars (sources);
 
     // build a PSF residual image
@@ -337,5 +271,5 @@
     }
 
-    // build a PSF residual image
+    // build the flux-to-magnitude conversion information
     if (!psphotMakeFluxScale (readout->image, recipe, try->psf)) {
         psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF");
@@ -345,72 +279,15 @@
     }
 
+    // build curve-of-growth vector for this psf
+    if (!psphotMakeGrowthCurve (readout, recipe, try->psf)) {
+        psError(PSPHOT_ERR_PSF, false, "Unable to construct flux scale for PSF");
+        psFree (models);
+        psFree(options);
+        return NULL;
+    }
+
     // XXX test dump of psf star data and psf-subtracted image
     if (psTraceGetLevel("psphot.psfstars") > 5) {
-        psphotSaveImage (NULL, readout->image,  "rawstars.fits");
-
-        for (int i = 0; i < try->sources->n; i++) {
-            // masked for: bad model fit, outlier in parameters
-            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL)
-                continue;
-
-            pmSource *source = try->sources->data[i];
-            float x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-            float y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-
-            // set the mask and subtract the PSF model
-            // XXX should we be using maskObj? should we be unsetting the mask?
-            // use pmModelSub because modelFlux has not been generated
-            assert (source->maskObj);
-            psImageKeepCircle (source->maskObj, x, y, options->radius, "OR", markVal);
-            pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
-            psImageKeepCircle (source->maskObj, x, y, options->radius, "AND", PS_NOT_U8(markVal));
-        }
-
-        FILE *f = fopen ("shapes.dat", "w");
-        for (int i = 0; i < try->sources->n; i++) {
-            psF32 inPar[10];  // must be psF32 to pmPSF_FitToModel
-
-            // masked for: bad model fit, outlier in parameters
-            if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-
-            pmSource *source = try->sources->data[i];
-            psF32 *outPar = source->modelEXT->params->data.F32;
-
-            psEllipseShape shape;
-
-            shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
-            shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
-            shape.sxy = outPar[PM_PAR_SXY];
-
-            psEllipsePol pol = pmPSF_ModelToFit (outPar);
-            inPar[PM_PAR_E0] = pol.e0;
-            inPar[PM_PAR_E1] = pol.e1;
-            inPar[PM_PAR_E2] = pol.e2;
-            pmPSF_FitToModel (inPar, 0.1);
-
-            psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
-            psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1);
-
-            psEllipsePol pol2 = psEllipseAxesToPol (axes1);
-
-            fprintf (f, "%3d  %7.2f %7.2f  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f  :  %7.4f %7.4f %7.4f  --  %7.4f %7.4f %7.4f : %7.4f %7.4f %6.1f : %7.4f %7.4f %6.1f\n",
-                     i, outPar[PM_PAR_XPOS], outPar[PM_PAR_YPOS],
-                     outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY],
-                     pol.e0, pol.e1, pol.e2,
-                     pol2.e0, pol2.e1, pol2.e2,
-                     inPar[PM_PAR_SXX], inPar[PM_PAR_SXY], inPar[PM_PAR_SYY],
-                     axes1.major, axes1.minor, axes1.theta*PM_DEG_RAD,
-                     axes2.major, axes2.minor, axes2.theta*PM_DEG_RAD
-                     );
-        }
-        fclose (f);
-
-        psphotSaveImage (NULL, readout->image,  "psfstars.fits");
-        pmSourcesWritePSFs (try->sources, "psfstars.dat");
-        pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
-        // XXX need alternative output function
-        // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
-        // psMetadataConfigWrite (psfData, "psfmodel.dat");
-        psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
+	psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal);
     }