Index: /trunk/ppSim/src/Makefile.am
===================================================================
--- /trunk/ppSim/src/Makefile.am	(revision 14666)
+++ /trunk/ppSim/src/Makefile.am	(revision 14667)
@@ -9,4 +9,5 @@
 	ppSimLoadStars.c        \
 	ppSimMakeStars.c        \
+	ppSimMakeGalaxies.c     \
 	ppSimMakeBiassec.c      \
 	ppSimMakeBias.c         \
@@ -14,4 +15,5 @@
 	ppSimMakeSky.c          \
 	ppSimInsertSources.c    \
+	ppSimInsertGalaxies.c   \
 	ppSimAddOverscan.c      \
 	ppSimAddNoise.c         \
Index: /trunk/ppSim/src/ppSim.h
===================================================================
--- /trunk/ppSim/src/ppSim.h	(revision 14666)
+++ /trunk/ppSim/src/ppSim.h	(revision 14667)
@@ -48,6 +48,21 @@
 } ppSimStar;
 
+typedef struct {
+    double ra;
+    double dec;
+    float mag;
+    float x;
+    float y;
+    float flux;
+
+    float peak;
+    float Rmaj;
+    float Rmin;
+    float theta;
+    float index;
+} ppSimGalaxy;
 
 ppSimStar *ppSimStarAlloc ();
+ppSimGalaxy *ppSimGalaxyAlloc ();
 
 /// Parse command-line arguments
@@ -85,5 +100,5 @@
 		     );
 
-psArray *ppSimLoadStars (pmFPA *fpa, pmConfig *config);
+bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config);
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng);
 
@@ -117,3 +132,6 @@
 bool ppSimSetPSF (pmChip *chip, pmConfig *config);
 
+bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng);
+bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config);
+
 #endif
Index: /trunk/ppSim/src/ppSimArguments.c
===================================================================
--- /trunk/ppSim/src/ppSimArguments.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimArguments.c	(revision 14667)
@@ -71,6 +71,8 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist");
 
