Changeset 41173

Timestamp:

Nov 27, 2019, 12:05:20 PM (6 years ago)

Author:

eugene

Message:

merge fixes from EAM dev branch ipp-20191011

Location:

trunk/ppSim/src

Files:

• Makefile.am (modified) (1 diff)
• ppSim.h (modified) (1 diff)
• ppSimCreate.c (modified) (1 diff)
• ppSimLoadStars.c (modified) (2 diffs)
• ppSimLoop.c (modified) (1 diff)
• ppSimMakeGalaxies.c (modified) (3 diffs)
• ppSimMakeStarCluster.c (copied) (copied from branches/eam_branches/ipp-20191011/ppSim/src/ppSimMakeStarCluster.c )
• ppSimMakeStarGrid.c (modified) (1 diff)
• ppSimMakeStars.c (modified) (4 diffs)
• ppSimUtils.c (modified) (2 diffs)

trunk/ppSim/src/Makefile.am

@@ -15 +15 @@
         ppSimLoadStars.c          \
         ppSimMakeStars.c          \
+        ppSimMakeStarCluster.c    \
         ppSimMakeStarGrid.c       \
         ppSimMakeGalaxies.c       \

trunk/ppSim/src/ppSim.h

@@ -103 +103 @@
 bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config);
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng);
+bool ppSimMakeStarCluster(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng);
 bool ppSimMakeStarGrid(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng);
 bool ppSimInsertStars (pmReadout *readout, psImage *expCorr, psArray *stars, pmConfig *config);

trunk/ppSim/src/ppSimCreate.c

@@ -159 +159 @@
         return false;
     }
-    simSources->save = true;
+
+    // XXXX TEST this is causing trouble for unknown reasons -- output step is blowing up.
+    // simSources->save = true;
+    simSources->save = false;
 
     // if we have loaded an input image, we derive certain values from the image, if possible

trunk/ppSim/src/ppSimLoadStars.c

@@ -24 +24 @@
     float pa      = psMetadataLookupF32(NULL, recipe, "PA");        // Position angle (radians)
     float seeing  = psMetadataLookupF32(NULL, recipe, "SEEING");    // Seeing SIGMA (pixels)
-    float scale   = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+    float scale   = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel)
+    float scaleRad = scale * M_PI / 3600.0 / 180.0; // Plate scale in radians/pixel for WCS below
     float expTime = psMetadataLookupF32(NULL, recipe, "EXPTIME");   // Exposure time (sec)
 
     // Size of FPA
     psRegion *bounds = ppSimFPABounds (fpa);
-    float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;
+    float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scaleRad;
 
     float x0fpa = 0.5*(bounds->x0 + bounds->x1);
@@ -53 +54 @@
     stars = psArrayRealloc (stars, refStars->n);
 
-    psProjection *proj = psProjectionAlloc(ra0, dec0, scale, scale, PS_PROJ_TAN); // Projection
+    psProjection *proj = psProjectionAlloc(ra0, dec0, scaleRad, scaleRad, PS_PROJ_TAN); // Projection
 
     // Conversion loop

trunk/ppSim/src/ppSimLoop.c

@@ -51 +51 @@
 
         // Add random stars
-        if (!ppSimMakeStarGrid (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make random stars");
+        if (!ppSimMakeStarGrid (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make star grid");
+
+        // Add random stars
+        if (!ppSimMakeStarCluster (stars, fpa, config, rng)) ESCAPE (PS_ERR_UNKNOWN, "failed to make cluster stars");
 
         // Add random galaxies

trunk/ppSim/src/ppSimMakeGalaxies.c

@@ -45 +45 @@
     float zp              = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point
     float seeing          = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix)
-    float scale           = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+    float scale           = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel)
     float skyRate         = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate
     if (isnan(skyRate)) {
         float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
         skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
+        // skyMags is in mags / square arcsec so scale must be arcsec / pixel
     }
 
@@ -94 +95 @@
 
     // Normalisation, set by the specified stellar density at the specified bright magnitude
-    float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI);
+    float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale / 3600.0);
     float normLum =  refSum / (galaxyAlpha * logf(10.0) * powf(10.0, (galaxyAlpha * brightMag)));
     float normScale = normLum * galaxyAlpha * logf(10.0);
