IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 14816


Ignore:
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.

Location:
trunk/ppSim/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSim/src/ppSimInsertStars.c

    r14668 r14816  
    11# include "ppSim.h"
     2
     3// Reset a pointer: free and set to NULL
     4#define RESET(PTR) \
     5    psFree(PTR); \
     6    PTR = NULL;
     7
    28
    39bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config) {
    … …  
    2632    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    2733    if (isnan(readnoise)) {
    28         psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
    29         readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
    30         if (!mdok) {
    31             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    32             return false;
    33         }
     34        psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
     35        readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
     36        if (!mdok) {
     37            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
     38            return false;
     39        }
    3440    }
    3541
    … …  
    5864
    5965    // add sources to the readout image & weight
     66    psTrace("ppSim", 1, "Inserting %ld stars...\n", stars->n);
    6067    for (long i = 0; i < stars->n; i++) {
    61         ppSimStar *star = stars->data[i];
     68        ppSimStar *star = stars->data[i];
     69        psTrace("ppSim", 10, "Inserting star at %.1f,%.1f --> %.2f\n", star->x, star->y, star->flux);
    6270
    63         // star->x,y are in fpa coordinates
     71        // star->x,y are in fpa coordinates
    6472
    65         // Position on the cell and peak flux
    66         float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
    67         float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
     73        // Position on the cell and peak flux
     74        float xChip = PM_FPA_TO_CHIP(star->x, x0Chip, xParityChip);
     75        float yChip = PM_FPA_TO_CHIP(star->y, y0Chip, yParityChip);
    6876
    69         // Position on the cell and peak flux
    70         float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
    71         float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
     77        // Position on the cell and peak flux
     78        float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
     79        float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
    7280
    73         if (xCell < 0) continue;
    74         if (yCell < 0) continue;
    75         if (xCell > readout->image->numCols) continue;
    76         if (yCell > readout->image->numRows) continue;
    77         // XXX need to apply col0, row0 if readout is a subarray
     81        // XXX Note, the below does not put the edges of stars on the readout if they fall slightly off
     82        // This will be visible as cut-off stars at amplifier boundaries (e.g., Megacam)
     83        if (xCell < 0) continue;
     84        if (yCell < 0) continue;
     85        if (xCell > readout->image->numCols) continue;
     86        if (yCell > readout->image->numRows) continue;
     87        // XXX need to apply col0, row0 if readout is a subarray
    7888
    79         // XXX apply the expCorr to the star->flux before setting the model flux
     89        // Apply the expCorr to the star->flux before setting the model flux
     90        float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];
    8091
    81         // instantiate a model for the PSF at this location, set desired flux
    82         pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
    83         pmModelSetFlux (model, star->flux);
     92        // instantiate a model for the PSF at this location, set desired flux
     93        pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
     94        pmModelSetFlux (model, flux);
    8495
    85         // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
    86         float radius = model->modelRadius (model->params, roughNoise);
    87         radius = PS_MAX (radius, 1.0);
     96        // XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
     97        float radius = model->modelRadius (model->params, roughNoise);
     98        radius = PS_MAX (radius, 1.0);
    8899
    89         // construct a source, with model flux pixels set, based on the model
    90         pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
     100        // construct a source, with model flux pixels set, based on the model
     101        pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
    91102
    92         // XXX set the mag & err values (should this be done in pmSourceFromModel?)
    93         // XXX i should be applying the gain and the correct effective area
    94         psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
    95         psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
     103        // XXX set the mag & err values (should this be done in pmSourceFromModel?)
     104        // XXX i should be applying the gain and the correct effective area
     105        psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
     106        psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
    96107
    97         source->psfMag = -2.5*log10(star->flux);
    98         source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux;       
    99        
    100         // XXX add the sources to a source array
     108        source->psfMag = -2.5*log10(flux);
     109        source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
    101110
    102         // insert the source flux in the image
    103         pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
    104         psArrayAdd (sources, 100,source);
     111        // XXX add the sources to a source array
     112
     113        // insert the source flux in the image
     114        pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
     115        psArrayAdd (sources, 100,source);
     116        psFree(source);                 // Drop reference
     117
     118        // Blow away the image parts of the source, which makes the memory explode
     119        RESET(source->pixels);
     120        RESET(source->weight);
     121        RESET(source->maskObj);
     122        RESET(source->maskView);
     123        RESET(source->modelFlux);
     124        RESET(source->psfFlux);
     125        RESET(source->blends);
    105126    }
    106127
    107128    // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
    108129    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
     130    psFree(sources);
    109131
    110     // XXX many leaks in here, i think 
     132    // XXX many leaks in here, i think
    111133    return true;
    112134}

  • trunk/ppSim/src/ppSimMakeSky.c

    r14657 r14816  
    88
    99    pmCell *cell = readout->parent;
    10     pmChip *chip = cell->parent; 
     10    pmChip *chip = cell->parent;
    1111    pmFPA  *fpa  = chip->parent;
    1212
    … …  
    3838    for (int y = 0; y < signal->numRows; y++) {
    3939
    40         float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
    41             (bounds->y1 - bounds->y0); // Relative y position in FPA
     40        float yFPA = PPSIM_CELL_TO_FPA(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
     41            (bounds->y1 - bounds->y0); // Relative y position in FPA
    4242
    43         for (int x = 0; x < signal->numCols; x++) {
    44             float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
    45                 (bounds->x1 - bounds->x0); // Relative x position in FPA
     43        for (int x = 0; x < signal->numCols; x++) {
     44            float xFPA = PPSIM_CELL_TO_FPA(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
     45                (bounds->x1 - bounds->x0); // Relative x position in FPA
    4646
    47             // Shutter: adjust exposure time
    48             float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
     47            // Shutter: adjust exposure time
     48            float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
    4949
    50             // Gaussian flat-field over the FPA
    51             float flatValue = exp(-0.5 / PS_SQR(flatSigma) *
    52                                   (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrt(M_PI);
     50            // Gaussian flat-field over the FPA
     51            float flatValue = expf(-0.5 / PS_SQR(flatSigma) * (PS_SQR(yFPA) + PS_SQR(xFPA))) /
     52                flatSigma / sqrtf(2.0 * M_PI);
    5353
    54             if (type == PPSIM_TYPE_FLAT) {
    55                 float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field
    56                 signal->data.F32[y][x] += flatFlux;
    57                 variance->data.F32[y][x] += flatFlux;
    58                 continue;
    59             }
     54            if (type == PPSIM_TYPE_FLAT) {
     55                float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field
     56                signal->data.F32[y][x] += flatFlux;
     57                variance->data.F32[y][x] += flatFlux;
     58                continue;
     59            }
    6060
    61             expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
     61            expCorr->data.F32[y][x] = realExpTime / expTime;
    6262
    63             // Sky background
    64             float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky
    65             signal->data.F32[y][x] += skyFlux;
    66             variance->data.F32[y][x] += skyFlux;
     63            // Sky background
     64            float skyFlux = skyRate * flatValue * realExpTime; // Flux from sky
     65            signal->data.F32[y][x] += skyFlux;
     66            variance->data.F32[y][x] += skyFlux;
    6767
    68             // TO DO: Add fringes
     68            // TO DO: Add fringes
    6969
    70         }
     70        }
    7171    }
    7272    psFree(bounds);

  • trunk/ppSim/src/ppSimMakeStars.c

    r14667 r14816  
    11# include "ppSim.h"
     2
     3#define FAINT_FUDGE_FACTOR   1.0        // Fraction of the noise in a (boxcar) PSF for the faint limit
     4
    25
    36bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
    … …  
    3639    float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
    3740    if (!mdok) {
    38         psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    39         psFree(bounds);
    40         return false;
     41        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
     42        psFree(bounds);
     43        return false;
    4144    }
    4245
    43     // Peak fluxes: faintest and brightest levels for random stars
    44     float faint = 0.1 * 10.0 * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
    45     float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
     46    // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars
     47    float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) *
     48        seeing; // Faint limit is related to the noise in a (boxcar) PSF
     49    float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag
     50    psTrace("ppSim", 6, "Faint limit: %f\n", faint);
     51    psTrace("ppSim", 6, "Bright limit: %f\n", bright);
    4652    if (bright < faint) {
    47         psLogMsg("ppSim", PS_LOG_INFO,
    48                 "Image noise is above brightest random star --- no random stars added.");
    49         return true;
     53        psLogMsg("ppSim", PS_LOG_INFO,
     54                "Image noise is above brightest random star --- no random stars added.");
     55        return true;
    5056    }
    5157
    5258    // Normalisation, set by the specified stellar density at the specified bright magnitude
    5359    float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
    54         powf(bright, starsLum);
     60        powf(bright, starsLum);
    5561
    5662    // Total number of stars down to the faint flux end
    5763    long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
    5864
    59     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
    60              -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
     65    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n",
     66             num, -2.5 * log10(faint / expTime) + zp);
    6167
    6268    long oldSize = stars->n;
    … …  
    6470
    6571    for (long i = 0; i < num; i++) {
    66         ppSimStar *star = ppSimStarAlloc ();
     72        ppSimStar *star = ppSimStarAlloc ();
    6773
    68         // make fpa center of distribution
    69         star->x    = psRandomUniform(rng) * xSize; // x position
    70         star->y    = psRandomUniform(rng) * ySize; // y position
     74        // make fpa center of distribution
     75        star->x    = psRandomUniform(rng) * xSize; // x position
     76        star->y    = psRandomUniform(rng) * ySize; // y position
    7177
    72         star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
    73         star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
     78        star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
     79        star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
    7480
    75         stars->data[oldSize + i] = star;
     81        stars->data[oldSize + i] = star;
    7682    }
    77     stars->n = stars->nalloc;
     83    stars->n = oldSize + num;
    7884
    7985    return true;
Note: See TracChangeset for help on using the changeset viewer.