Index: trunk/pswarp/src/pswarpLoop.c
===================================================================
--- trunk/pswarp/src/pswarpLoop.c	(revision 20115)
+++ trunk/pswarp/src/pswarpLoop.c	(revision 20348)
@@ -3,4 +3,6 @@
 
 #define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
+#define PSPHOT_FIND_PSF 0               // Use psphot's findPSF function?
+//#define TESTING 1                       // Testing output?
 
 // Lists of file rules for the detectors and skycells, and an independent list of everything else
@@ -8,5 +10,9 @@
 static char *detectorFiles[] = { "PSWARP.INPUT", "PSWARP.MASK", "PSWARP.WEIGHT", NULL };
 static char *skycellFiles[] = { "PSWARP.OUTPUT", "PSWARP.OUTPUT.MASK", "PSWARP.OUTPUT.WEIGHT", NULL };
-static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.INPUT.CMF", "PSPHOT.RESID",
+static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT",
+#if PSPHOT_FIND_PSF
+                             "PSPHOT.INPUT.CMF",
+#endif
+                             "PSPHOT.RESID",
                              "PSPHOT.BACKMDL", "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB",
                              "PSPHOT.PSF.SAVE", "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL",
@@ -350,5 +356,5 @@
                           "Don't allow psphot to fit a fake PSF!", false);
 
-#if 0
+#if !PSPHOT_FIND_PSF
         // Get rid of the transformed sources so that we can measure the PSF
         psMetadataRemoveKey(output->analysis, "PSPHOT.SOURCES");
@@ -360,8 +366,42 @@
         pmFPAview *view = pmFPAviewAlloc(0); // View into skycell
         view->chip = view->cell = view->readout = 0;
+
+#if PSPHOT_FIND_PSF
+        // Need to adjust the weight --- the main operation in psphotReadoutFindPSF is not determining the
+        // signficance of sources, but looking at the significance of individual pixels.
+        // We can adjust the weight directly since this is a deep copy
+        pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa);
+        float vf = psMetadataLookupF32(NULL, photRO->parent->concepts, "CELL.VARFACTOR"); // Variance factor
+        psBinaryOp(photRO->weight, photRO->weight, "*", psScalarAlloc(vf, PS_TYPE_F32));
+
         if (!psphotReadoutFindPSF(config, view)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF for warped image.");
             return false;
         }
+#else
+        if (!psphotReadout(config, view)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF for warped image.");
+            return false;
+        }
+#endif
+
+#ifdef TESTING
+        {
+            #define PSF_SIZE 20         // Half-size of PSF
+            #define PSF_FLUX 10000      // Central flux for PSF
+            pmChip *photChip = pmFPAviewThisChip(view, photFile->fpa);
+            pmPSF *psf = psMetadataLookupPtr(NULL, photChip->analysis, "PSPHOT.PSF");
+            psImage *image = psImageAlloc(2 * PSF_SIZE + 1, 2 * PSF_SIZE + 1, PS_TYPE_F32);
+            psImageInit(image, 0);
+            pmModel *model = pmModelFromPSFforXY(psf, PSF_SIZE, PSF_SIZE, PSF_FLUX);
+            pmModelAdd(image, NULL, model, PM_MODEL_OP_FULL, 0);
+            psFree(model);
+            psFits *fits = psFitsOpen("psf.fits", "w");
+            psFitsWriteImage(fits, NULL, image, 0, NULL);
+            psFitsClose(fits);
+            psFree(image);
+        }
+#endif
+
         psFree(view);
     }