Changeset 29060 for branches/sc_branches/trunkTest/ppSim/src/ppSimMakeGalaxies.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• ppSim/src/ppSimMakeGalaxies.c (modified) (5 diffs)

branches/sc_branches/trunkTest/ppSim/src/ppSimMakeGalaxies.c

@@ -18 +18 @@
     int galaxyGridDX      = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes
     int galaxyGridDY      = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
+    bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes
     
     float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes
@@ -43 +44 @@
 
     if (galaxyDensity <= 0) return true;
+
+    if (galaxyGridRandom) {
+        long A, B;
+        A = time(NULL);
+        for (B = 0; A == time(NULL); B++);
+        srand48(B);
+    }
 
     // Size of FPA
@@ -81 +89 @@
         psLogMsg("ppSim", PS_LOG_INFO, "Adding grid of %ld galaxies\n", num);
 
-        float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin) / num;
-        float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin) / num;
-        float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin) / num;
-        float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin) / num;
+        float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin);
+        float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin);
+        float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin);
+        float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin);
 
         int i = 0;
@@ -96 +104 @@
                 galaxy->y    = iy;
 
-                galaxy->peak = 20000;
+                galaxy->peak = 1000;
 
-                galaxy->index = (galaxyIndexMin  + i * galaxyIndexSlope);
-                galaxy->Rmaj  = (galaxyRmajorMin + i * galaxyRmajorSlope);
-                galaxy->Rmin  = (galaxyARatioMin + i * galaxyARatioSlope) * galaxy->Rmaj;
-                galaxy->theta = (galaxyThetaMin  + i * galaxyThetaSlope);
+                // galaxyIndex from user should be for function of this form: exp(-r^(1/n))
+                // galaxy->index = (1/2n)
+
+                float rndValue;
+
+                rndValue = galaxyGridRandom ? drand48() : i / num;
+                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 / num;
+                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 / num;
+                galaxy->Rmin  = (galaxyARatioMin + rndValue * galaxyARatioSlope) * galaxy->Rmaj;
+
+                rndValue = galaxyGridRandom ? drand48() : i / num;
+                galaxy->theta = (galaxyThetaMin  + rndValue * galaxyThetaSlope);
+
+                // galaxy->peak *= Io;
+
+                if (0) {
+                    fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io);
+                }
 
                 psArrayAdd (galaxies, 100, galaxy);
@@ -124 +158 @@
             galaxy->y    = psRandomUniform(rng) * ySize; // y position
 
-            galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum));
+            float index = (sqrt(psRandomUniform(rng)) * galaxyIndexSlope  + galaxyIndexMin);
+            galaxy->index = 0.5/index; // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n)
 
-            galaxy->index = (psRandomUniform(rng) * galaxyIndexSlope  + galaxyIndexMin);
+            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);
-            galaxy->Rmaj  = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin);
-            galaxy->Rmin  = (psRandomUniform(rng) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj;
 
+            // 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));
-            // galaxy->peak = 5000;
+
+            if (0) {
+              fprintf (stderr, "Rmaj: %f, scale: %f, index: %f, bn: %f, Ro: %f, Io: %f\n", galaxy->Rmaj, scale, index, bn, fR, Io);
+            }
 
             psArrayAdd (galaxies, 100, galaxy);