@@ -143 +144 @@
                 rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal;
                 float scale = (galaxyRmajorMin + rndValue * galaxyRmajorSlope);
-
                 galaxy->Rmaj  = scale;
 

trunk/ppSim/src/ppSimMakeStarGrid.c

@@ -22 +22 @@
     float zp           = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point
     float seeing       = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
-    float scale        = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+    float scale        = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel)
 
     if (isnan(expTime))  ESCAPE("EXPTIME is not defined");

trunk/ppSim/src/ppSimMakeStars.c

@@ -27 +27 @@
     float zp           = psMetadataLookupF32(&status, recipe, "ZEROPOINT"); // Photometric zero point
     float seeing       = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
-    float scale        = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+    float scale        = psMetadataLookupF32(&status, recipe, "PIXEL.SCALE"); // Plate scale (arcsec/pixel)
+    // there has been some confusion over pixel scale: PIXEL.SCALE is supplied in arcsec / pixel.
+    // In some other places (ppSimLoadStars.c, ppSimUtils.c) it is needed in radius for WCS conversion.
+    // in the past, it was converted to radians on load (but applied inconsisently elsewhere)
+    // Now the radian version is carried as scaleRad to be more explicit
 
     if (isnan(darkRate)) darkRate = 0.0;
@@ -42 +46 @@
         float skyMags = psMetadataLookupF32(&status, recipe, "SKY.MAGS");  assert (status);
         skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
+        // skyMags is in mags / square arcsec, so scale must be in arcsec / pixel
     }
 
@@ -61 +66 @@
         if (!status) {
             refMag = brightMag;
-            refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI);
+            refSum = starsDensity * xSize * ySize * PS_SQR(scale / 3600.0);
         }
     } else {
         refMag = brightMag;
-        refSum = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI);
+        refSum = starsDensity * xSize * ySize * PS_SQR(scale / 3600.0);
     }
     psTrace("ppSim", 6, "refMag: %f, refSum: %f\n", refMag, refSum);
@@ -120 +125 @@
     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars between %f and %f mag\n", nTotal, brightMag, faintMag);
 
+    if (nTotal > 1e7) {
+      psLogMsg("ppSim", PS_LOG_INFO, "Trying to add %d stars (far more than 10M stars!), giving up\n", (int) nTotal);
+      exit (2);
+    }
+
     long oldSize = stars->n;
     psArrayRealloc (stars, stars->n + nTotal);

trunk/ppSim/src/ppSimUtils.c

@@ -16 +16 @@
     float pa    = psMetadataLookupF32(NULL, recipe, "PA");  // Position angle (radians)
     float scale = psMetadataLookupF32(NULL, recipe, "PIXEL.SCALE"); // plate scale in arcsec / pixel
-    scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel
+    float scaleRad = scale * M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel
 
     int binning = psMetadataLookupS32(NULL, recipe, "BINNING"); // Binning in x and y
@@ -61 +61 @@
     psMetadata *header = psMetadataAlloc(); // Header, to return
     pmAstromWCS *wcs = pmAstromWCSAlloc(1, 1); // WCS structure
-    wcs->toSky = psProjectionAlloc(ra0, dec0, scale * xParity, scale * yParity, PS_PROJ_TAN);
-    wcs->cdelt1 = scale * PM_DEG_RAD * xParity;
-    wcs->cdelt2 = scale * PM_DEG_RAD * yParity;
+    wcs->toSky = psProjectionAlloc(ra0, dec0, scaleRad * xParity, scaleRad * yParity, PS_PROJ_TAN);
+    wcs->cdelt1 = scaleRad * PM_DEG_RAD * xParity;
+    wcs->cdelt2 = scaleRad * PM_DEG_RAD * yParity;
     wcs->crpix1 = x0;
     wcs->crpix2 = y0;