-    if (!config->camera) {
-        psErrorStackPrint(stderr, "A camera name must be specified using the -camera option.");
+    // only one of -camera and -file is needed
+    bool status = pmConfigFileSetsMD (config->arguments, &argc, argv, "INPUT", "-file", "-list");
+    if (!config->camera && !status) {
+        psErrorStackPrint(stderr, "A camera name (-camera NAME) or an image (-file NAME) must be specified");
         usage(argv[0], arguments, config);
     }
@@ -82,4 +84,5 @@
     psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
     if (formatName) {
+	// XXX delay the config below until ppSimCreate?
         config->formatName = psMemIncrRefCounter(formatName);
 
Index: /trunk/ppSim/src/ppSimCreate.c
===================================================================
--- /trunk/ppSim/src/ppSimCreate.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimCreate.c	(revision 14667)
@@ -9,10 +9,26 @@
 pmFPAfile *ppSimCreate(pmConfig *config)
 {
+    bool status;
+    bool simImage = false;
     PS_ASSERT_PTR_NON_NULL(config, NULL);
+    pmFPA *fpa = NULL;
 
-    pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the observation
-    if (!fpa) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
-        return NULL;
+    // the input image defines the camera.  if it is not supplied, the user must have
+    // supplied a camera and other metadata on the command line
+    pmFPAfile *input = pmFPAfileDefineFromArgs (&status, config, "PPIMAGE.INPUT", "INPUT");
+    if (!input) {
+	simImage = true;
+	fpa = pmFPAConstruct(config->camera); // FPA to contain the observation
+	if (!fpa) {
+	    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to construct an FPA from camera configuration.");
+	    return NULL;
+	}
+    } else {
+	simImage = false;
+	if (input->type != PM_FPA_FILE_IMAGE) {
+	    psError(PS_ERR_IO, true, "PPIMAGE.INPUT is not of type IMAGE");
+	    return NULL;
+	}
+	fpa = input->fpa;
     }
 
@@ -56,4 +72,9 @@
     }
     simSources->save = true;
+
+    // if we have loaded an input image, we do not need to populate the fpa
+    if (!simImage) {
+	return file;
+    }
 
     pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes
Index: /trunk/ppSim/src/ppSimInsertGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 14667)
+++ /trunk/ppSim/src/ppSimInsertGalaxies.c	(revision 14667)
@@ -0,0 +1,150 @@
+# include "ppSim.h"
+static char *defaultModel = "PS_MODEL_SERSIC";
+
+bool ppSimInsertGalaxies (pmReadout *readout, psImage *expCorr, psArray *galaxies, pmConfig *config) {
+
+    bool mdok;
+
+    assert (readout);
+    assert (galaxies);
+    
+    if (!galaxies->n) { return true; }
+
+    // XXX is this needed?
+    // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
+
+    pmCell *cell = readout->parent;
+    pmChip *chip = cell->parent;
+
+    // XXX this is an estimate of the sky noise based on the inputs to the image simulation.
+    // XXX update this to allow the estimate based on the measured sky background
+    // XXX this is missing the gain.
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+
+    float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
+    float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
+    float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
+    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
+    if (isnan(readnoise)) {
+	psWarning("CELL.READNOISE is not set; reverting to recipe value READNOISE.");
+	readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE");
+	if (!mdok) {
+	    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+	    return false;
+	}
+    }
+
+    // Rough noise estimate, appropriate for entire cell (use for source radius?)
+    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
+    int x0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.X0");
+    int y0Chip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.Y0");
+    int xParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.XPARITY");
+    int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY");
+
+    int x0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
+    int y0Cell = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+    int xParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.XPARITY");
+    int yParityCell = psMetadataLookupS32(NULL, cell->concepts, "CELL.YPARITY");
+
+    int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
+
+    // determine the galaxy model
+    char *modelName = psMetadataLookupStr(&mdok, config->arguments, "GALAXY.MODEL"); // Seeing sigma (pix)
+    if (modelName == NULL) {
+	modelName = defaultModel;
+    }
+    pmModelType type = pmModelClassGetType (modelName);
+    if (type == -1) {
+	psError (PS_ERR_UNKNOWN, false, "invalid model name");
+        return false;
+    }
+
+    int nParam = pmModelClassParameterCount (type);
+
+    pmPSF *psf = psMetadataLookupPtr (&mdok, chip->analysis, "PSPHOT.PSF");
+    assert (psf);
+
+    int dX = PM_CELL_TO_CHIP (0.0, x0Cell, xParityCell, binning);
+    int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning);
+
+    // psMetadataLookupPtr (readout->analysis, "PSPHOT.SOURCES", 0);
+
+    psArray *sources = psArrayAllocEmpty (galaxies->n);
+
+
+    // add sources to the readout image & weight
+    for (long i = 0; i < galaxies->n; i++) {
+	ppSimGalaxy *galaxy = galaxies->data[i];
+
+	// galaxy->x,y are in fpa coordinates
+
+	// Position on the cell and peak flux
+	float xChip = PM_FPA_TO_CHIP(galaxy->x, x0Chip, xParityChip);
+	float yChip = PM_FPA_TO_CHIP(galaxy->y, y0Chip, yParityChip);
+
+	// Position on the cell and peak flux
+	float xCell = PM_CHIP_TO_CELL(xChip, x0Cell, xParityCell, binning);
+	float yCell = PM_CHIP_TO_CELL(yChip, y0Cell, yParityCell, binning);
+
+	if (xCell < 0) continue;
+	if (yCell < 0) continue;
+	if (xCell > readout->image->numCols) continue;
+	if (yCell > readout->image->numRows) continue;
+	// XXX need to apply col0, row0 if readout is a subarray
+
+	// XXX apply the expCorr to the galaxy->flux before setting the model flux
+
+	// instantiate a model for the PSF at this location, set desired flux
+	pmModel *model = pmModelAlloc (type);
+	
+	psF32 *PAR = model->params->data.F32;
+
+	PAR[PM_PAR_I0]   = galaxy->peak;
+	PAR[PM_PAR_SKY]  = 0.0;
+	PAR[PM_PAR_XPOS] = xChip;
+	PAR[PM_PAR_YPOS] = yChip;
+
+	psEllipseAxes axes;
+        axes.major       = galaxy->Rmaj;
+        axes.minor       = galaxy->Rmin;
+        axes.theta       = galaxy->theta;
+	pmPSF_AxesToModel (PAR, axes);
+	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+
+	if (nParam == 8) {
+	    PAR[PM_PAR_7] = galaxy->index;
+	}
+
+	// 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, 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); 
+
+	// construct a source, with model flux pixels set, based on the model
+	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_EXTENDED);
+
+	// XXX set the mag & err values (should this be done in pmSourceFromModel?)
+	// XXX i should be applying the gain and the correct effective area
+	source->psfMag = -2.5*log10(galaxy->flux);
+	source->extMag = -2.5*log10(galaxy->flux);
+	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
+	pmSourceAddWithOffset (source, PM_MODEL_OP_FULL, 0xff, dX, dY);
+	psArrayAdd (sources, 100,source);
+    }
+
+    // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT"
+    // psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sources);
+
+    // XXX many leaks in here, i think 
+    return true;
+}
+
Index: /trunk/ppSim/src/ppSimInsertSources.c
===================================================================
--- /trunk/ppSim/src/ppSimInsertSources.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimInsertSources.c	(revision 14667)
@@ -5,7 +5,9 @@
     bool mdok;
 
