Index: branches/tap_branches/ppSub/src/ppSubMakePSF.c
===================================================================
--- branches/tap_branches/ppSub/src/ppSubMakePSF.c	(revision 25900)
+++ branches/tap_branches/ppSub/src/ppSubMakePSF.c	(revision 27838)
@@ -21,4 +21,6 @@
 
 #include "ppSub.h"
+
+psArray *ppSubSelectPSFSources(psArray *sources);
 
 bool ppSubMakePSF(ppSubData *data)
@@ -57,17 +59,39 @@
     pmFPAfile *photFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.INPUT"); // Photometry file
     if (!pmFPACopy(photFile->fpa, minuendFile->fpa)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to copy FPA for photometry");
+        psError(PPSUB_ERR_CONFIG, false, "Unable to copy FPA for photometry");
         psFree(view);
         return false;
     }
+
     pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout to photometer
-    if (psMetadataLookup(photRO->analysis, "PSPHOT.SOURCES")) {
-        psMetadataRemoveKey(photRO->analysis, "PSPHOT.SOURCES");
+
+    // we need to remove any existing PSPHOT.DETECTIONS (why not do this in psphot?)
+    if (psMetadataLookup(photRO->analysis, "PSPHOT.DETECTIONS")) {
+        psMetadataRemoveKey(photRO->analysis, "PSPHOT.DETECTIONS");
     }
+    if (psMetadataLookup(photRO->parent->parent->analysis, "PSPHOT.PSF")) {
+        psMetadataRemoveKey(photRO->parent->parent->analysis, "PSPHOT.PSF");
+    }
+
+# ifdef TESTING
+    // XXX for testing, dump these images:
+    psphotSaveImage (NULL, photRO->image, "findpsf.im.fits");
+    psphotSaveImage (NULL, photRO->variance, "findpsf.wt.fits");
+    psphotSaveImage (NULL, photRO->mask, "findpsf.mk.fits");
+# endif
 
     // Extract the loaded sources from the associated readout, and generate PSF
     // Here, we assume the image is background-subtracted
-    psArray *sources = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.SOURCES");
-    if (!psphotReadoutFindPSF(config, view, sources)) {
+    pmDetections *detections = psMetadataLookupPtr(&mdok, minuend->analysis, "PSPHOT.DETECTIONS");
+    if (!detections || !detections->allSources) {
+        psError(PPSUB_ERR_CONFIG, true, "No sources from which to determine PSF.");
+        return false;
+    }
+    psArray *sources = detections->allSources;
+
+    // XXX filter sources?  limit the total number and return only brighter objects?
+    // use flags to toss totally bogus entries?
+    psArray *goodSources = ppSubSelectPSFSources (sources);
+    if (!psphotReadoutFindPSF(config, view, goodSources)) {
         // This is likely a data quality issue
         // XXX Split into multiple cases using error codes?
@@ -76,5 +100,14 @@
         ppSubDataQuality(data, PSPHOT_ERR_PSF, PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
         psFree(view);
+        psFree(goodSources);
         return true;
+    }
+
+    // save the resulting PSF information on the pmFPAfile PSPHOT.PSF.LOAD
+    pmFPAfile *psfFile = psMetadataLookupPtr(&mdok, config->files, "PSPHOT.PSF.LOAD"); // PSF file
+    if (!ppSubCopyPSF(psfFile, photFile, view)) {
+        psErrorStackPrint(stderr, "PSF was not generated");
+        psWarning("PSF was not generated --- suspect bad data quality.");
+        ppSubDataQuality(data, psErrorCodeLast(), PPSUB_FILES_PHOT_SUB | PPSUB_FILES_PHOT_INV);
     }
 
@@ -82,11 +115,70 @@
     psMetadataRemoveKey(photRO->analysis, "PSPHOT.HEADER");
 
-    if (!data->quality) {
-        data->psf = psMemIncrRefCounter(psMetadataLookupPtr(NULL, photRO->parent->parent->analysis,
-                                                            "PSPHOT.PSF"));
-    }
-
+    psFree(goodSources);
     psFree(view);
 
     return true;
 }
+
+bool ppSubCopyPSF(pmFPAfile *output, pmFPAfile *input, pmFPAview *view)
+{
+    pmChip *inputChip   = pmFPAviewThisChip(view, input->fpa); // Chip with PSF info
+    pmChip *outputChip  = pmFPAviewThisChip(view, output->fpa); // Chip to store PSF info
+
+    pmReadout *inputRO  = pmFPAviewThisReadout(view, input->fpa); // Readout with PSF info
+    pmReadout *outputRO = pmFPAviewThisReadout(view, output->fpa); // Readout to store PSF info
+
+    if (!outputRO) {
+        pmCell *outputCell  = pmFPAviewThisCell(view, output->fpa);
+        outputRO = pmReadoutAlloc(outputCell);
+        outputRO->image = psMemIncrRefCounter(inputRO->image);
+    }
+
+    // Copy the PSF-related data
+    psMetadataIterator *iter = psMetadataIteratorAlloc(inputRO->analysis, PS_LIST_HEAD, "^PSF\\.CLUMP.*");
+    psMetadataItem *item;               // Item from iteration
+    while ((item = psMetadataGetAndIncrement(iter))) {
+        psMetadataAddItem(outputRO->analysis, item, PS_LIST_TAIL, PS_META_REPLACE);
+    }
+    psFree(iter);
+
+    // copy the PSF model data
+    pmPSF *psf = psMetadataLookupPtr(NULL, inputChip->analysis, "PSPHOT.PSF"); // PSF for photometry
+    if (!psf) {
+        psErrorStackPrint(stderr, "No PSF available");
+        return false;
+    }
+
+    psMetadataAddPtr(outputChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN | PS_META_REPLACE,
+                     "PSF from ppSubMakePSF", psf);
+
+    return true;
+}
+
+
+// XXX hardwired MAX for now
+# define MAX_NPSF 500
+
+psArray *ppSubSelectPSFSources(psArray *sources){
+
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    psArray *subset = psArrayAllocEmpty(MAX_NPSF);
+
+    int nPSF = 0;
+    for (int i = 0; (nPSF < MAX_NPSF) && (i < sources->n); i++) {
+
+        pmSource *source = sources->data[i];
+        if (!source) continue;
+
+        // skip non-astronomical objects (very likely defects)
+        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+
+        psArrayAdd (subset, 100, source);
+        nPSF++;
+    }
+
+    return subset;
+}