Index: /trunk/ppSim/src/ppSimInsertGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 29002)
@@ -1,4 +1,14 @@
 # include "ppSim.h"
 static char *defaultModel = "PS_MODEL_SERSIC";
+
+float imageSum (psImage *image) {
+    float sum = 0.0;
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    sum += image->data.F32[iy][ix];
+	}
+    }
+    return sum;
+}
 
 bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config) {
@@ -84,4 +94,11 @@
     psArray *sources = sources = detections->allSources;
 
+    // output filename
+    char outname[1024];
+    char *outroot = psMetadataLookupStr(&mdok, config->arguments, "OUTPUT");
+    sprintf (outname, "%s.dat", outroot);
+
+    FILE *outfile = fopen (outname, "a");
+
     // add sources to the readout image & weight
     for (long i = 0; i < galaxies->n; i++) {
@@ -131,11 +148,13 @@
 
 	// XXX let the flux limit be a user-defined number of sky sigmas (not just 1.0)
-	float radius = model->modelRadius (model->params, roughNoise);
+	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, 150.0); 
+	radius = PS_MIN (radius, 300.0); 
 
 	// construct a source, with model flux pixels set, based on the model
 	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_EXTENDED);
+
+	galaxy->flux = model->modelFlux (model->params);
 
 	// XXX set the mag & err values (should this be done in pmSourceFromModel?)
@@ -145,44 +164,18 @@
 	source->errMag = sqrt(Area*PS_SQR(roughNoise) + galaxy->flux) / galaxy->flux;	
 	
-	// XXX add the sources to a source array
+	// insert the source flux in the image
+	float sum1 = imageSum(source->pixels);
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	float sum2 = imageSum(source->pixels);
+	float flux = sum2 - sum1;
 
-# if (0)	
-	if (CONVOLVED_FIT) {
+	// fprintf (stderr, "flux: %f, sum1: %f, sum2: %f, flux2: %f, dflux: %f\n", galaxy->flux, sum1, sum2, sum2 - sum1, galaxy->flux - sum2 + sum1);
 
-	    // select the PSF (XXX : move out of loop)
-	    pmPSF *psfModel = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
-	    assert (psf);
-	    
-	    // supply the psf to the source (normalized?)
-	    source->modelPSF = pmModelFromPSFforXY (psfModel, xChip, yChip, 1.0);
+	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->errMag, axes.major, axes.minor, axes.theta, par8);
 
-	    // instantiate the psf flux
-	    pmSourceCachePSF (source, maskVal);
-
-	    // convert the cached cached psf model for this source to a psKernel
-	    // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)
-	    // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)
-	    psKernel *psf = psphotKernelFromPSF (source, psfSize);
-
-	    // instantiate the source model flux
-	    // XXX this may need to handle the offset?
-	    pmSourceCacheModel (source, maskVal);
-
-	    // make a storage buffer for the output image
-	    psImage *buffer = psImageCopy (source->modelFlux, PS_TYPE_F32);
-
-	    // convolve the psf image with the model image
-	    psImageConvolveDirect (buffer, source->modelFlux, psf);
-
-	    // save the result back in the source
-	    psFree (source->modelFlux);
-	    source->modelFlux = buffer;
-	}
-# endif
-
-	// insert the source flux in the image
-	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
 	psArrayAdd (sources, 100,source);
     }
+    fclose (outfile);
 
     // XXX many leaks in here, i think 
@@ -191,3 +184,2 @@
     return true;
 }
-
Index: /trunk/ppSim/src/ppSimInsertStars.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertStars.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimInsertStars.c	(revision 29002)
@@ -105,11 +105,4 @@
         float flux = star->flux * expCorr->data.F32[(int)yCell][(int)xCell];
 
-	// if psfConvolve is TRUE, we will (elsewhere) convolve the image we a PSF
-	// in this case, simply place delta functions in the image
-	if (psfConvolve) {
-	    readout->image->data.F32[(int)(yCell)][(int)(xCell)] += flux;
-	    continue;
-	}
-
 	// instantiate a model for the PSF at this location, set desired flux
 	pmModel *model = pmModelFromPSFforXY (psf, xChip, yChip, 1.0);
@@ -125,26 +118,35 @@
 	// XXX i should be applying the gain and the correct effective area
 	psEllipseAxes axes = pmPSF_ModelToAxes (model->params->data.F32, 20.0);
-	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+	float Area = 2.0 * M_PI * axes.major * axes.minor;
 
-	// this value is the pure (input) flux, and is saved in the output source cmf files
 	source->psfMag = -2.5*log10(star->flux);
 	source->errMag = sqrt(Area*PS_SQR(roughNoise) + flux) / flux;
