Index: /branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- /branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c	(revision 35969)
+++ /branches/eam_branches/ipp-20130711/ppSim/src/ppSimMakeGalaxies.c	(revision 35970)
@@ -204,4 +204,44 @@
 }
 
+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,16 +259,25 @@
     }
 
+    // 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;
     }