+    assert (readout);
     assert (stars);
-    assert (readout);
-    
+
+    if (!stars->n) { return true; }
+
     // XXX is this needed?
     // pmFPAfile *simSources = psMetadataLookupPtr(NULL, config->files, "PPSIM.SOURCES"); // Output sources
@@ -88,4 +90,12 @@
 	pmSource *source = pmSourceFromModel (model, readout, radius, PM_SOURCE_TYPE_STAR);
 
+	// 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);
+	psF64 Area = 2.0 * M_PI * axes.major * axes.minor;
+
+	source->psfMag = -2.5*log10(star->flux);
+	source->errMag = sqrt(Area*PS_SQR(roughNoise) + star->flux) / star->flux;	
+	
 	// XXX add the sources to a source array
 
Index: /trunk/ppSim/src/ppSimLoadStars.c
===================================================================
--- /trunk/ppSim/src/ppSimLoadStars.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimLoadStars.c	(revision 14667)
@@ -1,5 +1,13 @@
 #include "ppSim.h"
 
-psArray *ppSimLoadStars (pmFPA *fpa, pmConfig *config) {
+bool ppSimLoadStars (psArray *stars, pmFPA *fpa, pmConfig *config) {
+
+    bool mdok;
+    assert (stars);
+
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
+
+    bool starsReal = psMetadataLookupBool(&mdok, recipe, "STARS.REAL"); // Density of fakes
+    if (!starsReal) return true;
 
     // Read catalogue stars using psastro
@@ -38,5 +46,6 @@
     psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld reference stars", refStars->n);
 
-    psArray *stars = psArrayAlloc (refStars->n);
+    long oldSize = stars->n;
+    stars = psArrayRealloc (stars, refStars->n);
 
     // Conversion loop
@@ -61,5 +70,5 @@
 	star->flux = powf(10.0, -0.4 * (star->mag - zp)) * expTime;
 	star->peak = star->flux / sqrt(2.0*M_PI) / seeing;
-	stars->data[i] = star;
+	stars->data[oldSize + i] = star;
     }
     psFree(refStars);
Index: /trunk/ppSim/src/ppSimLoop.c
===================================================================
--- /trunk/ppSim/src/ppSimLoop.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimLoop.c	(revision 14667)
@@ -17,7 +17,7 @@
 
     // Load catalogue stars
-    psArray *stars = NULL;
+    psArray *stars = psArrayAllocEmpty (1);
     if (type == PPSIM_TYPE_OBJECT) {
-	stars = ppSimLoadStars (fpa, config);
+	ppSimLoadStars (stars, fpa, config);
     }
 
@@ -25,4 +25,10 @@
     if (type == PPSIM_TYPE_OBJECT) {
 	ppSimMakeStars (stars, fpa, config, rng);
+    }
+
+    // Add random galaxies
+    psArray *galaxies = psArrayAllocEmpty (1);
+    if (type == PPSIM_TYPE_OBJECT) {
+	ppSimMakeGalaxies (galaxies, fpa, config, rng);
     }
 
@@ -62,8 +68,5 @@
         while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
 
-            // Size, position and orientation of cell
-            int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
-            int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
-
+	    // check that we are able to work with this cell (readdir must be 1)
             int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction
             if (readdir != 1) {
@@ -73,8 +76,12 @@
             }
 
-	    psVector *biasCols = ppSimMakeBiassec (cell, config);
+	    // if no readout exists, generate a single readout (upgrade to allow multiple
+	    // readouts?)
+	    if (cell->readouts->n == 0) {
+		// Size, position and orientation of cell
+		int numCols = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE") / binning;
+		int numRows = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE") / binning;
 
-	    int numReadouts = 1;
-            for (int i = 0; i < numReadouts; i++) {
+		// generate a new readout for this cell
                 pmReadout *readout = pmReadoutAlloc(cell); // Readout within cell
 
@@ -83,12 +90,32 @@
                 readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
 
+		psImageInit (readout->image, 0.0);
+		psImageInit (readout->weight, 0.0);
+	    }
+
+	    psVector *biasCols = ppSimMakeBiassec (cell, config);
+
+            for (int i = 0; i < cell->readouts->n; i++) {
+
+		pmReadout *readout = cell->readouts->data[i];
+		assert (readout);
+
+		// if we have not read in a weight or generated a fake image above, we need to
+		// build one here
+		if (!readout->weight) {
+		    if (!pmReadoutGenerateWeight(readout, true)) {
+			psError (PS_ERR_UNKNOWN, false, "trouble creating weight");
+			return false;
+		    }
+		}
+
                 psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
                 if (type == PPSIM_TYPE_OBJECT) {
-                    expCorr = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                    expCorr = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_F32);
                 }
 
 		psVector *biasRows = ppSimMakeBias (readout, config, rng);
 		if (type == PPSIM_TYPE_BIAS) goto done;
-		
+	
 		ppSimMakeDark (readout, config);
 		if (type == PPSIM_TYPE_DARK) goto done;
@@ -99,4 +126,8 @@
 		if (type == PPSIM_TYPE_OBJECT) {
 		    ppSimInsertSources (readout, expCorr, stars, config);
+		}
+
+		if (type == PPSIM_TYPE_OBJECT) {
+		    ppSimInsertGalaxies (readout, expCorr, galaxies, config);
 		}
 
@@ -111,5 +142,4 @@
                 readout->parent->data_exists = true;
                 readout->parent->parent->data_exists = true;
-                psFree(readout);        // Drop reference
             }
             psFree(biasCols);
Index: /trunk/ppSim/src/ppSimMakeBias.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeBias.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimMakeBias.c	(revision 14667)
@@ -46,6 +46,6 @@
 
 	    // Bias level
-	    signal->data.F32[y][x] = biasRows->data.F32[y];
-	    variance->data.F32[y][x] = PS_SQR(readnoise);
+	    signal->data.F32[y][x] += biasRows->data.F32[y];
+	    variance->data.F32[y][x] += PS_SQR(readnoise);
 
 	}
