Changeset 35970 for branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c

Timestamp:

Aug 19, 2013, 6:42:58 AM (13 years ago)

Author:

eugene

Message:

add some sersic support functions

File:

• branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c (modified) (2 diffs)

branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c

@@ -204 +204 @@
 }
 
+float pmSersicNorm (float index) {
+
+    float C0 = NAN;
+    float C1 = NAN;
+    float C2 = NAN;
+
+    // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 
+    // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 
+    // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 
+    // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 
+    // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 
+
+    if ((index >= 0.0) && (index < 1.0)) { 
+        C0 = 0.201545; C1 = -0.950965; C2 = -0.315248;
+        // y = 0.201545 x^0 -0.950965 x^1 -0.315248 x^2 
+    }
+    if ((index >= 1.0) && (index < 2.0)) { 
+        C0 = 0.402084; C1 = -1.357775; C2 = -0.105102;
+        // y = 0.402084 x^0 -1.357775 x^1 -0.105102 x^2 
+    }
+    if ((index >= 2.0) && (index < 3.0)) { 
+        C0 = 0.619093; C1 = -1.591674; C2 = -0.041576;
+        // y = 0.619093 x^0 -1.591674 x^1 -0.041576 x^2 
+    }
+    if ((index >= 3.0) && (index < 4.0)) { 
+        C0 = 0.770263; C1 = -1.696421; C2 = -0.023363;
+        // y = 0.770263 x^0 -1.696421 x^1 -0.023363 x^2 
+    }
+    if ((index >= 4.0) && (index < 5.5)) { 
+        C0 = 0.885891; C1 = -1.755684; C2 = -0.015753; 
+        // y = 0.885891 x^0 -1.755684 x^1 -0.015753 x^2 
+    }
+
+    if (isnan(C0)) return NAN;
+
+    float lnorm = C0 + C1*index + C2*index*index;
+    float norm = exp(lnorm);
+    return norm;
+}
+
 bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing) {
 
@@ -219 +259 @@
     }
 
+    // exponential
+    // f = I0 exp (sqrt(z)) (index = 1.0)
+
+    // f = I0 exp (-0.5*z) 
+
     if (isSersicType) {
         // for a sersic model, 
-        float bn = 1.9992*index - 0.3271;
-        float Io = exp(bn);
+        // 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;
+        // 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;
                     
+        // find the flux of a sersic with peak = 1.0:
+        float norm = pmSersicNorm (index);
+        float flux = 2.0 * M_PI * galaxy->Rmaj * galaxy->Rmin * norm;
+
         // XXX probably should limit the allowed thinness of a galaxy
-        galaxy->peak = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj);
+        galaxy->peak = galaxy->flux / flux;
         return true;
     }