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

Timestamp:

Sep 11, 2007, 1:01:03 PM (19 years ago)

Author:

Paul Price

Message:

Ensuring consistency between flux scales (peak vs integrated) --- stars now appear on the output image! Some memory fixing too.

File:

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

trunk/ppSim/src/ppSimInsertStars.c

@@ -1 +1 @@
 # include "ppSim.h"
+
+// Reset a pointer: free and set to NULL
+#define RESET(PTR) \
+    psFree(PTR); \
+    PTR = NULL;
+
 
 bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
@@ -26 +32 @@
     float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
     if (isnan(readnoise)) {
-        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
-        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
-        if (!mdok) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
-            return false;
-        }
+        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
+        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
+        if (!mdok) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+            return false;
+        }
     }
 
@@ -58 +64 @@
 
     // add sources to the readout image & weight
+    psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
     for (long i = 0; i < stars->n; i++) {
-        ppSimStar *star = stars->data[i];
+        ppSimStar *star = stars->data[i];
+        psTrace("ppSim", 10, "Inserting star at %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux);
 
-        // star->x,y are in fpa coordinates
+        // star->x,y are in fpa coordinates
 
-        // Position on the cell and peak flux
-        float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
-        float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
+        // Position on the cell and peak flux
+        float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
+        float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
 
-        // Position on the cell and peak flux
-        float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
-        float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
+        // Position on the cell and peak flux
+        float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
+        float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
 
-        if (xCell < 0) continue;
-        if (yCell < 0) continue;
-        if (xCell > readout->image->numCols) continue;
-        if (yCell > readout->image->numRows) continue;
-        // XXX need to apply col0, row0 if readout is a subarray
+        // XXX Note, the below does not put the edges of stars on the readout if they fall slightly off
+        // This will be visible as cut-off stars at amplifier boundaries (e.g., Megacam)
+        if (xCell < 0) continue;
+        if (yCell < 0) continue;
+        if (xCell > readout->image->numCols) continue;
+        if (yCell > readout->image->numRows) continue;
+        // XXX need to apply col0, row0 if readout is a subarray
 
-        // XXX apply the expCorr to the star->flux before setting the model flux
+        // Apply the expCorr to the star->flux before setting the model flux
+        float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];
 
-        // instantiate a model for the PSF at this location, set desired flux
-        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
-        pmModelSetFlux (model, star->flux);
+        // instantiate a model for the PSF at this location, set desired flux
+        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
+        pmModelSetFlux (model, flux);
 
-        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
-        float radius = model->modelRadius (model->params, roughNoise);
-        radius = PS_MAX (radius, 1.0);
+        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
+        float radius = model->modelRadius (model->params, roughNoise);
+        radius = PS_MAX (radius, 1.0);
 
-        // construct a source, with model flux pixels set, based on the model
-        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
+        // construct a source, with model flux pixels set, based on the model
+        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
 
-        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
-        // XXX i should be applying the gain and the correct effective area
-        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
-        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
+        // XXX i should be applying the gain and the correct effective area
+        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 
-        source->psfMag = -2.5*log10(star->flux);
-        source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux;       
-        
-        // XXX add the sources to a source array
+        source->psfMag = -2.5*log10(flux);
+        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
 
-        // insert the source flux in the image
-        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
-        psArrayAdd (sources, 100,source);
+        // XXX add the sources to a source array
+
+        // insert the source flux in the image
+        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+        psArrayAdd (sources, 100,source);
+        psFree(source);                 // Drop reference
+
+        // Blow away the image parts of the source, which makes the memory explode
+        RESET(source->pixels);
+        RESET(source->weight);
+        RESET(source->maskObj);
+        RESET(source->maskView);
+        RESET(source->modelFlux);
+        RESET(source->psfFlux);
+        RESET(source->blends);
     }
 
     // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
     psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
+    psFree(sources);
 
-    // XXX many leaks in here, i think 
+    // XXX many leaks in here, i think
     return true;
 }