IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 9, 2007, 6:03:19 PM (19 years ago)
Author:
eugene
Message:

splitting ppSimLoop into a number of stand-alone functions

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppSim/src/ppSimLoop.c

    r14367 r14463  
    1 #ifdef HAVE_CONFIG_H
    2 #include <config.h>
    3 #endif
    4 
    5 #include <stdio.h>
    6 #include <pslib.h>
    7 #include <psmodules.h>
    8 #include <psastro.h>
    9 
    101#include "ppSim.h"
    11 
    12 #define STAR_RANGE_MIN 4.5              // Minimum pixel range for adding stars, sigma
    13 #define STAR_BG_FRAC 0.1                // Fraction of background error to go out to when adding stars
    14 
    15 #define LOG10(X) (log(X)/log(10))       // Base 10 logarithm
    16 
    17 // Return FPA position, given a cell position; calculations are all done in pixel units
    18 static inline float cell2fpa(float pos, // Cell position
    19                              int cell0, // Cell offset
    20                              int cellParity, // Cell parity
    21                              int binning, // Binning factor
    22                              int chip0, // Chip offset
    23                              int chipParity // Chip parity
    24     )
    25 {
    26     return (chip0 + binning * chipParity * (cell0 + cellParity * pos));
    27 }
    28 
    29 // Return cell position, given an FPA position; calculations are all done in pixel units
    30 static inline float fpa2cell(float pos, // FPA position
    31                              int cell0, // Cell offset
    32                              int cellParity, // Cell parity
    33                              int binning, // Binning factor
    34                              int chip0, // Chip offset
    35                              int chipParity // Chip parity
    36     )
    37 {
    38     return ((pos - chip0) * chipParity - cell0) * cellParity / binning;
    39 }
    40 
    41 // Compare a value with minimum and maximum values, replacing where required.
    42 #define COMPARE(VALUE,MIN,MAX) { \
    43         if (VALUE < MIN) { \
    44             MIN = VALUE; \
    45         } \
    46         if (VALUE > MAX) { \
    47             MAX = VALUE; \
    48         } \
    49     }
    50 
    51 // Return bounds of a chip, based on the concepts
    52 static psRegion *chipBounds(const pmChip *chip, // Chip for which to determine size
    53                             pmFPAview *view // View for chip
    54     )
    55 {
    56     assert(chip);
    57     assert(view);
    58 
    59     // Bounds of chip
    60     int xMin = INT_MAX;
    61     int xMax = - INT_MAX;
    62     int yMin = INT_MAX;
    63     int yMax = - INT_MAX;
    64 
    65     int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
    66     int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
    67     int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY");
    68     int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY");
    69 
    70     if ((xParityChip != 1 && xParityChip != -1) || (yParityChip != 1 && yParityChip != -1)) {
    71         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Chip parities are not set.");
    72         psFree(view);
    73         return NULL;
    74     }
    75 
    76     pmCell *cell;                   // Cell from chip
    77     while ((cell = pmFPAviewNextCell(view, chip->parent, 1))) {
    78         int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
    79         int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");
    80 
    81         if (xSize == 0 || ySize == 0) {
    82             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell sizes are not set.");
    83             psFree(view);
    84             return NULL;
    85         }
    86 
    87         int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    88         int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    89         int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    90         int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    91 
    92         if ((xParityCell != 1 && xParityCell != -1) || (yParityCell != 1 && yParityCell != -1)) {
    93             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell parities are not set.");
    94             psFree(view);
    95             return NULL;
    96         }
    97 
    98         // "Left", "Right", "Bottom" and "Top" don't take into account the parity
    99         int cellLeft = cell2fpa(0, x0Cell, xParityCell, 1, x0Chip, xParityChip);
    100         int cellRight = cell2fpa(xSize, x0Cell, xParityCell, 1, x0Chip, xParityChip);
    101         int cellBottom = cell2fpa(0, y0Cell, yParityCell, 1, y0Chip, yParityChip);
    102         int cellTop = cell2fpa(ySize, y0Cell, yParityCell, 1, y0Chip, yParityChip);
    103 
    104         COMPARE(cellLeft, xMin, xMax);
    105         COMPARE(cellRight, xMin, xMax);
    106         COMPARE(cellBottom, yMin, yMax);
    107         COMPARE(cellTop, yMin, yMax);
    108     }
    109 
    110     return psRegionAlloc(xMin, xMax, yMin, yMax);
    111 }
    112 
    113 // Return bounds of an FPA, based on the concepts
    114 static psRegion *fpaBounds(const pmFPA *fpa       // FPA for which to determine size
    115     )
    116 {
    117     assert(fpa);
    118 
    119     pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA
    120 
    121     // Bounds of focal plane
    122     int xMin = INT_MAX;
    123     int xMax = - INT_MAX;
    124     int yMin = INT_MAX;
    125     int yMax = - INT_MAX;
    126 
    127     pmChip *chip;                       // Chip from FPA
    128     while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
    129         psRegion *bounds = chipBounds(chip, view); // Bounds for chip
    130         COMPARE(bounds->x0, xMin, xMax);
    131         COMPARE(bounds->x1, xMin, xMax);
    132         COMPARE(bounds->y0, yMin, yMax);
    133         COMPARE(bounds->y1, yMin, yMax);
    134         psFree(bounds);
    135     }
    136 
    137     psFree(view);
    138     return psRegionAlloc(xMin, xMax, yMin, yMax);
    139 }
    140 
    141 
    142 // Add noise to an image
    143 static psImage *addNoise(psImage *signal, // Signal image, modified and returned
    144                          const psImage *variance, // Variance image
    145                          psRandom *rng, // Random number generator
    146                          float gain     // CCD gain (e/ADU)
    147     )
    148 {
    149     assert(signal->type.type == PS_TYPE_F32);
    150     assert(signal->numCols == variance->numCols && signal->numRows == variance->numRows);
    151 
    152     // Add the noise into the image
    153     for (int y = 0; y < signal->numRows; y++) {
    154         for (int x = 0; x < signal->numCols; x++) {
    155             signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng);
    156             signal->data.F32[y][x] /= gain; // Converting to ADU
    157         }
    158     }
    159 
    160     return signal;
    161 }
    162 
    163 // Apply saturation limit to image
    164 static void saturate(psImage *image, // Image to apply saturation
    165                      float saturation // Saturation level
    166     )
    167 {
    168     for (int y = 0; y < image->numRows; y++) {
    169         for (int x = 0; x < image->numCols; x++) {
    170             if (image->data.F32[y][x] > saturation) {
    171                 image->data.F32[y][x] = saturation;
    172             }
    173         }
    174     }
    175     return;
    176 }
    177 
    178 // Return star flux at a position
    179 static inline float starFlux(float x, float y, // Position of interest
    180                              float x0, float y0, // Centre of star
    181                              float seeing // Seeing FWHM (pixels)
    182     )
    183 {
    184 #ifdef STAR_GAUSSIAN
    185     // Gaussian star
    186     return exp(-0.5 * (PS_SQR(x - x0) + PS_SQR(y - y0)) / PS_SQR(seeing));
    187 #else
    188     // Waussian star
    189     float z = (PS_SQR(x - x0) + PS_SQR(y - y0)) / (2.0 * PS_SQR(seeing));
    190     return 1.0 / (1.0 + z + PS_SQR(z));
    191 #endif
    192 }
    193 
    194 // Return size of star in pixels
    195 static inline int starSize(float noise, float peak, float seeing)
    196 {
    197 #ifdef STAR_GAUSSIAN
    198     // Gaussian star (solving Gaussian)
    199     float target = STAR_BG_FRAC * seeing * sqrt(M_2_PI) * noise / peak;
    200     return target < 1.0 ? 2.0 * seeing * sqrtf(-logf(target)) + 0.5 : seeing * STAR_RANGE_MIN;
    201 #else
    202     // Waussian star (solving Waussian where z >> 1 and peak/sqrt(bg) >> 1)
    203     float target = STAR_BG_FRAC * noise / peak;
    204     return PS_MAX(sqrtf(1.0 / target) + 0.5, seeing * STAR_RANGE_MIN);
    205 #endif
    206 }
    207 
    208 // Add a star into the signal and variance images
    209 static void star(psImage *signal,       // Signal image, to which to add star
    210                  psImage *variance,     // Variance image, to which to add star
    211                  float x0, float y0,    // Position of star
    212                  float peak,            // Peak flux of star
    213                  float noise,           // Rough noise estimate
    214                  float seeing,          // Seeing for star
    215                  const psImage *correction // Exposure correction as a function of position
    216     )
    217 {
    218     int size = starSize(noise, peak, seeing); // Size of star
    219 
    220     // Range in image pixels on which to add star
    221     int xMin = PS_MAX(0, x0 - size);
    222     int xMax = PS_MIN(signal->numCols - 1, x0 + size);
    223     int yMin = PS_MAX(0, y0 - size);
    224     int yMax = PS_MIN(signal->numRows - 1, y0 + size);
    225 
    226     for (int y = yMin; y <= yMax; y++) {
    227         for (int x = xMin; x <= xMax; x++) {
    228             float star = peak * correction->data.F32[y][x] * starFlux(x, y, x0, y0, seeing); // Star flux
    229             signal->data.F32[y][x] += star;
    230             variance->data.F32[y][x] += star;
    231         }
    232     }
    233     return;
    234 }
    235 
    236 // Generate a header containing WCS keywords
    237 static psMetadata *wcsHeader(float ra0, float dec0, // Boresight (radians)
    238                              float pa,// Position angle (radians)
    239                              float scale, // Pixel scale (radians/pix)
    240                              float x0, float y0, // Position of 0,0 on the FPA
    241                              int xParity, int yParity // Parity in x and y
    242                              )
    243 {
    244     psMetadata *header = psMetadataAlloc(); // Header, to return
    245     pmAstromWCS *wcs = pmAstromWCSAlloc(1, 1); // WCS structure
    246     wcs->toSky = psProjectionAlloc(ra0, dec0, scale * xParity, scale * yParity, PS_PROJ_TAN);
    247     wcs->cdelt1 = scale * PM_DEG_RAD * xParity;
    248     wcs->cdelt2 = scale * PM_DEG_RAD * yParity;
    249     wcs->crpix1 = x0;
    250     wcs->crpix2 = y0;
    251     wcs->trans->x->coeff[1][0] = cos(pa) * wcs->cdelt1;
    252     wcs->trans->x->coeff[0][1] = -sin(pa) * wcs->cdelt1;
    253     wcs->trans->y->coeff[1][0] = sin(pa) * wcs->cdelt2;
    254     wcs->trans->y->coeff[0][1] = cos(pa) * wcs->cdelt2;
    255     // These aren't used by pmAstromWCStoHeader, but set them anyway
    256     wcs->trans->x->coeff[0][0] = ra0;
    257     wcs->trans->y->coeff[0][0] = dec0;
    258     wcs->trans->x->coeff[1][1] = 0.0;
    259     wcs->trans->y->coeff[1][1] = 0.0;
    260 
    261     pmAstromWCStoHeader(header, wcs);
    262     psFree(wcs);
    263 
    264     return header;
    265 }
    266 
    2672
    2683psExit ppSimLoop(pmConfig *config)
    … …  
    27510    assert(fpa);
    27611
    277     bool mdok;                          // Status of MD lookups
     12    ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate
    27813
    279     psRegion *bounds = fpaBounds(fpa);  // Bounds of FPA
    280     if (!bounds || (bounds->x0 == 0.0 && bounds->x1 == 0.0) || (bounds->y0 == 0.0 && bounds->y1 == 0.0)) {
    281         psError(PS_ERR_UNKNOWN, false, "Unable to determine bounds of FPA");
    282         return PS_EXIT_CONFIG_ERROR;
    283     }
     14    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    28415
    285     float biasLevel = psMetadataLookupF32(NULL, config->arguments, "BIAS.LEVEL"); // Bias level
    286     float biasRange = psMetadataLookupF32(NULL, config->arguments, "BIAS.RANGE"); // Bias range
    287     int biasOrder = psMetadataLookupS32(NULL, config->arguments, "BIAS.ORDER"); // Bias order
    288     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
    289     float flatSigma = psMetadataLookupF32(NULL, config->arguments, "FLAT.SIGMA"); // Flat-field coefficient
    290     float flatRate = psMetadataLookupF32(NULL, config->arguments, "FLAT.RATE"); // Flat-field rate
    291     float shutterTime = psMetadataLookupF32(NULL, config->arguments, "SHUTTER.TIME"); // Shutter time
    292     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
    293     float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
    294     float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
    295     float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes
    29616    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
    29717
    298     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
    299     const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name
    300     if (!filter) {
    301         filter = "NONE";
    302     }
    303     ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate
    304 
    305     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
    306     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
    307 
    308     psVector *ra = NULL;                // Star RAs
    309     psVector *dec = NULL;               // Star Declinations
    310     psVector *mag = NULL;               // Star magnitudes
    311     float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
    312     float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") *
    313         M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
    314     float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
    315     float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
    316     float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
    317     float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
    318 
    319     // Add catalogue stars
     18    // Load catalogue stars
     19    psArray *stars = NULL;
    32020    if (type == PPSIM_TYPE_OBJECT) {
    321         // Read catalogue stars using psastro
    322         psMetadata *astroRecipe = psMetadataLookupPtr(NULL, config->recipes, PSASTRO_RECIPE);
    323         if (!astroRecipe) {
    324             psError(PSASTRO_ERR_CONFIG, false, "Can't find recipe %s", PSASTRO_RECIPE);
    325             psFree(rng);
    326             psFree(bounds);
    327             return PS_EXIT_CONFIG_ERROR;
    328         }
    329 
    330         // Size of FPA
    331         float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;
    332 
    333         psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
    334         psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MAX",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 + radius);
    335         psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MIN", PS_DATA_F32 | PS_META_REPLACE, "", dec0 - radius);
    336         psMetadataAdd(astroRecipe, PS_LIST_TAIL, "DEC_MAX", PS_DATA_F32 | PS_META_REPLACE, "", dec0 + radius);
    337         psArray *refStars = psastroLoadRefstars(config);
    338         if (!refStars || refStars->n == 0) {
    339             psError(PS_ERR_UNKNOWN, false, "Unable to find reference stars.");
    340             psFree(refStars);
    341             psFree(rng);
    342             psFree(bounds);
    343             return PS_EXIT_CONFIG_ERROR;
    344         }
    345         psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);
    346 
    347         ra = psVectorAlloc(refStars->n, PS_TYPE_F32);
    348         dec = psVectorAlloc(refStars->n, PS_TYPE_F32);
    349         mag = psVectorAlloc(refStars->n, PS_TYPE_F32);
    350 
    351         // Conversion loop
    352         for (long i = 0; i < refStars->n; i++) {
    353             pmAstromObj *ref = refStars->data[i]; // Reference star
    354             float raStar = ref->sky->r; // RA of star
    355             float decStar = ref->sky->d; // Dec of star
    356             float magStar = ref->Mag;       // Magnitude of star
    357 
    358 
    359             // Convert to x,y position on tangent plane, in pixels
    360             float div = sin(decStar) * sin(dec0) + cos(decStar) * cos(dec0) * cos(raStar - ra0); // Divisor
    361             float xi = cos(decStar) * sin(raStar - ra0) / div / scale;
    362             float eta = (sin(decStar) * cos(dec0) - cos(decStar) * sin(dec0) * cos(raStar - ra0)) /
    363                 div / scale;
    364 
    365             // Apply rotation
    366             ra->data.F32[i] = cos(pa) * xi - sin(pa) * eta;
    367             dec->data.F32[i] = sin(pa) * xi + cos(pa) * eta;
    368 
    369             // Convert magnitude to peak flux
    370             mag->data.F32[i] = powf(10.0, -0.4 * (magStar - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
    371         }
    372         psFree(refStars);
    373 
    374         psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0);
    375         psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0);
     21        stars = ppSimLoadStars (fpa, config);
    37622    }
    37723
    37824    // Add random stars
    379     if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) {
    380 
    381         // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
    382         // shouldn't be too bad.
    383         float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
    384         if (!mdok) {
    385             psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    386             psFree(bounds);
    387             psFree(rng);
    388             return PS_EXIT_CONFIG_ERROR;
    389         }
    390 
    391         // Peak fluxes: faintest and brightest levels for random stars
    392         float faint = 0.1 * 10.0 * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
    393         float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
    394         if (bright < faint) {
    395             psLogMsg("ppSim", PS_LOG_INFO,
    396                      "Image noise is above brightest random star --- no random stars added.");
    397         } else {
    398             // Size of focal plane
    399             int xSize = bounds->x1 - bounds->x0;
    400             int ySize = bounds->y1 - bounds->y0;
    401 
    402             // Normalisation, set by the specified stellar density at the specified bright magnitude
    403             float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
    404                 powf(bright, starsLum);
    405 
    406             // Total number of stars down to the faint flux end
    407             long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
    408 
    409             psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
    410                      -2.5 * LOG10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
    411 
    412             // Extend the vectors
    413             long oldSize;               // Old size of ra, dec, mag vectors
    414             if (ra && dec && mag) {
    415                 oldSize = ra->n;
    416                 ra = psVectorRealloc(ra, oldSize + num);
    417                 dec = psVectorRealloc(dec, oldSize + num);
    418                 mag = psVectorRealloc(mag, oldSize + num);
    419                 ra->n = dec->n = mag->n = oldSize + num;
    420             } else {
    421                 oldSize = 0;
    422                 ra = psVectorAlloc(num, PS_TYPE_F32);
    423                 dec = psVectorAlloc(num, PS_TYPE_F32);
    424                 mag = psVectorAlloc(num, PS_TYPE_F32);
    425             }
    426 
    427             for (long i = 0; i < num; i++) {
    428                 ra->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * xSize ; // x position
    429                 dec->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * ySize; // y position
    430                 mag->data.F32[oldSize + i] = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
    431             }
    432         }
     25    // XXX put this in a wrapper (add to array of stars)
     26    if (type == PPSIM_TYPE_OBJECT) {
     27        ppSimMakeStars (stars, fpa, config, rng);
    43328    }
    43429
    43530    pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA
    43631
    437     // Update FPA concepts
    438     const char *typeStr;                // Exposure type String
    439     switch (type) {
    440       case PPSIM_TYPE_BIAS:   typeStr = "BIAS";   break;
    441       case PPSIM_TYPE_DARK:   typeStr = "DARK";   break;
    442       case PPSIM_TYPE_FLAT:   typeStr = "FLAT";   break;
    443       case PPSIM_TYPE_OBJECT: typeStr = "OBJECT"; break;
    444       default:
    445         psAbort("Should never get here.");
    446     }
    447     psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBSTYPE", PS_META_REPLACE,
    448                      "Observation type", typeStr);
    449     psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.OBJECT", PS_META_REPLACE,
    450                      "Observation name", typeStr);
    451     psMetadataAddF32(fpa->concepts, PS_LIST_TAIL, "FPA.EXPTIME", PS_META_REPLACE,
    452                      "Exposure time (sec)", expTime);
    453     psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTERID", PS_META_REPLACE, "Filter name", filter);
    454     psMetadataAddStr(fpa->concepts, PS_LIST_TAIL, "FPA.FILTER", PS_META_REPLACE, "Filter name", filter);
    455 
     32    ppSimUpdateConceptsFPA (fpa, config);
    45633
    45734    pmChip *chip;                       // Chip from FPA
    45835    while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
    459 
    460         int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
    461         int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
    462         int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY");
    463         int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY");
    464 
    465         // Correct chip offsets so that boresight is in the middle of the FPA
    466         x0Chip -= 0.5 * (bounds->x1 - bounds->x0);
    467         y0Chip -= 0.5 * (bounds->y1 - bounds->y0);
    46836
    46937        pmCell *cell;                   // Cell from chip
    … …  
    47341            int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
    47442            int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
    475             int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
    476             int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
    477             int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
    478             int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
    479 
    480             float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // CCD gain, e/ADU
    481             if (isnan(gain)) {
    482                 psWarning("CELL.GAIN is not set; reverting to recipe value GAIN.");
    483                 gain = psMetadataLookupF32(&mdok, recipe, "GAIN");
    484                 if (!mdok) {
    485                     psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe.");
    486                     psFree(bounds);
    487                     psFree(rng);
    488                     return PS_EXIT_CONFIG_ERROR;
    489                 }
    490             }
    491             float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
    492             if (isnan(readnoise)) {
    493                 psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
    494                 readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
    495                 if (!mdok) {
    496                     psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
    497                     psFree(bounds);
    498                     psFree(rng);
    499                     return PS_EXIT_CONFIG_ERROR;
    500                 }
    501             }
    502 
    503             float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
    504             if (isnan(saturation)) {
    505                 psWarning("CELL.SATURATION is not set; reverting to recipe value SATURATION.");
    506                 readnoise = psMetadataLookupF32(&mdok, recipe, "SATURATION");
    507                 if (!mdok) {
    508                     psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find SATURATION in recipe.");
    509                     psFree(bounds);
    510                     psFree(rng);
    511                     return PS_EXIT_CONFIG_ERROR;
    512                 }
    513             }
    51443
    51544            int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction
    51645            if (readdir != 1) {
    51746                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.READDIR = 1 is the only supported mode.");
    518                 psFree(bounds);
    51947                psFree(rng);
    52048                return PS_EXIT_CONFIG_ERROR;
    52149            }
    52250
    523             psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions
    524             int numBias = psListLength(biassec); // Number of bias sections
    525             psVector *biasCols;
    526             if (numBias > 0) {
    527                 biasCols = psVectorAlloc(numBias, PS_TYPE_S32);
    528                 psListIterator *iter = psListIteratorAlloc(biassec, PS_LIST_HEAD, false); // Iterator
    529                 psRegion *bias;         // Bias region, from iteration
    530                 int i = 0;              // Counter
    531                 while ((bias = psListGetAndIncrement(iter))) {
    532                     biasCols->data.S32[i++] = bias->x1 = bias->x0;
    533                 }
    534             } else {
    535                 biasCols = psVectorAlloc(1, PS_TYPE_S32);
    536                 biasCols->data.S32[0] = psMetadataLookupS32(&mdok, recipe, "OVERSCAN.SIZE");
    537                 if (!mdok) {
    538                     psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find OVERSCAN.SIZE in recipe.");
    539                     psFree(biasCols);
    540                     psFree(bounds);
    541                     psFree(rng);
    542                     return PS_EXIT_CONFIG_ERROR;
    543                 }
    544                 psRegion *newBias = psRegionAlloc(NAN, NAN, NAN, NAN); // New bias region, to be set later
    545                 biassec = psListAlloc(newBias);
    546                 psFree(newBias);
    547                 psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE,
    548                               "Bias sections", biassec);
    549                 psFree(biassec);
    550                 numBias = 1;
    551             }
     51            psVector *biasCols = ppSimMakeBiassec (cell, config);
    55252
    553             // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately
    554 
    555             // TO DO: Decide if cell is to be video readout
    556             int numReadouts = 1;        // Number of readouts in cell
    557 
     53            int numReadouts = 1;
    55854            for (int i = 0; i < numReadouts; i++) {
    55955                pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell
    56056
     57                // TO DO: Decide if cell is to be windowed, reduce numCols, numRows appropriately
    56158                psImage *signal = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Signal in pixels
    56259                psImage *variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
    563 
    564                 // Polynomial for bias
    565                 psPolynomial1D *biasPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, biasOrder);
    566                 for (int j = 0; j < biasOrder + 1; j++) {
    567                     biasPoly->coeff[j] = biasRange * psRandomGaussian(rng);
    568                 }
    569                 psVector *rowBias = psVectorAlloc(numRows, PS_TYPE_F32); // Bias value, per row
    570                 int biasOffset = 0.5 * numRows * psRandomUniform(rng); // Offset to prevent common pattern
    57160
    57261                psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
    … …  
    57564                }
    57665
    577                 for (int y = 0; y < numRows; y++) {
    578                     // Adjust bias level for this row
    579                     rowBias->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) /
    580                                                               (float)numRows - 0.5) + biasLevel;
    581                     float yFPA = cell2fpa(y, y0Cell, yParityCell, binning, y0Chip, yParityChip) * 2.0 /
    582                         (bounds->y1 - bounds->y0); // Relative y position in FPA
     66                psVector *biasRows = ppSimMakeBias (signal, variance, cell, config, rng);
     67                if (type == PPSIM_TYPE_BIAS) goto done;
     68               
     69                ppSimMakeDark (signal, variance, config);
     70                if (type == PPSIM_TYPE_DARK) goto done;
    58371
    584                     for (int x = 0; x < numCols; x++) {
    585                         float xFPA = cell2fpa(x, x0Cell, xParityCell, binning, x0Chip, xParityChip) * 2.0 /
    586                             (bounds->x1 - bounds->x0); // Relative x position in FPA
     72                ppSimMakeSky (signal, variance, expCorr, type, config, fpa, chip, cell);
     73                if (type == PPSIM_TYPE_FLAT) goto done;
    58774
    588                         // Bias level
    589                         signal->data.F32[y][x] = rowBias->data.F32[y];
    590                         variance->data.F32[y][x] = PS_SQR(readnoise);
     75                if (type == PPSIM_TYPE_OBJECT) {
     76                    ppSimInsertStars (signal, variance, expCorr, stars, config, chip, cell);
     77                }
    59178
    592                         if (type == PPSIM_TYPE_BIAS) {
    593                             continue;
    594                         }
     79            done:
     80                ppSimAddNoise(signal, variance, config, cell, rng);
     81                ppSimSaturate(signal, config, cell);
    59582
    596                         // Dark current
    597                         float darkCurrent = darkRate * expTime; // Dark current accumulated
    598                         signal->data.F32[y][x] += darkCurrent;
    599                         variance->data.F32[y][x] += darkCurrent;
    600 
    601                         if (type == PPSIM_TYPE_DARK) {
    602                             continue;
    603                         }
    604 
    605                         // Shutter: adjust exposure time
    606                         float realExpTime = expTime + shutterTime * (xFPA + yFPA + 2.0) / 4.0;
    607 
    608                         // Gaussian flat-field over the FPA
    609                         float flatValue = exp(-0.5 / PS_SQR(flatSigma) *
    610                                               (PS_SQR(yFPA) + PS_SQR(xFPA))) / flatSigma / sqrt(M_PI);
    611                         if (type == PPSIM_TYPE_FLAT) {
    612                             float flatFlux = flatRate * flatValue * realExpTime; // Flux from flat-field
    613                             signal->data.F32[y][x] += flatFlux;
    614                             variance->data.F32[y][x] += flatFlux;
    615                             continue;
    616                         }
    617                         expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
    618                         // Sky background
    619                         float skyFlux = skyRate * realExpTime * flatValue; // Flux from sky
    620                         signal->data.F32[y][x] += skyFlux;
    621                         variance->data.F32[y][x] += skyFlux;
    622 
    623                         // TO DO: Add fringes
    624 
    625                     }
    626                 }
    627                 psFree(biasPoly);
    628 
    629                 // Add stars
    630                 if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) {
    631                     // Rough noise estimate, appropriate for entire cell
    632                     float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
    633 
    634                     for (long j = 0; j < ra->n; j++) {
    635                         // Position on the cell and peak flux
    636                         float x = fpa2cell(ra->data.F32[j], x0Cell, xParityCell, binning,
    637                                            x0Chip, xParityChip);
    638                         float y = fpa2cell(dec->data.F32[j], y0Cell, yParityCell, binning,
    639                                            y0Chip, yParityChip);
    640                         float flux = mag->data.F32[j];
    641                         star(signal, variance, x, y, flux, roughNoise, seeing, expCorr);
    642                     }
    643                 }
    644                 readout->image = addNoise(signal, variance, rng, gain);
    645                 saturate(readout->image, saturation);
     83                readout->image = signal;
    64684                psFree(variance);
    647 
    64885                psFree(expCorr);
    64986
    650                 // Add overscan
    651                 for (int j = 0; j < numBias; j++) {
    652                     psImage *signal = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Overscan
    653                     for (int y = 0; y < signal->numRows; y++) {
    654                         for (int x = 0; x < signal->numCols; x++) {
    655                             signal->data.F32[y][x] = rowBias->data.F32[y];
    656                         }
    657                     }
    658                     psImage *variance = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Variance
    659                     psImageInit(variance, PS_SQR(readnoise));
    660                     addNoise(signal, variance, rng, gain);
    661                     psFree(variance);
    662 
    663                     psListAdd(readout->bias, PS_LIST_TAIL, signal);
    664                     psFree(signal);     // Drop reference
    665                 }
    666 
    667                 psFree(rowBias);
     87                ppSimAddOverscan (readout, config, biasCols, biasRows, rng);
     88                psFree(biasRows);
    66889
    66990                readout->data_exists = true;
    … …  
    67495            psFree(biasCols);
    67596
    676 
    677             psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,
    678                              "Exposure time (sec)", expTime);
    679             psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.DARKTIME", PS_META_REPLACE,
    680                              "Dark time (sec)", expTime);
    681             psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.XBIN", PS_META_REPLACE,
    682                              "Binning in x", binning);
    683             psMetadataAddS32(cell->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE,
    684                              "Binning in y", binning);
     97            ppSimUpdateConceptsCell (cell, config);
    68598
    68699            if (cell->hdu) {
    687                 cell->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning,
    688                                               fpa2cell(0.0, x0Cell, xParityCell, binning,
    689                                                        x0Chip, xParityChip),
    690                                               fpa2cell(0.0, y0Cell, yParityCell, binning,
    691                                                        y0Chip, yParityChip),
    692                                               xParityCell * xParityChip, yParityCell * yParityChip);
     100                ppSimInitHeader(config, NULL, NULL, cell);
    693101            }
    694             cell->data_exists = true;
    695102
    696103            if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    … …  
    698105                psFree(rng);
    699106                psFree(view);
    700                 psFree(bounds);
    701107                return PS_EXIT_SYS_ERROR;
    702108            }
    … …  
    704110
    705111        if (chip->hdu) {
    706             chip->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning,
    707                                           fpa2cell(0.0, 0, 1, binning, x0Chip, xParityChip),
    708                                           fpa2cell(0.0, 0, 1, binning, y0Chip, yParityChip),
    709                                           xParityChip, yParityChip);
     112            ppSimInitHeader(config, NULL, chip, NULL);
    710113        }
    711         chip->data_exists = true;
    712114
    713115        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    … …  
    715117            psFree(rng);
    716118            psFree(view);
    717             psFree(bounds);
    718119            return PS_EXIT_SYS_ERROR;
    719120        }
    … …  
    722123
    723124    if (fpa->hdu) {
    724         fpa->hdu->header = wcsHeader(ra0, dec0, pa, scale * binning, 0.5 * (bounds->x1 - bounds->x0),
    725                                      0.5 * (bounds->y1 - bounds->y0), 1, 1);
     125        ppSimInitHeader(config, fpa, NULL, NULL);
    726126    }
    727127
    … …  
    730130        psFree(rng);
    731131        psFree(view);
    732         psFree(bounds);
    733132        return PS_EXIT_SYS_ERROR;
    734133    }
    735134
    736     psFree(ra);
    737     psFree(dec);
    738     psFree(mag);
    739135    psFree(rng);
    740136    psFree(view);
    741     psFree(bounds);
    742137
    743138    return 0;
Note: See TracChangeset for help on using the changeset viewer.