Index: trunk/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 18011)
+++ trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 29002)
@@ -18,4 +18,5 @@
     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,4 +44,11 @@
 
     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,8 +89,8 @@
 	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,10 +104,36 @@
 		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 factor = galaxyGridRandom ? drand48() : i / num;
+		float index = (galaxyIndexMin  + factor * galaxyIndexSlope); // factor of 0.5 because the Sersic model creates exp(-z^n), not exp(-r^n)
+		galaxy->index = 0.5/index;
+
+		factor = galaxyGridRandom ? drand48() : i / num;
+		float scale = (galaxyRmajorMin + factor * 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 * fR;
+
+		galaxy->Rmaj  = scale;
+
+		factor = galaxyGridRandom ? drand48() : i / num;
+		galaxy->Rmin  = (galaxyARatioMin + factor * galaxyARatioSlope) * galaxy->Rmaj;
+
+		factor = galaxyGridRandom ? drand48() : i / num;
+		galaxy->theta = (galaxyThetaMin  + factor * 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,14 +158,15 @@
 	    galaxy->y    = psRandomUniform(rng) * ySize; // y position
 
-	    galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum));
+	    // galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum));
 
-	    galaxy->index = (psRandomUniform(rng) * galaxyIndexSlope  + galaxyIndexMin);
+	    galaxy->index = (sqrt(psRandomUniform(rng)) * galaxyIndexSlope  + galaxyIndexMin);
 	    galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope  + galaxyThetaMin);
 	    galaxy->Rmaj  = (psRandomUniform(rng) * galaxyRmajorSlope + galaxyRmajorMin);
-	    galaxy->Rmin  = (psRandomUniform(rng) * galaxyARatioSlope + galaxyARatioMin) * galaxy->Rmaj;
+	    galaxy->Rmin  = (PS_SQR(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;
 
 	    psArrayAdd (galaxies, 100, galaxy);