Index: trunk/ppSim/src/ppSimArguments.c
===================================================================
--- trunk/ppSim/src/ppSimArguments.c	(revision 14463)
+++ trunk/ppSim/src/ppSimArguments.c	(revision 14531)
@@ -8,5 +8,5 @@
 {
     fprintf(stderr, "\nPan-STARRS data simulator\n\n");
-    fprintf(stderr, "Usage: %s -camera CAMERA_NAME\n", program);
+    fprintf(stderr, "Usage: %s -camera CAMERA_NAME (output)\n", program);
     fprintf(stderr, "\n");
     psArgumentHelp(arguments);
@@ -68,12 +68,14 @@
     psMetadataAddS32(arguments, PS_LIST_TAIL, "-bin", 0, "Binning in x and y", 1);
 
+    pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf", "-psflist");
+
     if (!config->camera) {
         psErrorStackPrint(stderr, "A camera name must be specified using the -camera option.");
         usage(argv[0], arguments, config);
     }
-
-    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) {
-        usage(argv[0], arguments, config);
-    }
+ 
+    if (!psArgumentParse(arguments, &argc, argv)) { usage (argv[0], arguments, config); }
+
+    if (argc != 2) { usage(argv[0], arguments, config); }
 
     psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
@@ -113,5 +115,5 @@
         type = PPSIM_TYPE_OBJECT;
     } else {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised exposure type: %s", typeStr);
+        psErrorStackPrint(stderr, "Unrecognised exposure type: %s", typeStr);
         usage(argv[0], arguments, config);
     }
@@ -179,6 +181,5 @@
 
         if (isnan(ra0) || isnan(dec0) || isnan(pa) || isnan(scale) || isnan(zp) || isnan(seeing)) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, false,
-                    "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type");
+            psErrorStackPrint(stderr, "-ra, -dec, -pa, -scale, -zp, -seeing must be specified for OBJECT type");
             usage(argv[0], arguments, config);
         }
Index: trunk/ppSim/src/ppSimCreate.c
===================================================================
--- trunk/ppSim/src/ppSimCreate.c	(revision 14463)
+++ trunk/ppSim/src/ppSimCreate.c	(revision 14531)
@@ -24,4 +24,14 @@
     }
     file->save = true;
+
+    // have we supplied a psf model?
+    if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) {
+	bool status = false;
+        pmFPAfileBindFromArgs(&status, file, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF");
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to find/build PSPHOT.PSF.LOAD");
+            return status;
+        }
+    }
 
     pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes
Index: trunk/ppSim/src/ppSimInsertSource.c
===================================================================
--- trunk/ppSim/src/ppSimInsertSource.c	(revision 14531)
+++ trunk/ppSim/src/ppSimInsertSource.c	(revision 14531)
@@ -0,0 +1,55 @@
+# include "ppSim.h"
+
+bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
+
+    bool mdok;
+
+    assert (stars);
+    
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+
+    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
+    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;
+	}
+    }
+
+    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
+
+    // Rough noise estimate, appropriate for entire cell
+    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
+    pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSF");
+
+    for (long i = 0; i < stars->n; i++) {
+	ppSimStar *star = stars->data[i];
+
+	// Position on the cell and peak flux
+	float x = PPSIM_FPA_TO_CELL(star->x, x0Cell, xParityCell, binning, x0Chip, xParityChip);
+	float y = PPSIM_FPA_TO_CELL(star->y, y0Cell, yParityCell, binning, y0Chip, yParityChip);
+
+	pmModel *model = pmModelFromPSFforXY (psf, x, y, star->flux);
+
+	ppSimInsertSource(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr);
+    }
+    return true;
+}
+
Index: trunk/ppSim/src/ppSimInsertSources.c
===================================================================
--- trunk/ppSim/src/ppSimInsertSources.c	(revision 14531)
+++ trunk/ppSim/src/ppSimInsertSources.c	(revision 14531)
@@ -0,0 +1,60 @@
+# include "ppSim.h"
+
+bool ppSimInsertSources (psImage *signal, psImage *variance, psImage *expCorr, psArray *stars, pmConfig *config, pmChip *chip, pmCell *cell) {
+
+    bool mdok;
+
+    assert (stars);
+    
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
+
+    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
+    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;
+	}
+    }
+
+    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
+
+    // Rough noise estimate, appropriate for entire cell
+    float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
+    pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSF");
+
+    for (long i = 0; i < stars->n; i++) {
+	ppSimStar *star = stars->data[i];
+
+	// Position on the cell and peak flux
+	float x = PPSIM_FPA_TO_CELL(star->x, x0Cell, xParityCell, binning, x0Chip, xParityChip);
+	float y = PPSIM_FPA_TO_CELL(star->y, y0Cell, yParityCell, binning, y0Chip, yParityChip);
+
+	// instantiate a model for the PSF at this location (normalized or not?)
+	// XXX 
+	pmModel *model = pmModelFromPSFforXY (psf, star->x, star->y, star->flux);
+
+	// construct a source, with model flux pixels set, based on the model
+	pmSource *source = pmSourceFromModel (model);
+
+	ppSimInsertSource(signal, variance, x, y, star->flux, roughNoise, seeing, expCorr);
+    }
+    return true;
+}
+
Index: trunk/ppSim/src/ppSimLoadStars.c
===================================================================
--- trunk/ppSim/src/ppSimLoadStars.c	(revision 14463)
+++ trunk/ppSim/src/ppSimLoadStars.c	(revision 14531)
@@ -22,4 +22,7 @@
     float radius = 0.5 * PS_MAX(bounds->x1 - bounds->x0, bounds->y1 - bounds->y0) * scale;
     psFree(bounds);
