Index: trunk/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 32350)
+++ trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 34261)
@@ -1,3 +1,5 @@
 # include "ppSim.h"
+
+bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing);
 
 bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
@@ -19,5 +21,5 @@
     int galaxyGridDY 	  = psMetadataLookupS32 (&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
     bool galaxyGridRandom = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID.RANDOM"); // Density of fakes
-    float galaxyGridPeak  = psMetadataLookupF32 (&mdok, recipe, "GALAXY.GRID.PEAK"); // peak flux 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
     
@@ -42,5 +44,5 @@
     float expTime  	  = psMetadataLookupF32(&mdok, recipe, "EXPTIME"); // Exposure time
     float zp       	  = psMetadataLookupF32(&mdok, recipe, "ZEROPOINT"); // Photometric zero point
-    // float seeing   	  = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix)
+    float seeing   	  = psMetadataLookupF32(&mdok, recipe, "SEEING"); // Seeing sigma (pix)
     float scale    	  = psMetadataLookupF32(&mdok, recipe, "PIXEL.SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
     float skyRate  	  = psMetadataLookupF32(&mdok, recipe, "SKY.RATE"); // Sky rate
@@ -150,26 +152,6 @@
 		galaxy->theta = (galaxyThetaMin  + rndValue * galaxyThetaSlope);
 		
-		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;
-		}
+		galaxy->flux = galaxyGridFlux;
+		ppSimSetGalaxyPeak (galaxy, type, index, seeing);
 
 		// if (0) {
@@ -206,27 +188,9 @@
 	    galaxy->theta = (psRandomUniform(rng) * galaxyThetaSlope  + galaxyThetaMin);
 
-	    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
-		float galaxyMag =  PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / galaxyAlpha)));
-
-		galaxy->flux = ppSimMagToFlux (galaxyMag, zp) * expTime;
-		galaxy->peak = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj);
-	    }
+	    // XXX probably should limit the allowed thinness of a galaxy
+	    float galaxyMag =  PS_MIN (faintMag, PS_MAX (brightMag, (log10((i + 1.0) / normScale) / galaxyAlpha)));
+	    
+	    galaxy->flux = ppSimMagToFlux (galaxyMag, zp) * expTime;
+	    ppSimSetGalaxyPeak (galaxy, type, index, seeing);
 
 	    if (1) {
@@ -239,2 +203,44 @@
     return true;
 }
+
+bool ppSimSetGalaxyPeak (ppSimGalaxy *galaxy, pmModelType type, float index, float seeing) {
+
+    bool isSersicType = false;
+    if (type == pmModelClassGetType ("PS_MODEL_DEV")) {
+	index = 4.0;
+	isSersicType = true;
+    }
+    if (type == pmModelClassGetType ("PS_MODEL_EXP")) {
+	index = 1.0;
+	isSersicType = true;
+    }
+    if (type == pmModelClassGetType ("PS_MODEL_SERSIC")) {
+	isSersicType = true;
+    }
+
+    if (isSersicType) {
+	// 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 = galaxy->flux / norm / (galaxy->Rmin / galaxy->Rmaj);
+	return true;
+    }
+
+    if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) {
+	// galaxy->Rmaj -> PM_PAR_LENGTH
+	// seeing -> PM_PAR_SIGMA
+	galaxy->Rmin = seeing; // Rmin is used for sigma
+	galaxy->peak = galaxy->flux / (galaxy->Rmaj * galaxy->Rmin * 2.0 * sqrt(2.0 * M_PI));
+	return true;
+    }
+
+    galaxy->peak = galaxy->flux / (2 * M_PI * galaxy->Rmaj * galaxy->Rmin);
+    return true;
+}