+	
+	// set the expected model errors
+	model->dparams->data.F32[PM_PAR_I0] = source->errMag * model->params->data.F32[PM_PAR_I0];
 
-	fprintf (outfile, "%8.3f %8.3f  %10.2f  %7.3f %5.3f\n", star->x, star->y, star->flux, source->psfMag, source->errMag);
+	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", star->x, star->y, star->flux, 0, source->psfMag, source->errMag, axes.major, axes.minor, axes.theta, par8);
 
-	// insert the source flux in the image
-	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
-
-	// insert the source flux in the noise image
-	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);
-
-	// Blow away the image parts of the source, which makes the memory explode
-	RESET(source->pixels);
-	RESET(source->variance);
-	RESET(source->maskObj);
-	RESET(source->maskView);
-	RESET(source->modelFlux);
-	RESET(source->psfImage);
-	RESET(source->blends);
+	// if psfConvolve is TRUE, we will (elsewhere) convolve the image we a PSF
+	// in this case, simply place delta functions in the image
+	if (psfConvolve) {
+	    readout->image->data.F32[(int)(yCell)][(int)(xCell)] += flux;
+	} else {
+	    // insert the source flux in the image
+	    pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	    
+	    // insert the source flux in the noise image
+	    pmSourceAddWithOffset (source, PM_MODEL_OP_FULL | PM_MODEL_OP_NOISE, 0xff, dX, dY);
+	    
+	    // Blow away the image parts of the source, which makes the memory explode
+	    RESET(source->pixels);
+	    RESET(source->variance);
+	    RESET(source->maskObj);
+	    RESET(source->maskView);
+	    RESET(source->modelFlux);
+	    RESET(source->psfImage);
+	    RESET(source->blends);
+	}
 
 	// add the sources to the source array
Index: /trunk/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 29001)
+++ /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);
Index: /trunk/ppSim/src/ppSimPhotomReadout.c
===================================================================
--- /trunk/ppSim/src/ppSimPhotomReadout.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimPhotomReadout.c	(revision 29002)
@@ -69,5 +69,5 @@
     psArray *fakeSources = psArrayAlloc (injectedSources->n);
     for (int i = 0; i < injectedSources->n; i++) {
-      fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]);
+      fakeSources->data[i] = pmSourceCopyData (injectedSources->data[i]);
     }
 
Index: /trunk/ppSim/src/ppSimPhotomReadoutFake.c
===================================================================
--- /trunk/ppSim/src/ppSimPhotomReadoutFake.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimPhotomReadoutFake.c	(revision 29002)
@@ -40,5 +40,5 @@
     psArray *fakeSources = psArrayAlloc (injectedSources->n);
     for (int i = 0; i < injectedSources->n; i++) {
-      fakeSources->data[i] = pmSourceCopy (injectedSources->data[i]);
+      fakeSources->data[i] = pmSourceCopyData (injectedSources->data[i]);
     }
 
@@ -47,5 +47,5 @@
     psArray *realSources = psArrayAlloc (realMeasuredSources->n);
     for (int i = 0; i < realMeasuredSources->n; i++) {
-        realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
+        realSources->data[i] = pmSourceCopyData (realMeasuredSources->data[i]);
     }
 
Index: /trunk/ppSim/src/ppSimPhotomReadoutForce.c
===================================================================
--- /trunk/ppSim/src/ppSimPhotomReadoutForce.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimPhotomReadoutForce.c	(revision 29002)
@@ -46,5 +46,5 @@
     psArray *realSources = psArrayAlloc (realMeasuredSources->n);
     for (int i = 0; i < realMeasuredSources->n; i++) {
-        realSources->data[i] = pmSourceCopy (realMeasuredSources->data[i]);
+        realSources->data[i] = pmSourceCopyData (realMeasuredSources->data[i]);
     }
 
Index: /trunk/ppSim/src/ppSimSmoothReadout.c
===================================================================
--- /trunk/ppSim/src/ppSimSmoothReadout.c	(revision 29001)
+++ /trunk/ppSim/src/ppSimSmoothReadout.c	(revision 29002)
@@ -7,5 +7,5 @@
     // XXX use these defaults?
     // float minGauss = 0.1;
-    float nSigma = 3.0;
+    float nSigma = 4.0;
     float sigma = psMetadataLookupF32(&status, recipe, "SEEING"); // Seeing SIGMA (pixels)
 
@@ -28,2 +28,6 @@
 // 4*pi*SIGMA_SMTH^2, but for measurements based on apertures comparable to or larger than the
 // smoothing kernel, the effective per-pixel variance is maintained.
+
+
+// XXX this only allows for Gaussian PSFs.  To extend this, we need to convolve with an image
+// of the PSF model