Changeset 17557 for trunk/ppSim/src/ppSimMakeStars.c

Timestamp:

May 7, 2008, 10:55:21 AM (18 years ago)

Author:

eugene

Message:

fix errors with star density normalizations and flux consistency

File:

• trunk/ppSim/src/ppSimMakeStars.c (modified) (5 diffs)

trunk/ppSim/src/ppSimMakeStars.c

@@ -1 +1 @@
 # 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 +38 @@
     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 +63 @@
 
     // 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 +82 @@
     }
 
-    // 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 +107 @@
         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;