Index: /trunk/ppSim/src/ppSimMakeGalaxies.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 14667)
+++ /trunk/ppSim/src/ppSimMakeGalaxies.c	(revision 14667)
@@ -0,0 +1,139 @@
+# include "ppSim.h"
+
+bool ppSimMakeGalaxies(psArray *galaxies, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
+
+    bool mdok;
+    assert (galaxies);
+
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
+
+    bool galaxyFake = psMetadataLookupBool(&mdok, recipe, "GALAXY.FAKE"); // Density of fakes
+    if (!galaxyFake) return NULL;
+
+    float galaxyLum = psMetadataLookupF32(&mdok, recipe, "GALAXY.LUM"); // Galaxy luminosity func slope
+    float galaxyMag = psMetadataLookupF32(&mdok, recipe, "GALAXY.MAG"); // Galaxy brightest magnitude
+    float galaxyDensity = psMetadataLookupF32(&mdok, recipe, "GALAXY.DENSITY"); // Density of fakes
+
+    bool galaxyGrid = psMetadataLookupBool(&mdok, recipe, "GALAXY.GRID"); // Density of fakes
+    int galaxyGridDX = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DX"); // Density of fakes
+    int galaxyGridDY = psMetadataLookupS32(&mdok, recipe, "GALAXY.GRID.DY"); // Density of fakes
+    
+    float galaxyRmajorMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MAX"); // Density of fakes
+    float galaxyRmajorMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.RMAJOR.MIN"); // Density of fakes
+
+    float galaxyARatioMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MAX"); // Density of fakes
+    float galaxyARatioMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.ARATIO.MIN"); // Density of fakes
+
+    float galaxyThetaMax = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MAX"); // Density of fakes
+    float galaxyThetaMin = psMetadataLookupF32(&mdok, recipe, "GALAXY.THETA.MIN"); // Density of fakes
+
+    float galaxyIndexMin  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MIN"); // Density of fakes
+    float galaxyIndexMax  = psMetadataLookupF32(&mdok, recipe, "GALAXY.INDEX.MAX"); // Density of fakes
+
+    // XXX push these into the recipe...
+    float skyRate = psMetadataLookupF32(&mdok, config->arguments, "SKY.RATE"); // Sky rate
+    float darkRate = psMetadataLookupF32(&mdok, config->arguments, "DARK.RATE"); // Dark rate
+    float expTime = psMetadataLookupF32(&mdok, config->arguments, "EXPTIME"); // Exposure time
+
+    float zp = psMetadataLookupF32(&mdok, config->arguments, "ZEROPOINT"); // Photometric zero point
+    float seeing = psMetadataLookupF32(&mdok, config->arguments, "SEEING"); // Seeing sigma (pix)
+    float scale = psMetadataLookupF32(&mdok, config->arguments, "SCALE") * M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
+
+    if (galaxyDensity <= 0) return NULL;
+
+    // Size of FPA
+    psRegion *bounds = ppSimFPABounds (fpa);
+    // Size of focal plane (can't I get this from the config?)
+    int xSize = bounds->x1 - bounds->x0;
+    int ySize = bounds->y1 - bounds->y0;
+    psFree(bounds);
+
+    // Grabbing read noise from the recipe rather than the cell, which is a potential danger, but it
+    // shouldn't be too bad.
+    float readnoise = psMetadataLookupF32(&mdok, recipe, "READNOISE"); // Default read noise
+    if (!mdok) {
+	psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find READNOISE in recipe.");
+	psFree(bounds);
+	return false;
+    }
+
+    // Peak fluxes: faintest and brightest levels for random galaxy
+    float faint = 0.1 * 10.0 * sqrt(2.0*M_PI) * seeing * sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+    float bright = powf(10.0, -0.4 * (galaxyMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
+    if (bright < faint) {
+	psLogMsg("ppSim", PS_LOG_INFO,
+		 "Image noise is above brightest random galaxy --- no random galaxy added.");
+	return NULL;
+    }
+
+    // Normalisation, set by the specified stellar density at the specified bright magnitude
+    float norm = galaxyDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
+	powf(bright, galaxyLum);
+
+    // Total number of galaxy down to the faint flux end
+    long num = expf(logf(norm) + galaxyLum * logf(faint)) + 0.5;
+
+    if (galaxyGrid) {
+	num = (int)(xSize / galaxyGridDX) * (int)(ySize / galaxyGridDY);
+    
+	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;
+
+	int i = 0;
+
+	for (long iy = 0.5*galaxyGridDY; iy < ySize; iy += galaxyGridDY) {
+	    for (long ix = 0.5*galaxyGridDX; ix < xSize; ix += galaxyGridDX) {
+		ppSimGalaxy *galaxy = ppSimGalaxyAlloc ();
+
+		// make fpa center of distribution
+		galaxy->x    = ix;
+		galaxy->y    = iy;
+
+		galaxy->peak = 20000;
+
+		galaxy->index = (galaxyIndexMin  + i * galaxyIndexSlope);
+		galaxy->Rmaj  = (galaxyRmajorMin + i * galaxyRmajorSlope);
+		galaxy->Rmin  = (galaxyARatioMin + i * galaxyARatioSlope) * galaxy->Rmaj;
+		galaxy->theta = (galaxyThetaMin  + i * galaxyThetaSlope);
+
+		psArrayAdd (galaxies, 100, galaxy);
+		i++;
+	    }
+	}
+    } else {    
+
+	float galaxyIndexSlope = (galaxyIndexMax - galaxyIndexMin);
+	float galaxyThetaSlope = (galaxyThetaMax - galaxyThetaMin);
+	float galaxyARatioSlope = (galaxyARatioMax - galaxyARatioMin);
+	float galaxyRmajorSlope = (galaxyRmajorMax - galaxyRmajorMin);
+
+	psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld galaxies down to %f mag\n", num,
+		 -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
+
+	for (long i = 0; i < num; i++) {
+	    ppSimGalaxy *galaxy = ppSimGalaxyAlloc ();
+
+	    // make fpa center of distribution
+	    galaxy->x    = psRandomUniform(rng) * xSize; // x position
+	    galaxy->y    = psRandomUniform(rng) * ySize; // y position
+
+	    galaxy->flux = pow ((double)((i + 1) / norm), (double) (1.0 / galaxyLum));
+
+	    galaxy->index = (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->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);
+	}
+    }
+    return galaxies;
+}
Index: /trunk/ppSim/src/ppSimMakeStars.c
===================================================================
--- /trunk/ppSim/src/ppSimMakeStars.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimMakeStars.c	(revision 14667)
@@ -3,9 +3,11 @@
 bool ppSimMakeStars(psArray *stars, pmFPA *fpa, pmConfig *config, const psRandom *rng) {
 
+    bool mdok;
     assert (stars);
 
-    bool mdok;
+    psMetadata *recipe = psMetadataLookupMetadata(&mdok, config->recipes, PPSIM_RECIPE); // Recipe
 
-    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+    bool starsFake = psMetadataLookupBool(&mdok, recipe, "STARS.FAKE"); // Density of fakes
+    if (!starsFake) return true;
 
     float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
@@ -45,31 +47,34 @@
 	psLogMsg("ppSim", PS_LOG_INFO,
 		 "Image noise is above brightest random star --- no random stars added.");
-    } else {
+	return true;
+    }
 
-	// Normalisation, set by the specified stellar density at the specified bright magnitude
-	float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
-	    powf(bright, starsLum);
+    // Normalisation, set by the specified stellar density at the specified bright magnitude
+    float norm = starsDensity * xSize * ySize * PS_SQR(scale * 180.0 / M_PI) /
+	powf(bright, starsLum);
 
-	// Total number of stars down to the faint flux end
-	long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
+    // Total number of stars down to the faint flux end
+    long num = expf(logf(norm) + starsLum * logf(faint)) + 0.5;
 
-	psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
-		 -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
+    psLogMsg("ppSim", PS_LOG_INFO, "Adding %ld stars down to %f mag\n", num,
+	     -2.5 * log10(faint * sqrt(2.0*M_PI) * seeing / expTime) + zp);
 
-	long oldSize = stars->n;
-	psArrayRealloc (stars, stars->n + num);
+    long oldSize = stars->n;
+    psArrayRealloc (stars, stars->n + num);
 
-	for (long i = 0; i < num; i++) {
-	    ppSimStar *star = ppSimStarAlloc ();
+    for (long i = 0; i < num; i++) {
+	ppSimStar *star = ppSimStarAlloc ();
 
-	    // make fpa center of distribution
-	    star->x    = psRandomUniform(rng) * xSize; // x position
-	    star->y    = psRandomUniform(rng) * ySize; // y position
-	    star->peak = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
-	    star->flux = star->peak * sqrt(2.0*M_PI)*seeing;
-	    stars->data[oldSize + i] = star;
-	}
-	stars->n = stars->nalloc;
+	// make fpa center of distribution
+	star->x    = psRandomUniform(rng) * xSize; // x position
+	star->y    = psRandomUniform(rng) * ySize; // y position
+
+	star->flux = expf((logf(i + 1) - logf(norm)) / starsLum); // Peak flux
+	star->peak = star->flux / (2.0*M_PI*PS_SQR(seeing));
+
+	stars->data[oldSize + i] = star;
     }
+    stars->n = stars->nalloc;
+
     return true;
 }
Index: /trunk/ppSim/src/ppSimStars.c
===================================================================
--- /trunk/ppSim/src/ppSimStars.c	(revision 14666)
+++ /trunk/ppSim/src/ppSimStars.c	(revision 14667)
@@ -1,64 +1,3 @@
 # include "ppSim.h"
-
-#define STAR_RANGE_MIN 4.5              // Minimum pixel range for adding stars, sigma
-#define STAR_BG_FRAC 0.1                // Fraction of background error to go out to when adding stars
-
-// Return star flux at a position
-static inline float starFlux(float x, float y, // Position of interest
-                             float x0, float y0, // Centre of star
-                             float seeing // Seeing FWHM (pixels)
-    )
-{
-#ifdef STAR_GAUSSIAN
-    // Gaussian star
-    return exp(-0.5 * (PS_SQR(x - x0) + PS_SQR(y - y0)) / PS_SQR(seeing));
-#else
-    // Waussian star
-    float z = (PS_SQR(x - x0) + PS_SQR(y - y0)) / (2.0 * PS_SQR(seeing));
-    return 1.0 / (1.0 + z + PS_SQR(z));
-#endif
-}
-
-// Return size of star in pixels
-static inline int starSize(float noise, float peak, float seeing)
-{
-#ifdef STAR_GAUSSIAN
-    // Gaussian star (solving Gaussian)
-    float target = STAR_BG_FRAC * seeing * sqrt(M_2_PI) * noise / peak;
-    return target < 1.0 ? 2.0 * seeing * sqrtf(-logf(target)) + 0.5 : seeing * STAR_RANGE_MIN;
-#else
-    // Waussian star (solving Waussian where z >> 1 and peak/sqrt(bg) >> 1)
-    float target = STAR_BG_FRAC * noise / peak;
-    return PS_MAX(sqrtf(1.0 / target) + 0.5, seeing * STAR_RANGE_MIN);
-#endif
-}
-
-// Add a star into the signal and variance images
-bool ppSimInsertStar(psImage *signal,       // Signal image, to which to add star
-		     psImage *variance,     // Variance image, to which to add star
-		     float x0, float y0,    // Position of star
-		     float peak,            // Peak flux of star
-		     float noise,           // Rough noise estimate
-		     float seeing,          // Seeing for star
-		     const psImage *correction // Exposure correction as a function of position
-    )
-{
-    int size = starSize(noise, peak, seeing); // Size of star
-
-    // Range in image pixels on which to add star
-    int xMin = PS_MAX(0, x0 - size);
-    int xMax = PS_MIN(signal->numCols - 1, x0 + size);
-    int yMin = PS_MAX(0, y0 - size);
-    int yMax = PS_MIN(signal->numRows - 1, y0 + size);
-
-    for (int y = yMin; y <= yMax; y++) {
-        for (int x = xMin; x <= xMax; x++) {
-            float star = peak * correction->data.F32[y][x] * starFlux(x, y, x0, y0, seeing); // Star flux
-            signal->data.F32[y][x] += star;
-            variance->data.F32[y][x] += star;
-        }
-    }
-    return true;
-}
 
 void ppSimStarFree(ppSimStar *star)
@@ -74,2 +13,15 @@
     return star;
 }
+
+void ppSimGalaxyFree(ppSimGalaxy *galaxy)
+{
+    return;
+}
+
+ppSimGalaxy *ppSimGalaxyAlloc () {
+
+    ppSimGalaxy *galaxy = (ppSimGalaxy *) psAlloc(sizeof(ppSimGalaxy));
+    psMemSetDeallocator(galaxy, (psFreeFunc) ppSimGalaxyFree);
+
+    return galaxy;
+}