Changeset 32350 for trunk/ppSim/src/ppSimMakeGalaxies.c

Timestamp:

Sep 6, 2011, 1:34:52 PM (15 years ago)

Author:

eugene

Message:

fix the mags of input stars and galaxies (prove we inserted the amount of flux we expected); realistic galaxy luminosity function

File:

• trunk/ppSim/src/ppSimMakeGalaxies.c (modified) (11 diffs)

trunk/ppSim/src/ppSimMakeGalaxies.c

@@ -11 +11 @@
     if (!galaxyFake) return true;
 
-    float galaxyLum       = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
-    float galaxyMag       = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
+    float galaxyAlpha     = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
+    float brightMag       = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
     float galaxyDensity   = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes
 
@@ -20 +20 @@
     bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes
     float galaxyGridPeak  = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.PEAK"); // peak flux of fakes
+    float galaxyGridMag   = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.MAG"); // peak flux of fakes
     
+    float galaxyGridXoff  = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.XOFF"); // centroid offset
+    float galaxyGridYoff  = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.YOFF"); // centroid offset
+
     float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes
     float galaxyRmajorMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MIN"); // Density of fakes
@@ -35 +39 @@
     float galaxyIndexMax  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes
 
-    float darkRate        = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate
+    // float darkRate     = psMetadataLookupF32(&mdok, recipe, "DARK.RATE"); // Dark rate
     float expTime         = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time
     float zp              = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point
-    float seeing          = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix)
+    // 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 skyRate         = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate
@@ -44 +48 @@
         float skyMags = psMetadataLookupF32(&mdok, recipe, "SKY.MAGS");  assert (mdok);
         skyRate = scale * scale * ppSimMagToFlux (skyMags, zp);
+    }
+
+    // determine the galaxy model
+    char *modelName = psMetadataLookupStr(&mdok, recipe, "GALAXY.MODEL"); // galaxy model name
+    pmModelType type = pmModelClassGetType (modelName);
+    if (type == -1) {
+        psError (PS_ERR_UNKNOWN, false, "invalid model name");
+        return false;
     }
 
@@ -64 +76 @@
     // 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);
-        return false;
-    }
-
-    // Peak fluxes: faintest and brightest levels for random galaxy
-    float faint = 0.1 * 10.0 * sqrt(2.0*M_PI) * seeing * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
-    float bright = powf(10.0, -0.4 * (galaxyMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
-    if (bright < faint) {
-        psLogMsg("ppSim", PS_LOG_INFO,
-                 "Image noise is above brightest random galaxy --- no random galaxy added.");
+    // 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);
+    //  return false;
+    // }
+
+    // faintest and brightest magnitudes for random galaxies
+    // float brightFlux = ppSimMagToFlux (brightMag, zp) * expTime;
+    float faintMag = brightMag + 5.0;
+    if (brightMag > faintMag) {
+        psLogMsg("ppSim", PS_LOG_INFO, "Image noise is above brightest random galaxy --- no random galaxy added.");
         return true;
     }
 
     // Normalisation, set by the specified stellar density at the specified bright magnitude
-    float norm = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
-        powf(bright, galaxyLum);
+    float refSum = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI);
+    float normLum =  refSum / (galaxyAlpha * logf(10.0) * powf(10.0, (galaxyAlpha * brightMag)));
+    float normScale = normLum * galaxyAlpha * logf(10.0);
 
     // Total number of galaxy down to the faint flux end
-    long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5;
+    long nTotal = 0;
+    if (galaxyAlpha < 0.1) {
+        nTotal = 100;
+    } else {
+        nTotal = normScale * powf (10.0, (galaxyAlpha * faintMag));
+    }
 
     if (galaxyGrid) {
-        num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
+        fprintf (stderr, "galaxy grid inst mag: %f\n", galaxyGridMag - zp - 2.5*log10(expTime));
+        nTotal = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
     
-        psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", num);
+        psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", nTotal);
 
         float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin);
@@ -98 +116 @@
 
         int i = 0;
-        float refNorm = 1.0;
-        float ourNorm = 1.0;
+        // float refNorm = 1.0;
+        // float ourNorm = 1.0;
 
         for (long iy = 0.5*galaxyGridDY; iy < ySize; iy += galaxyGridDY) {
@@ -106 +124 @@
 
                 // make fpa center of distribution
-                galaxy->x    = ix;
-                galaxy->y    = iy;
-                galaxy->peak = galaxyGridPeak;
+                // center the galaxy peak on the center of a pixel
+                galaxy->x    = ix + 0.5 + galaxyGridXoff;
+                galaxy->y    = iy + 0.5 + galaxyGridYoff;
+
+                float galaxyGridFlux = expTime * powf(10.0, -0.4 * (galaxyGridMag - zp));
 
                 // galaxyIndex from user should be for function of this form: exp(-r^(1/n))
@@ -115 +135 @@
                 float rndValue;
 
-                rndValue = galaxyGridRandom ? drand48() : i / (float) num;
+                rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal;
                 float index = (galaxyIndexMin  + rndValue * galaxyIndexSlope);
                 galaxy->index = 0.5/index; // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n)
 
-                rndValue = galaxyGridRandom ? drand48() : i / (float) num;
+                rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal;
                 float scale = (galaxyRmajorMin + rndValue * galaxyRmajorSlope);
 
-                // for a sersic model, 
-                float bn = 1.9992*index - 0.3271;
-                float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
-                float Io = exp(bn);
-
                 galaxy->Rmaj  = scale;
 
-                rndValue = galaxyGridRandom ? drand48() : i / (float) num;
+                rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal;
                 galaxy->Rmin  = (galaxyARatioMin + rndValue * galaxyARatioSlope) * galaxy->Rmaj;
 
-                rndValue = galaxyGridRandom ? drand48() : i / (float) num;
+                rndValue = galaxyGridRandom ? drand48() : i / (float) nTotal;
                 galaxy->theta = (galaxyThetaMin  + rndValue * galaxyThetaSlope);
                 
-                if (i == 0) {
-                    refNorm = Io*scale*scale;
-                } 
-                ourNorm = refNorm / (Io*scale*scale);
-                galaxy->peak *= ourNorm;
-
-                if (0) {
-                    fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f, theta: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io, galaxy->theta);
+                bool sersicType = false;
+                sersicType |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
+                sersicType |= (type == pmModelClassGetType ("PS_MODEL_EXP"));
+                sersicType |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
+                if (sersicType) {
+                    if (type == pmModelClassGetType ("PS_MODEL_DEV")) index = 4.0;
+                    if (type == pmModelClassGetType ("PS_MODEL_EXP")) index = 1.0;
+                        
+                    // for a sersic model, 
+                    float bn = 1.9992*index - 0.3271;
+                    float Io = exp(bn);
+                    
+                    // the integral of a Sersic has an analytical form as follows:
+                    float logGamma = lgamma(2.0*index);
+                    float bnFactor = pow(bn, 2.0*index);
+                    float norm = 2.0 * M_PI * PS_SQR(galaxy->Rmaj) * index * Io * exp(logGamma) / bnFactor;
+                    
+                    // XXX probably should limit the allowed thinness of a galaxy
+                    galaxy->peak = galaxyGridFlux / norm / (galaxy->Rmin / galaxy->Rmaj);
+                } else {
+                    galaxy->peak = galaxyGridPeak;
                 }
+
+                // if (0) {
+                //     fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f, theta: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io, galaxy->theta);
+                // }
 
                 psArrayAdd (galaxies, 100, galaxy);
@@ -151 +183 @@
     } else {    
 
+        fprintf (stderr, "galaxy ref inst mag: %f\n", brightMag - zp - 2.5*log10(expTime));
         float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin);
         float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin);
@@ -156 +189 @@
         float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin);
 
