Index: trunk/psphot/src/psphotOutput.c
===================================================================
--- trunk/psphot/src/psphotOutput.c	(revision 20082)
+++ trunk/psphot/src/psphotOutput.c	(revision 20083)
@@ -226,2 +226,75 @@
     return header;
 }
+
+// XXX add args as needed
+bool psphotDumpPSFStars (pmReadout *readout, pmPSFtry *try, float radius, psMaskType maskVal, psMaskType markVal) {
+
+    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, radius, "OR", markVal);
+	pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
+	psImageKeepCircle (source->maskObj, x, y, 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");
+
+    return true;
+}