Index: /trunk/ppSim/src/ppSimArguments.c
===================================================================
--- /trunk/ppSim/src/ppSimArguments.c	(revision 34260)
+++ /trunk/ppSim/src/ppSimArguments.c	(revision 34261)
@@ -225,5 +225,5 @@
         ppSimArgToRecipeStr(&status, options, "GALAXY.MODEL", arguments, "-galmodel"); // Galaxy model name
 
-	// the user supplies FWHM in arcsec; here we convert to Sigma in pixels
+	// the user supplies FWHM in arcsec; we convert to Sigma in pixels (in ppSimCreate.c:103)
         psMetadataAddF32(options, PS_LIST_TAIL, "SEEING", 0, "Seeing FWHM (arcsec)", seeing);
         psMetadataAddF32(options, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)", ra0 * M_PI / 180.0);
Index: /trunk/ppSim/src/ppSimInsertGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 34260)
+++ /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 34261)
@@ -76,5 +76,5 @@
     if (type == -1) {
 	psError (PS_ERR_UNKNOWN, false, "invalid model name");
-        return false;
+	return false;
     }
 
@@ -136,10 +136,16 @@
 	PAR[PM_PAR_YPOS] = yChip;
 
-	psEllipseAxes axes;
-        axes.major       = galaxy->Rmaj;
-        axes.minor       = galaxy->Rmin;
-        axes.theta       = galaxy->theta;
-	pmPSF_AxesToModel (PAR, axes, type);
-	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+	if (type == pmModelClassGetType ("PS_MODEL_TRAIL")) {
+	    PAR[PM_PAR_LENGTH] = galaxy->Rmaj;
+	    PAR[PM_PAR_SIGMA]  = galaxy->Rmin;
+	    PAR[PM_PAR_THETA]  = galaxy->theta;
+	} else {
+	    psEllipseAxes axes;
+	    axes.major       = galaxy->Rmaj;
+	    axes.minor       = galaxy->Rmin;
+	    axes.theta       = galaxy->theta;
+	    pmPSF_AxesToModel (PAR, axes, type);
+	}
+	psF64 Area = 2.0 * M_PI * galaxy->Rmaj * galaxy->Rmin;
 
 	if (nParam == 8) {
@@ -147,10 +153,8 @@
 	}
 
-	// set the normalization to get the desired flux
-	// XXX for now, use peak: pmModelSetFlux (model, galaxy->flux);
-
 	// XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
 	float radius = model->modelRadius (model->params, 0.1*roughNoise);
 	radius = PS_MAX (radius, 1.0);
+
 	// XXX the exp(-r^0.25) models can go way out if allowed...
 	radius = PS_MIN (radius, 300.0); 
@@ -176,5 +180,5 @@
 
 	float par8 = (model->params->n == 8) ? model->params->data.F32[7] : 0.0;
-	fprintf (outfile, "%8.3f %8.3f %10.2f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->psfMagErr, axes.major, axes.minor, axes.theta, par8);
+	fprintf (outfile, "%8.3f %8.3f %10.2f  %2d  %7.3f %5.3f  %5.3f %5.3f %5.3f %5.3f\n", galaxy->x, galaxy->y, flux, 1, source->psfMag, source->psfMagErr, galaxy->Rmaj, galaxy->Rmin, galaxy->theta, par8);
 
 	psArrayAdd (sources, 100,source);
Index: /trunk/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 34260)
+++ /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;
+}