+
+    float x0fpa = 0.5*(bounds->x0 + bounds->x1);
+    float y0fpa = 0.5*(bounds->y0 + bounds->y1);
 
     psMetadataAdd(astroRecipe, PS_LIST_TAIL, "RA_MIN",  PS_DATA_F32 | PS_META_REPLACE, "", ra0 - radius);
@@ -47,11 +50,11 @@
 
 	// Convert to x,y position on tangent plane, in pixels
-	float div = sin(star->ra) * sin(dec0) + cos(star->dec) * cos(dec0) * cos(star->ra - ra0) / scale;
+	float div = (sin(star->ra) * sin(dec0) + cos(star->dec) * cos(dec0) * cos(star->ra - ra0)) * scale;
 	float xi = cos(star->dec) * sin(star->ra - ra0) / div;
 	float eta = (sin(star->dec) * cos(dec0) - cos(star->dec) * sin(dec0) * cos(star->ra - ra0)) / div;
 
-	// Apply rotation
-	star->x = cos(pa) * xi - sin(pa) * eta;
-	star->y = sin(pa) * xi + cos(pa) * eta;
+	// Apply rotation, make FPA center of boresite
+	star->x = cos(pa) * xi - sin(pa) * eta + x0fpa;
+	star->y = sin(pa) * xi + cos(pa) * eta + y0fpa;
 
 	// Convert magnitude to peak flux
Index: trunk/ppSim/src/ppSimLoop.c
===================================================================
--- trunk/ppSim/src/ppSimLoop.c	(revision 14463)
+++ trunk/ppSim/src/ppSimLoop.c	(revision 14531)
@@ -23,5 +23,4 @@
 
     // Add random stars
-    // XXX put this in a wrapper (add to array of stars)
     if (type == PPSIM_TYPE_OBJECT) {
 	ppSimMakeStars (stars, fpa, config, rng);
@@ -30,8 +29,27 @@
     pmFPAview *view = pmFPAviewAlloc(0);// View for iterating over FPA
 
+    // load any needed files (eg, input image, PSF)
+    if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+	psError(PSPHOT_ERR_DATA, false, "failed IO for fpa in ppSim\n");
+	psFree(view);
+	return false;
+    }
+
     ppSimUpdateConceptsFPA (fpa, config);
 
     pmChip *chip;                       // Chip from FPA
     while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
+
+	// load any needed files (eg, input image, PSF)
+	if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) {
+            psError(PSPHOT_ERR_DATA, false, "failed IO for chip %d in ppSim\n", view->chip);
+	    psFree (view);
+	    return false;
+	}
+
+	pmPSF *psf = NULL;
+	if (type == PPSIM_TYPE_OBJECT) {
+	    psf = ppSimSetPSF (config);
+	}
 
         pmCell *cell;                   // Cell from chip
Index: trunk/ppSim/src/ppSimMakeStars.c
===================================================================
--- trunk/ppSim/src/ppSimMakeStars.c	(revision 14463)
+++ trunk/ppSim/src/ppSimMakeStars.c	(revision 14531)
@@ -63,9 +63,12 @@
 	    ppSimStar *star = ppSimStarAlloc ();
 
-	    star->x    = (psRandomUniform(rng) - 0.5) * xSize ; // x position
-	    star->y    = (psRandomUniform(rng) - 0.5) * ySize; // y position
+	    // 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;
     }
     return true;
Index: trunk/ppSim/src/ppSimUtils.c
===================================================================
--- trunk/ppSim/src/ppSimUtils.c	(revision 14463)
+++ trunk/ppSim/src/ppSimUtils.c	(revision 14531)
@@ -2,4 +2,5 @@
 
 // Generate a header containing WCS keywords
+// this function is called with only one of fpa, chip, cell not NULL
 bool ppSimInitHeader(pmConfig *config,
 		     pmFPA *fpa,
@@ -11,16 +12,6 @@
     float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
     float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
-    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE") *
-        M_PI / 3600.0 / 180.0; // Plate scale (radians/pixel)
-
-    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");
+    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // plate scale in arcsec / pixel
+    scale *= M_PI / 3600.0 / 180.0; // convert plate scale to radians/pixel
 
     int binning = psMetadataLookupS32(NULL, config->arguments, "BINNING"); // Binning in x and y
@@ -29,4 +20,14 @@
     int xParity, yParity;
     if (cell) {
+	int x0Chip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.X0");
+	int y0Chip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.Y0");
+	int xParityChip = psMetadataLookupS32(NULL, cell->parent->concepts, "CHIP.XPARITY");
+	int yParityChip = psMetadataLookupS32(NULL, cell->parent->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");
+
 	x0 = PPSIM_FPA_TO_CELL(0.0, x0Cell, xParityCell, binning, x0Chip, xParityChip);
 	y0 = PPSIM_FPA_TO_CELL(0.0, y0Cell, yParityCell, binning, y0Chip, yParityChip);
@@ -35,4 +36,9 @@
     }
     if (chip) {
+	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");
+
 	x0 = PPSIM_FPA_TO_CELL(0.0, 0, 1, binning, x0Chip, xParityChip);
 	y0 = PPSIM_FPA_TO_CELL(0.0, 0, 1, binning, y0Chip, yParityChip);