Changeset 10096 for trunk/psphot/src/psphotChoosePSF.c

Timestamp:

Nov 19, 2006, 7:26:08 PM (20 years ago)

Author:

eugene

Message:

working on global linear fit with sky level

File:

• trunk/psphot/src/psphotChoosePSF.c (modified) (5 diffs)

trunk/psphot/src/psphotChoosePSF.c

@@ -26 +26 @@
     }
 
+    // supply the measured sky variance for optional constant errors (non-poissonian)
+    float SKY_STDEV = psMetadataLookupF32 (&status, recipe, "SKY_STDEV");
+    if (!status) {
+        SKY_STDEV = 1.0;
+        psWarning("SKY_STDEV is not set --- defaulting to %f\n", SKY_STDEV);
+    }
     // use poissonian errors or local-sky errors
     bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
@@ -32 +38 @@
         psWarning("POISSON_ERRORS is not set in the recipe --- defaulting to true.\n");
     }
+    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_STDEV), POISSON_ERRORS);
 
     // how to model the PSF variations across the field
@@ -48 +55 @@
     }
 
-    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
-
     stars = psArrayAllocEmpty (sources->n);
 
@@ -94 +99 @@
         modelName = item->data.V;
         models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask);
-    }
-
-    // XXX test dump of psf stars and model
-    if (1) { 
-        psphotSaveImage (NULL, readout->image,  "testsub.fits");
-        pmSourcesWritePSFs (stars, "psfstars.dat");
-        try = models->data[0];
-        psf = try->psf;
-        psMetadata *psfData = pmPSFtoMetadata (NULL, psf);
-        psMetadataConfigWrite (psfData, "psfmodel.dat");
-        psFree (psfData);
     }
 
@@ -138 +132 @@
     try = models->data[bestN];
 
+    // XXX test dump of psf star data and psf-subtracted image
+    if (0) { 
+        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
+            psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_MASK_MARK);
+            pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
+            psImageKeepCircle (source->mask, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
+        }
+
+        psphotSaveImage (NULL, readout->image,  "psfstars.fits");
+        pmSourcesWritePSFs (try->sources, "psfstars.dat");
+        psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
+        psMetadataConfigWrite (psfData, "psfmodel.dat");
+        psFree (psfData);
+        psLogMsg ("psphot.choosePSF", 3, "wrote out psf-subtracted image, psf data, exiting\n");
+        exit (0);
+    }
+
     // unset the PSFSTAR flag for stars not used for PSF model
     for (int i = 0; i < try->sources->n; i++) {