Index: /trunk/ppSim/src/ppSimLoop.c =================================================================== --- /trunk/ppSim/src/ppSimLoop.c (revision 12906) +++ /trunk/ppSim/src/ppSimLoop.c (revision 12907) @@ -208,5 +208,5 @@ float peak, // Peak flux of star float noise, // Rough noise estimate - float seeing, // Seeing for star + float seeing, // Seeing for star const psImage *correction // Exposure correction as a function of position ) @@ -247,39 +247,4 @@ psError(PS_ERR_UNKNOWN, false, "Unable to determine bounds of FPA"); return PS_EXIT_CONFIG_ERROR; - } - - psVector *ra = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.RA"); // Star RAs - psVector *dec = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.DEC"); // Star Declinations - psVector *mag = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.MAG"); // Star magnitudes - float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // Plate scale (arcsec/pixel) - scale *= M_PI / 3600.0 / 180.0; // Convert to radians/pixel - float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point - - if (ra && dec && mag) { - // Convert star ra,dec to pixels from FPA centre - float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians) - float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians) - float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians) - - // Conversion loop - for (long i = 0; i < ra->n; i++) { - float div = sin(dec->data.F32[i]) * sin(dec0) + - cos(dec->data.F32[i]) * cos(dec0) * cos(ra->data.F32[i] - ra0); // Divisor - - // Convert to x,y position on tangent plane, in pixels - float xi = cos(dec->data.F32[i]) * sin(ra->data.F32[i] - ra0) / div / scale; - float eta = (sin(dec->data.F32[i]) * cos(dec0) - - cos(dec->data.F32[i]) * sin(dec0) * cos(ra->data.F32[i] - ra0)) / div / scale; - - // Apply rotation - ra->data.F32[i] = cos(pa) * xi - sin(pa) * eta; - dec->data.F32[i] = sin(pa) * xi + cos(pa) * eta; - ra->data.F32[i] = xi; - dec->data.F32[i] = eta; - - // Convert magnitude to peak flux - mag->data.F32[i] = powf(10.0, -0.4 * (mag->data.F32[i] - zp)) / sqrt(2.0*M_PI) / seeing; - } } @@ -301,9 +266,103 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + psVector *ra = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.RA"); // Star RAs + psVector *dec = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.DEC"); // Star Declinations + psVector *mag = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.MAG"); // Star magnitudes + float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) + float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // Plate scale (arcsec/pixel) + scale *= M_PI / 3600.0 / 180.0; // Convert to radians/pixel + float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point + + // Add catalogue stars + if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) { + // Convert star ra,dec to pixels from FPA centre + float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians) + float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians) + float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians) + + // Conversion loop + for (long i = 0; i < ra->n; i++) { + float div = sin(dec->data.F32[i]) * sin(dec0) + + cos(dec->data.F32[i]) * cos(dec0) * cos(ra->data.F32[i] - ra0); // Divisor + + // Convert to x,y position on tangent plane, in pixels + float xi = cos(dec->data.F32[i]) * sin(ra->data.F32[i] - ra0) / div / scale; + float eta = (sin(dec->data.F32[i]) * cos(dec0) - + cos(dec->data.F32[i]) * sin(dec0) * cos(ra->data.F32[i] - ra0)) / div / scale; + + // Apply rotation + ra->data.F32[i] = cos(pa) * xi - sin(pa) * eta; + dec->data.F32[i] = sin(pa) * xi + cos(pa) * eta; + + // Convert magnitude to peak flux + mag->data.F32[i] = powf(10.0, -0.4 * (mag->data.F32[i] - zp)) / sqrt(2.0*M_PI) / seeing * expTime; + } + + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.RA", PS_META_REPLACE, "Right ascension", ra0); + psMetadataAddF64(fpa->concepts, PS_LIST_TAIL, "FPA.DEC", PS_META_REPLACE, "Declination", dec0); + } + + // Add random stars + if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) { + + // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it + // shouldn't be too bad. + float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise + if (!mdok) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); + psFree(bounds); + psFree(rng); + return PS_EXIT_CONFIG_ERROR; + } + + // Peak fluxes: faintest and brightest levels for random stars + float faint = 0.1 * 10.0 * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime; + if (bright < faint) { + psLogMsg("ppSim", PS_LOG_INFO, + "Image noise is above brightest random star --- no random stars added."); + } else { + // Size of focal plane + int xSize = bounds->x1 - bounds->x0; + int ySize = bounds->y1 - bounds->y0; + + // Normalisation, set by the specified stellar density at the specified bright magnitude + float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) / + powf(bright, starsLum); + + // Total number of stars down to the faint flux end + long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5; + + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num, + -2.5 * LOG10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp); + + // Extend the vectors + long oldSize; // Old size of ra, dec, mag vectors + if (ra && dec && mag) { + oldSize = ra->n; + ra = psVectorRealloc(ra, oldSize + num); + dec = psVectorRealloc(dec, oldSize + num); + mag = psVectorRealloc(mag, oldSize + num); + ra->n = dec->n = mag->n = oldSize + num; + } else { + oldSize = 0; + ra = psVectorAlloc(num, PS_TYPE_F32); + dec = psVectorAlloc(num, PS_TYPE_F32); + mag = psVectorAlloc(num, PS_TYPE_F32); + } + + for (long i = 0; i < num; i++) { + ra->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * xSize ; // x position + dec->data.F32[oldSize + i] = (psRandomUniform(rng) - 0.5) * ySize; // y position + mag->data.F32[oldSize + i] = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux + } + } + } + pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA - // Update concepts + // Update FPA concepts const char *typeStr; // Exposure type String switch (type) { @@ -334,39 +393,4 @@ x0Chip -= 0.5 * (bounds->x1 - bounds->x0); y0Chip -= 0.5 * (bounds->y1 - bounds->y0); - -#if 0 - // Add random stars - if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) { - // Peak fluxes: faintest and brightest levels for random stars - float faint = roughNoise; - float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime; - if (bright < faint) { - psLogMsg("ppSim", PS_LOG_INFO, - "Image noise is above brightest random star --- no random stars added."); - } else { - // Normalisation, set by the specified stellar density at the specified bright - // magnitude - float norm = starsDensity * numCols * numRows * PS_SQR(scale * 180.0 / M_PI) / - powf(bright, starsLum); - - // Total number of stars down to the faint flux end - long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5; - - printf("Adding %ld stars to %f\n", num, faint); - - for (long j = 0; j < num; j++) { - float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star - - // TO DO: put stars on chips, rather than cells, so they can overlap cell - // boundaries - - // Position on cell - float x = psRandomUniform(rng) * numCols; - float y = psRandomUniform(rng) * numRows; - star(signal, variance, x, y, flux, roughNoise, seeing, expCorr); - } - } - } -#endif pmCell *cell; // Cell from chip @@ -529,55 +553,17 @@ } - // Rough noise estimate, appropriate for entire cell - float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); - - - // XXX Not sure I've got the scaling right on the peak flux: it's multiplied by the expTime - // here, and the exposure up above. - - // Add specified (catalogue) stars + // Add stars if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) { + // Rough noise estimate, appropriate for entire cell + float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + for (long j = 0; j < ra->n; j++) { - // Position on the cell + // Position on the cell and peak flux float x = fpa2cell(ra->data.F32[j], x0Cell, xParityCell, x0Chip, xParityChip); float y = fpa2cell(dec->data.F32[j], y0Cell, yParityCell, y0Chip, yParityChip); - star(signal, variance, x, y, mag->data.F32[j] * expTime, - roughNoise, seeing, expCorr); + float flux = mag->data.F32[j]; + star(signal, variance, x, y, flux, roughNoise, seeing, expCorr); } } - - // Add random stars - if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) { - // Peak fluxes: faintest and brightest levels for random stars - float faint = roughNoise; - float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime; - if (bright < faint) { - psLogMsg("ppSim", PS_LOG_INFO, - "Image noise is above brightest random star --- no random stars added."); - } else { - // Normalisation, set by the specified stellar density at the specified bright - // magnitude - float norm = starsDensity * numCols * numRows * PS_SQR(scale * 180.0 / M_PI) / - powf(bright, starsLum); - - // Total number of stars down to the faint flux end - long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5; - - printf("Adding %ld stars to %f\n", num, faint); - - for (long j = 0; j < num; j++) { - float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star - - // TO DO: put stars on chips, rather than cells, so they can overlap cell - // boundaries - - // Position on cell - float x = psRandomUniform(rng) * numCols; - float y = psRandomUniform(rng) * numRows; - star(signal, variance, x, y, flux, roughNoise, seeing, expCorr); - } - } - } - readout->image = addNoise(signal, variance, rng, gain); saturate(readout->image, saturation);