Changeset 12907

Timestamp:

Apr 18, 2007, 12:49:57 PM (19 years ago)

Author:

Paul Price

Message:

Putting random stars over the entire focal plane, so they're not cut
off by cell borders, etc.

File:

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

trunk/ppSim/src/ppSimLoop.c

@@ -208 +208 @@
                  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 +247 @@
         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 +266 @@
 
     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 +393 @@
         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 +553 @@
                 }
 
-                // 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);