Index: /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertGalaxies.c
===================================================================
--- /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertGalaxies.c	(revision 32158)
+++ /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertGalaxies.c	(revision 32159)
@@ -137,5 +137,5 @@
         axes.minor       = galaxy->Rmin;
         axes.theta       = galaxy->theta;
-	pmPSF_AxesToModel (PAR, axes);
+	pmPSF_AxesToModel (PAR, axes, type);
 	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
 
Index: /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertStars.c
===================================================================
--- /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertStars.c	(revision 32158)
+++ /branches/eam_branches/ipp-20110710/ppSim/src/ppSimInsertStars.c	(revision 32159)
@@ -117,5 +117,5 @@
 	// XXX set the mag & err values (should this be done in pmSourceFromModel?)
 	// XXX i should be applying the gain and the correct effective area
-	psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
+	psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0, model->type);
 	float Area = 2.0 * M_PI * axes.major * axes.minor;
 
Index: /branches/eam_branches/ipp-20110710/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- /branches/eam_branches/ipp-20110710/ppSim/src/ppSimMakeGalaxies.c	(revision 32158)
+++ /branches/eam_branches/ipp-20110710/ppSim/src/ppSimMakeGalaxies.c	(revision 32159)
@@ -20,5 +20,9 @@
     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
@@ -46,4 +50,12 @@
     }
 
+    // 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;
+    }
+
     if (galaxyDensity <= 0) return true;
 
@@ -87,4 +99,6 @@
     long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5;
 
+    fprintf (stderr, "galaxy grid inst mag: %f\n", galaxyGridMag - zp - 2.5*log10(expTime));
+
     if (galaxyGrid) {
 	num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
@@ -98,6 +112,6 @@
 
 	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,7 +120,9 @@
 
 		// 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))
@@ -122,9 +138,4 @@
 		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;
 
@@ -135,13 +146,30 @@
 		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);