Changeset 34083 for trunk/ppSim/src/ppSimInsertStars.c

Timestamp:

Jun 26, 2012, 11:25:02 AM (14 years ago)

Author:

eugene

Message:

fix numerical round-off errors and aperture bias in ppSimInsertStars

File:

• trunk/ppSim/src/ppSimInsertStars.c (modified) (5 diffs)

trunk/ppSim/src/ppSimInsertStars.c

@@ -6 +6 @@
     PTR = NULL;
 
-float imageSum (psImage *image);
+double imageSum (psImage *image);
 
 bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
@@ -79 +79 @@
     FILE *outfile = fopen (outname, "w");
 
+    float radius = -1.0;
+
     // add sources to the readout image & variance
     psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
@@ -110 +112 @@
         pmModelSetFlux (model, flux);
 
-        float radius = model->modelRadius (model->params, nSigmaLim * roughNoise);
-        radius = PS_MAX (radius, 1.0);
-        radius = 20.0;
+        // we must set the radius on the first object, and leave it there for everyone else
+        // (inserting flux with a variable radius introduces a flux-dependent bias)
+        if (radius < 0) {
+          radius = model->modelRadius (model->params, nSigmaLim * roughNoise);
+          radius = PS_MAX (radius,  1.0);
+          radius = PS_MIN (radius, 50.0);
+        }
 
         // construct a source, with model flux pixels set, based on the model
@@ -136 +142 @@
             readout->image->data.F32[(int)(yCell)][(int)(xCell)] += flux;
         } else {
-            // insert the source flux in the image
-            float sum1 = imageSum(source->pixels);
+            // insert the source flux in the image and measure inserted flux
+            // this is a difference of 2 large numbers : must use 'double' or we get floating
+            // point errors and trends on a scale of ~0.005 mags at > -10 mags instrumental
+            double sum1 = imageSum(source->pixels);
             pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
-            float sum2 = imageSum(source->pixels);
+            double sum2 = imageSum(source->pixels);
             starFlux = sum2 - sum1;
             
@@ -154 +162 @@
             RESET(source->blends);
         }
-        fprintf (outfile, "%8.3f %8.3f %10.2f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f  :  %f %f\n", star->x, star->y, starFlux, 0, source->psfMag, source->psfMagErr, axes.major, axes.minor, axes.theta, par8, model->params->data.F32[PM_PAR_I0], radius);
+        fprintf (outfile, "%8.3f %8.3f %10.5f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f  :  %f %f\n", star->x, star->y, starFlux, 0, source->psfMag, source->psfMagErr, axes.major, axes.minor, axes.theta, par8, model->params->data.F32[PM_PAR_I0], radius);
 
         // add the sources to the source array