-        psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", num,
-                 -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
-
-        for (long i = 0; i < num; i++) {
+        psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", nTotal, faintMag);
+
+        for (long i = 0; i < nTotal; i++) {
             ppSimGalaxy *galaxy = ppSimGalaxyAlloc ();
 
@@ -170 +202 @@
 
             float scale = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin);
-
-            // for a sersic model, 
-            float bn = 1.9992*index - 0.3271;
-            float fR = 1.0 / (sqrt(2.0) * pow (bn, index));
-            float Io = exp(bn);
-
             galaxy->Rmaj  = scale;
-
             galaxy->Rmin  = (PS_SQR(psRandomUniform(rng)) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj;
-
             galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope  + galaxyThetaMin);
 
-            // the flux is only approximately correct (scaling of the peak is wrong)
-            // elsewhere (ppSimInsertGalaxies.c) we calculate the true galaxy flux
-            galaxy->flux = expf((logf(i + 1) - logf(norm)) / galaxyLum); // Peak flux
-            galaxy->peak = galaxy->flux / (2.0*M_PI*PS_SQR(seeing));
-
-            if (0) {
-              fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io);
+            bool sersicType = false;
+            sersicType |= (type == pmModelClassGetType ("PS_MODEL_SERSIC"));
+            sersicType |= (type == pmModelClassGetType ("PS_MODEL_EXP"));
+            sersicType |= (type == pmModelClassGetType ("PS_MODEL_DEV"));
+            if (sersicType) {
+                if (type == pmModelClassGetType ("PS_MODEL_DEV")) index = 4.0;
+                if (type == pmModelClassGetType ("PS_MODEL_EXP")) index = 1.0;
+                
+                // for a sersic model, 
+                float bn = 1.9992*index - 0.3271;
+                float Io = exp(bn);
+                
+                // the integral of a Sersic has an analytical form as follows:
+                float logGamma = lgamma(2.0*index);
+                float bnFactor = pow(bn, 2.0*index);
+                float norm = 2.0 * M_PI * PS_SQR(galaxy->Rmaj) * index * Io * exp(logGamma) / bnFactor;
+                    
+                // XXX probably should limit the allowed thinness of a galaxy
+                float galaxyMag =  PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / galaxyAlpha)));
+
+                galaxy->flux = ppSimMagToFlux (galaxyMag, zp) * expTime;
+                galaxy->peak = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj);
+            }
+
+            if (1) {
+                fprintf (stderr, "Rmaj: %f, Rmin: %f, theta: %f, scale: %f, index: %f, peak: %f, mag: %f\n", galaxy->Rmaj, galaxy->Rmin, galaxy->theta*180.0/M_PI, scale, index, galaxy->peak, -2.5*log10(galaxy->flux));
             }