Index: trunk/ppSim/src/ppSimMakeStars.c =================================================================== --- trunk/ppSim/src/ppSimMakeStars.c (revision 16498) +++ trunk/ppSim/src/ppSimMakeStars.c (revision 17557) @@ -1,35 +1,33 @@ # include "ppSim.h" - -#define FAINT_FUDGE_FACTOR 0.4 // Fraction of the noise in a (boxcar) PSF for the faint limit - bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) { - bool mdok; + bool status; assert (stars); - psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe + psMetadata *recipe = psMetadataLookupMetadata(&status, config->recipes, PPSIM_RECIPE); // Recipe - bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes + bool starsFake = psMetadataLookupBool(&status, recipe, "STARS.FAKE"); // Density of fakes if (!starsFake) return true; - float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope - float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude - float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes + bool starsReal = psMetadataLookupBool(&status, recipe, "STARS.REAL"); // were real stars generated? + bool matchDensity = psMetadataLookupBool(&status, recipe, "MATCH.DENSITY"); // match real star density? - float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate - float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time + float starsAlpha = psMetadataLookupF32(&status, recipe, "STARS.LUM"); // Star luminosity func slope + float starsDensity = psMetadataLookupF32(&status, recipe, "STARS.DENSITY"); // Density of fakes + float brightMag = psMetadataLookupF32(&status, recipe, "STARS.MAG"); // Star brightest magnitude + float nSigmaLim = psMetadataLookupF32(&status, recipe, "STARS.SIGMA.LIM"); // significance of faintest stars - float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point - float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix) - float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) + float darkRate = psMetadataLookupF32(&status, recipe, "DARK.RATE"); // Dark rate + float expTime = psMetadataLookupF32(&status, recipe, "EXPTIME"); // Exposure time + float zp = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point + float seeing = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels) + float scale = psMetadataLookupF32(&status, recipe, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel) - float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate + float skyRate = psMetadataLookupF32(&status, recipe, "SKY.RATE"); // Sky rate if (isnan(skyRate)) { - float skyMags = psMetadataLookupF32(&mdok, config->arguments, "SKY.MAGS"); assert (mdok); - skyRate = scale * scale * pow (10.0, -0.4*(skyMags - zp)); + float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS"); assert (status); + skyRate = scale * scale * ppSimMagToFlux (skyMags, zp); } - - if (starsDensity <= 0) return true; // Size of FPA @@ -40,8 +38,23 @@ psFree(bounds); + // choose reference magnitude & density to set normalization + float refMag = 0; + float refSum = 0; + if (starsReal && matchDensity) { + refMag = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.MAG.PEAK"); // Star brightest magnitude + refSum = psMetadataLookupF32(&status, fpa->concepts, "STARS.REAL.SUM.PEAK"); // Star brightest magnitude + assert (status); + } else { + refMag = brightMag; + refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI); + } + psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum); + + if (refSum <= 0) return true; + // 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) { + float readnoise = psMetadataLookupF32(&status, recipe, "READNOISE"); // Default read noise + if (!status) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe."); psFree(bounds); @@ -50,10 +63,18 @@ // Faintest and brightest (integrated) fluxes (actually fluence, since integrated) for random stars - float faint = FAINT_FUDGE_FACTOR * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime) * - (2.0 * sqrt(2.0 * log(2.0)) * seeing); // Faint limit is related to the noise in a (boxcar) PSF - float bright = powf(10.0, -0.4 * (starsMag - zp)) * expTime; // Bright limit is specified by user as mag - psTrace("ppSim", 6, "Faint limit: %f\n", faint); - psTrace("ppSim", 6, "Bright limit: %f\n", bright); - if (bright < faint) { + + // faint limit is set by detection limit, in turn set by sky noise + float skySigma = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime); + float skyNoise = ppSimStarSkyNoise (skySigma, seeing); + + // total flux of faintest star: + float faintCounts = nSigmaLim * skyNoise; + float faintMag = ppSimFluxToMag ((faintCounts / expTime), zp); + + float brightCounts = ppSimMagToFlux (brightMag, zp) * expTime; // Bright limit is specified by user as mag + + psTrace("ppSim", 6, "Faint limit: %f counts = %f mags\n", faintCounts, faintMag); + psTrace("ppSim", 6, "Bright limit: %f counts = %f mags\n", brightCounts, brightMag); + if (brightCounts < faintCounts) { psLogMsg("ppSim", PS_LOG_INFO, "Image noise is above brightest random star --- no random stars added."); @@ -61,18 +82,23 @@ } - // 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); + // given log_10 (dN / dmag) = alpha mag + beta + // or dN / dmag = No 10^(alpha mag), then: + // N(m < Mo) = alpha * No * log(10.0) * 10^(alpha*Mo) + + // XXX this needs to handle the case of a flat distribution for efficiency testing + + // Normalization, set by the specified stellar density at the specified bright magnitude + float norm = refSum / (starsAlpha * logf(10.0) * powf(10.0, (starsAlpha * refMag))); + float normScale = norm * starsAlpha * logf(10.0); // Total number of stars down to the faint flux end - long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5; + long nTotal = normScale * powf (10.0, (starsAlpha * faintMag)); - psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", - num, -2.5 * log10(faint / expTime) + zp); + psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", nTotal, faintMag); long oldSize = stars->n; - psArrayRealloc (stars, stars->n + num); + psArrayRealloc (stars, stars->n + nTotal); - for (long i = 0; i < num; i++) { + for (long i = 0; i < nTotal; i++) { ppSimStar *star = ppSimStarAlloc (); @@ -81,10 +107,13 @@ star->y = psRandomUniform(rng) * ySize; // y position - star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux - star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing)); + // XXX this needs to respect the bright limit + float starMag = log10((i + 1.0) / normScale) / starsAlpha; + + star->flux = ppSimMagToFlux (starMag, zp) * expTime; + star->peak = ppSimStarFluxToPeak (star->flux, seeing); stars->data[oldSize + i] = star; } - stars->n = oldSize + num; + stars->n = oldSize + nTotal; return true;