Changeset 12873

Timestamp:

Apr 17, 2007, 6:29:18 PM (19 years ago)

Author:

Paul Price

Message:

Fixing bias, adding catalog and random stars.

Location:

trunk/ppSim/src

Files:

• ppSimArguments.c (modified) (6 diffs)
• ppSimCreate.c (modified) (5 diffs)
• ppSimLoop.c (modified) (10 diffs)

trunk/ppSim/src/ppSimArguments.c

@@ -16 +16 @@
 {
     fprintf(stderr, "\nPan-STARRS data simulator\n\n");
-    fprintf(stderr, "Usage: %s -camera CAMERA_NAME -format FORMAT_NAME\n", config->argv[0]);
+    fprintf(stderr, "Usage: %s -camera CAMERA_NAME\n", config->argv[0]);
     fprintf(stderr, "\n");
     psArgumentHelp(arguments);
@@ -35 +35 @@
     float value = psMetadataLookupF32(NULL, arguments, argName); // Bias level
     if (isnan(value)) {
-        value = psMetadataLookupF32(NULL, recipe, recipeName);
-        if (isnan(value)) {
+        bool mdok;                      // Status of MD lookup
+        value = psMetadataLookupF32(&mdok, recipe, recipeName);
+        if (!mdok) {
             psErrorStackPrint(stderr, "Unable to find %s in recipe %s", recipeName, PPSIM_RECIPE);
             psFree((psPtr)arguments);
@@ -43 +44 @@
         }
     }
-
     return psMetadataAddF32(target, PS_LIST_TAIL, recipeName, 0, NULL, value);
 }
-
 
 void ppSimArguments(pmConfig *config)
@@ -64 +63 @@
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-shuttertime", 0, "Shutter time (s)", NAN);
     psMetadataAddF32(arguments, PS_LIST_TAIL, "-skyrate", 0, "Sky rate (e/s)", NAN);
+    psMetadataAddStr(arguments, PS_LIST_TAIL, "-stars", 0, "File with star RA, Dec, Mag", NULL);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-ra", 0, "RA (degrees)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-dec", 0, "Dec (degrees)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-pa", 0, "Position angle (degrees)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-scale", 0, "Plate scale (arcsec/pixel)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-zp", 0, "Photometric zero point", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-seeing", 0, "Seeing FWHM (arcsec)", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-starslum", 0, "Fake star luminosity function slope", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsmag", 0, "Brightest magnitude for fake stars", NAN);
+    psMetadataAddF32(arguments, PS_LIST_TAIL, "-starsdensity", 0, "Density of fake stars at magnitude", NAN);
 
     if (!config->camera) {
@@ -74 +83 @@
     }
 
-    const char *formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
-    if (!formatName) {
-        psErrorStackPrint(stderr, "A camera format must be specified using -format");
-        usage(arguments, config);
-    }
-    psMetadataAddStr(config->arguments, PS_LIST_TAIL, "FORMAT", 0, "Camera format name", formatName);
+    psString formatName = psMetadataLookupStr(NULL, arguments, "-format"); // Name of format
+    if (formatName) {
+        config->formatName = psMemIncrRefCounter(formatName);
+
+        psMetadata *formats = psMetadataLookupMetadata(NULL, config->camera, "FORMATS"); // The camera formats
+        if (!formats) {
+            psErrorStackPrint(stderr, "Unable to find FORMATS in camera configuration.");
+            psFree(arguments);
+            psFree(config);
+            exit(PS_EXIT_CONFIG_ERROR);
+        }
+        psMetadata *format = psMetadataLookupMetadata(NULL, formats, formatName); // Format of interest
+        if (!format) {
+            psErrorStackPrint(stderr, "Unable to find format %s in camera FORMATS.", formatName);
+            psFree(arguments);
+            psFree(config);
+            exit(PS_EXIT_CONFIG_ERROR);
+        }
+        config->format = psMemIncrRefCounter(format);
+    }
 
     const char *typeStr = psMetadataLookupStr(NULL, arguments, "-type"); // Exposure type
@@ -124 +147 @@
     }
 
-    valueArgRecipe(config, arguments, "-biaslevel",   recipe, "BIAS.LEVEL",   config->arguments);
-    valueArgRecipe(config, arguments, "-biasrange",   recipe, "BIAS.RANGE",   config->arguments);
-    valueArgRecipe(config, arguments, "-darkrate",    recipe, "DARK.RATE",    config->arguments);
-    valueArgRecipe(config, arguments, "-flatsigma",   recipe, "FLAT.SIGMA",   config->arguments);
-    valueArgRecipe(config, arguments, "-flatrate",    recipe, "FLAT.RATE",    config->arguments);
-    valueArgRecipe(config, arguments, "-shuttertime", recipe, "SHUTTER.TIME", config->arguments);
-    valueArgRecipe(config, arguments, "-skyrate",     recipe, "SKY.RATE",     config->arguments);
+    valueArgRecipe(config, arguments, "-biaslevel",    recipe, "BIAS.LEVEL",    config->arguments);
+    valueArgRecipe(config, arguments, "-biasrange",    recipe, "BIAS.RANGE",    config->arguments);
+    valueArgRecipe(config, arguments, "-darkrate",     recipe, "DARK.RATE",     config->arguments);
+    valueArgRecipe(config, arguments, "-flatsigma",    recipe, "FLAT.SIGMA",    config->arguments);
+    valueArgRecipe(config, arguments, "-flatrate",     recipe, "FLAT.RATE",     config->arguments);
+    valueArgRecipe(config, arguments, "-shuttertime",  recipe, "SHUTTER.TIME",  config->arguments);
+    valueArgRecipe(config, arguments, "-skyrate",      recipe, "SKY.RATE",      config->arguments);
+    valueArgRecipe(config, arguments, "-starslum",     recipe, "STARS.LUM",     config->arguments);
+    valueArgRecipe(config, arguments, "-starsmag",     recipe, "STARS.MAG",     config->arguments);
+    valueArgRecipe(config, arguments, "-starsdensity", recipe, "STARS.DENSITY", config->arguments);
+
+    bool mdok;                          // Status of MD lookup
+    int biasOrder = psMetadataLookupS32(&mdok, recipe, "BIAS.ORDER"); // Overscan polynomial order
+    if (!mdok) {
+        psWarning("BIAS.ORDER(S32) is not set in the recipe %s --- assuming %d", PPSIM_RECIPE, biasOrder);
+    }
+    psMetadataAddS32(config->arguments, PS_LIST_TAIL, "BIAS.ORDER", 0,
+                     "Overscan polynomial order", biasOrder);
+
 
     psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image",
                      config->argv[1]);
 
+
+    // Values common for stars
+    float scale = psMetadataLookupF32(NULL, arguments, "-scale"); // Plate scale
+    float zp = psMetadataLookupF32(NULL, arguments, "-zp"); // Zero point
+    float seeing = psMetadataLookupF32(NULL, arguments, "-seeing"); // Zero point
+    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SCALE", 0, "Plate scale (arcsec/pix)", scale);
+    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "ZEROPOINT", 0, "Photometric zeropoint", zp);
+    psMetadataAddF32(config->arguments, PS_LIST_TAIL, "SEEING", 0, "Seeing simga (pix)",
+                     seeing / 2.0 / sqrt(2.0 * log(2.0)) / scale);
+
+    const char *starName = psMetadataLookupStr(NULL, arguments, "-stars"); // Name of file with stars
+    if (starName) {
+        float ra0 = psMetadataLookupF32(NULL, arguments, "-ra"); // Right Ascension of boresight
+        float dec0 = psMetadataLookupF32(NULL, arguments, "-dec"); // Declination of boresight
+        float pa = psMetadataLookupF32(NULL, arguments, "-pa"); // Position angle
+
+        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 with -stars");
+            usage(arguments, config);
+        }
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "RA", 0, "Boresight RA (radians)",
+                         ra0 * M_PI / 180.0);
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "DEC", 0, "Boresight Declination (radians)",
+                         dec0 * M_PI / 180.0);
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "PA", 0, "Boresight position angle (radians)",
+                         pa * M_PI / 180.0);
+
+        psArray *raDecMag = psVectorsReadFromFile(starName, "%f %f %f"); // Vectors for RA Dec Mag
+        if (!raDecMag || raDecMag->n != 3) {
+            psError(PS_ERR_IO, false, "Unable to read RA Dec Mag from %s", starName);
+            psFree(raDecMag);
+            psFree(arguments);
+            psFree(config);
+            exit(PS_EXIT_CONFIG_ERROR);
+        }
+
+        psVector *ra = raDecMag->data[0]; // Star RAs
+        psVector *dec = raDecMag->data[1]; // Star Declinations
+        psVector *mag = raDecMag->data[2]; // Star magnitudes
+        for (long i = 0; i < ra->n; i++) {
+            ra->data.F32[i] *= M_PI / 180.0;
+            dec->data.F32[i] *= M_PI / 180.0;
+        }
+
+        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.RA", 0, "Star RAs (radians)", ra);
+        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.DEC", 0, "Star Decs (radians)", dec);
+        psMetadataAddVector(config->arguments, PS_LIST_TAIL, "CATALOG.MAG", 0, "Star magnitudes", mag);
+        psFree(raDecMag);
+    }
+
+    if (psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY") > 0 &&
+        (isnan(scale) || isnan(zp) || isnan(seeing))) {
+        psWarning("-scale, -zp, -seeing must be specified if fake stars are to be generated --- turned off.");
+        psMetadataAddS32(config->arguments, PS_LIST_TAIL, "STASR.NUM", PS_META_REPLACE,
+                         "Fake stars turned off", 0);
+    }
+
+
     psFree(arguments);
+    return;
 }

trunk/ppSim/src/ppSimCreate.c

@@ -13 +13 @@
     PS_ASSERT_PTR_NON_NULL(config, NULL);
 
-    psString formatName = psMetadataLookupStr(NULL, config->arguments, "FORMAT"); // Input format name
-    config->formatName = psMemIncrRefCounter(formatName);
-
-    psMetadata *formats = psMetadataLookupMetadata(NULL, config->camera, "FORMATS"); // The camera formats
-    if (!formats) {
-        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find FORMATS in camera configuration.");
-        return NULL;
-    }
-    psMetadata *format = psMetadataLookupMetadata(NULL, formats, formatName); // Format of interest
-    if (!format) {
-        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find format %s in camera FORMATS.", formatName);
-        return NULL;
-    }
-    pmFPALevel phuLevel = pmFPAPHULevel(format); // Level at which PHU goes
-
     pmFPA *fpa = pmFPAConstruct(config->camera); // FPA to contain the observation
     if (!fpa) {
@@ -34 +19 @@
     }
 
+    pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE);
+    psFree(fpa);
+    if (!file) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s.  "
+                "Did you forget to specify the format?", OUTPUT_FILE);
+        return NULL;
+    }
+    if (file->type != PM_FPA_FILE_IMAGE) {
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE);
+        psFree(file);
+        return NULL;
+    }
+
+    pmFPALevel phuLevel = pmFPAPHULevel(file->format); // Level at which PHU goes
+
     pmFPAview *view = pmFPAviewAlloc(0);// View for current level
 
     if (phuLevel == PM_FPA_LEVEL_FPA) {
-        if (!pmFPAAddSourceFromView(fpa, "Simulation", view, format)) {
+        if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
             psFree(fpa);
@@ -48 +48 @@
     while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
         if (phuLevel == PM_FPA_LEVEL_CHIP) {
-            if (!pmFPAAddSourceFromView(fpa, "Simulation", view, format)) {
+            if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
                 psFree(fpa);
@@ -59 +59 @@
         while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
             if (phuLevel == PM_FPA_LEVEL_CELL) {
-                if (!pmFPAAddSourceFromView(fpa, "Simulation", view, format)) {
+                if (!pmFPAAddSourceFromView(fpa, "Simulation", view, file->format)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to FPA.");
                     psFree(fpa);
@@ -70 +70 @@
     psFree(view);
 
-    pmFPAfile *file = pmFPAfileDefineOutput(config, fpa, OUTPUT_FILE);
-    psFree(fpa);
-    if (!file) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to create output file from %s", OUTPUT_FILE);
-        return NULL;
-    }
-    if (file->type != PM_FPA_FILE_IMAGE) {
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true, "%s type is not IMAGE", OUTPUT_FILE);
-        psFree(file);
-        return NULL;
-    }
-
     return file;
 }

trunk/ppSim/src/ppSimLoop.c

@@ -9 +9 @@
 #include "ppSim.h"
 
-#define COMPARE(VALUE,MIN,MAX) \
-    if (VALUE < MIN) { \
-        MIN = VALUE; \
-    } \
-    if (VALUE > MAX) { \
-        MAX = VALUE; \
-    }
+#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
+
+#define LOG10(X) (log(X)/log(10))       // Base 10 logarithm
+
+// Return FPA position, given a cell position; calculations are all done in pixel units
+static inline float cell2fpa(float pos, // Cell position
+                             int cell0, // Cell offset
+                             int cellParity, // Cell parity
+                             int chip0, // Chip offset
+                             int chipParity // Chip parity
+    )
+{
+    return (chip0 + chipParity * (cell0 + cellParity * pos));
+}
+
+// Return cell position, given an FPA position; calculations are all done in pixel units
+static inline float fpa2cell(float pos, // FPA position
+                             int cell0, // Cell offset
+                             int cellParity, // Cell parity
+                             int chip0, // Chip offset
+                             int chipParity // Chip parity
+    )
+{
+    return ((pos - chip0) * chipParity - cell0) * cellParity;
+}
+
+// Compare a value with minimum and maximum values, replacing where required.
+#define COMPARE(VALUE,MIN,MAX) { \
+        if (VALUE < MIN) { \
+            MIN = VALUE; \
+        } \
+        if (VALUE > MAX) { \
+            MAX = VALUE; \
+        } \
+    }
+
+// Return bounds of a chip, based on the concepts
+static psRegion *chipBounds(const pmChip *chip, // Chip for which to determine size
+                            pmFPAview *view // View for chip
+    )
+{
+    assert(chip);
+    assert(view);
+
+    // Bounds of chip
+    int xMin = INT_MAX;
+    int xMax = - INT_MAX;
+    int yMin = INT_MAX;
+    int yMax = - INT_MAX;
+
+    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");
+
+    if ((xParityChip != 1 && xParityChip != -1) || (yParityChip != 1 && yParityChip != -1)) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Chip parities are not set.");
+        psFree(view);
+        return NULL;
+    }
+
+    pmCell *cell;                   // Cell from chip
+    while ((cell = pmFPAviewNextCell(view, chip->parent, 1))) {
+        int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
+        int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");
+
+        if (xSize == 0 || ySize == 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell sizes are not set.");
+            psFree(view);
+            return NULL;
+        }
+
+        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");
+
+        if ((xParityCell != 1 && xParityCell != -1) || (yParityCell != 1 && yParityCell != -1)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell parities are not set.");
+            psFree(view);
+            return NULL;
+        }
+
+        // "Left", "Right", "Bottom" and "Top" don't take into account the parity
+        int cellLeft = cell2fpa(0, x0Cell, xParityCell, x0Chip, xParityChip);
+        int cellRight = cell2fpa(xSize, x0Cell, xParityCell, x0Chip, xParityChip);
+        int cellBottom = cell2fpa(0, y0Cell, yParityCell, y0Chip, yParityChip);
+        int cellTop = cell2fpa(ySize, y0Cell, yParityCell, y0Chip, yParityChip);
+
+        COMPARE(cellLeft, xMin, xMax);
+        COMPARE(cellRight, xMin, xMax);
+        COMPARE(cellBottom, yMin, yMax);
+        COMPARE(cellTop, yMin, yMax);
+    }
+
+    return psRegionAlloc(xMin, xMax, yMin, yMax);
+}
 
 // Return bounds of an FPA, based on the concepts
@@ -33 +124 @@
     pmChip *chip;                       // Chip from FPA
     while ((chip = pmFPAviewNextChip(view, fpa, 1))) {
-        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");
-
-        if ((xParityChip != 1 && xParityChip != -1) || (yParityChip != 1 && yParityChip != -1)) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Chip parities are not set.");
-            psFree(view);
-            return NULL;
-        }
-
-        pmCell *cell;                   // Cell from chip
-        while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
-            int xSize = psMetadataLookupS32(NULL, cell->concepts, "CELL.XSIZE");
-            int ySize = psMetadataLookupS32(NULL, cell->concepts, "CELL.YSIZE");
-
-            if (xSize == 0 || ySize == 0) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell sizes are not set.");
-                psFree(view);
-                return NULL;
-            }
-
-            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");
-
-            if ((xParityCell != 1 && xParityCell != -1) || (yParityCell != 1 && yParityCell != -1)) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cell parities are not set.");
-                psFree(view);
-                return NULL;
-            }
-
-            // "Left", "Right", "Bottom" and "Top" don't take into account the parity
-            int cellLeft = x0Chip + xParityChip * x0Cell;
-            int cellRight = x0Chip + xParityChip * (x0Cell + xParityCell * xSize);
-            int cellBottom = y0Chip + yParityChip * y0Cell;
-            int cellTop = y0Chip + yParityChip * (y0Cell + yParityCell * ySize);
-
-            COMPARE(cellLeft, xMin, xMax);
-            COMPARE(cellRight, xMin, xMax);
-            COMPARE(cellBottom, yMin, yMax);
-            COMPARE(cellTop, yMin, yMax);
-        }
+        psRegion *bounds = chipBounds(chip, view); // Bounds for chip
+        COMPARE(bounds->x0, xMin, xMax);
+        COMPARE(bounds->x1, xMin, xMax);
+        COMPARE(bounds->y0, yMin, yMax);
+        COMPARE(bounds->y1, yMin, yMax);
+        psFree(bounds);
     }
 
     return psRegionAlloc(xMin, xMax, yMin, yMax);
 }
+
+
+// Add noise to an image
+static psImage *addNoise(psImage *signal, // Signal image, modified and returned
+                         const psImage *variance, // Variance image
+                         psRandom *rng, // Random number generator
+                         float gain     // CCD gain (e/ADU)
+    )
+{
+    assert(signal->type.type == PS_TYPE_F32);
+    assert(signal->numCols == variance->numCols && signal->numRows == variance->numRows);
+
+    // Add the noise into the image
+    for (int y = 0; y < signal->numRows; y++) {
+        for (int x = 0; x < signal->numCols; x++) {
+            signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x]) * psRandomGaussian(rng);
+            signal->data.F32[y][x] /= gain; // Converting to ADU
+        }
+    }
+
+    return signal;
+}
+
+// Apply saturation limit to image
+static void saturate(psImage *image, // Image to apply saturation
+                     float saturation // Saturation level
+    )
+{
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (image->data.F32[y][x] > saturation) {
+                image->data.F32[y][x] = saturation;
+            }
+        }
+    }
+    return;
+}
+
+// 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
+static void star(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;
+}
+
 
 
@@ -92 +241 @@
     assert(fpa);
 
+    bool mdok;                          // Status of MD lookups
+
     psRegion *bounds = fpaBounds(fpa);  // Bounds of FPA
     if (!bounds || (bounds->x0 == 0.0 && bounds->x1 == 0.0) || (bounds->y0 == 0.0 && bounds->y1 == 0.0)) {
@@ -98 +249 @@
     }
 
+    psVector *ra = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.RA"); // Star RAs
+    psVector *dec = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.DEC"); // Star Declinations
+    psVector *mag = psMetadataLookupPtr(NULL, config->arguments, "CATALOG.MAG"); // Star magnitudes
+    float seeing = psMetadataLookupF32(NULL, config->arguments, "SEEING"); // Seeing sigma (pix)
+    float scale = psMetadataLookupF32(NULL, config->arguments, "SCALE"); // Plate scale (arcsec/pixel)
+    scale *= M_PI / 3600.0 / 180.0; // Convert to radians/pixel
+    float zp = psMetadataLookupF32(NULL, config->arguments, "ZEROPOINT"); // Photometric zero point
+
+    if (ra && dec && mag) {
+        // Convert star ra,dec to pixels from FPA centre
+        float ra0 = psMetadataLookupF32(NULL, config->arguments, "RA"); // Boresight RA (radians)
+        float dec0 = psMetadataLookupF32(NULL, config->arguments, "DEC"); // Boresight Dec (radians)
+        float pa = psMetadataLookupF32(NULL, config->arguments, "PA"); // Position angle (radians)
+
+        // Conversion loop
+        for (long i = 0; i < ra->n; i++) {
+            float div = sin(dec->data.F32[i]) * sin(dec0) +
+                cos(dec->data.F32[i]) * cos(dec0) * cos(ra->data.F32[i] - ra0); // Divisor
+
+            // Convert to x,y position on tangent plane, in pixels
+            float xi = cos(dec->data.F32[i]) * sin(ra->data.F32[i] - ra0) / div / scale;
+            float eta = (sin(dec->data.F32[i]) * cos(dec0) -
+                cos(dec->data.F32[i]) * sin(dec0) * cos(ra->data.F32[i] - ra0)) / div / scale;
+
+            // Apply rotation
+            ra->data.F32[i] = cos(pa) * xi - sin(pa) * eta;
+            dec->data.F32[i] = sin(pa) * xi + cos(pa) * eta;
+            ra->data.F32[i] = xi;
+            dec->data.F32[i] = eta;
+
+            // Convert magnitude to peak flux
+            mag->data.F32[i] = powf(10.0, -0.4 * (mag->data.F32[i] - zp)) / sqrt(2.0*M_PI) / seeing;
+        }
+    }
+
     float biasLevel = psMetadataLookupF32(NULL, config->arguments, "BIAS.LEVEL"); // Bias level
     float biasRange = psMetadataLookupF32(NULL, config->arguments, "BIAS.RANGE"); // Bias range
+    int biasOrder = psMetadataLookupS32(NULL, config->arguments, "BIAS.ORDER"); // Bias order
     float darkRate = psMetadataLookupF32(NULL, config->arguments, "DARK.RATE"); // Dark rate
     float flatSigma = psMetadataLookupF32(NULL, config->arguments, "FLAT.SIGMA"); // Flat-field coefficient
@@ -105 +292 @@
     float shutterTime = psMetadataLookupF32(NULL, config->arguments, "SHUTTER.TIME"); // Shutter time
     float skyRate = psMetadataLookupF32(NULL, config->arguments, "SKY.RATE"); // Sky rate
+    float starsLum = psMetadataLookupF32(NULL, config->arguments, "STARS.LUM"); // Star luminosity func slope
+    float starsMag = psMetadataLookupF32(NULL, config->arguments, "STARS.MAG"); // Star brightest magnitude
+    float starsDensity = psMetadataLookupF32(NULL, config->arguments, "STARS.DENSITY"); // Density of fakes
 
     float expTime = psMetadataLookupF32(NULL, config->arguments, "EXPTIME"); // Exposure time
     const char *filter = psMetadataLookupStr(NULL, config->arguments, "FILTER"); // Filter name
     ppSimType type = psMetadataLookupS32(NULL, config->arguments, "TYPE"); // Type of image to simulate
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSIM_RECIPE); // Recipe
 
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
@@ -139 +331 @@
         int yParityChip = psMetadataLookupS32(NULL, chip->concepts, "CHIP.YPARITY");
 
+        // Correct chip offsets so that boresight is in the middle of the FPA
+        x0Chip -= 0.5 * (bounds->x1 - bounds->x0);
+        y0Chip -= 0.5 * (bounds->y1 - bounds->y0);
+
+#if 0
+        // Add random stars
+        if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) {
+            // Peak fluxes: faintest and brightest levels for random stars
+            float faint = roughNoise;
+            float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
+            if (bright < faint) {
+                psLogMsg("ppSim", PS_LOG_INFO,
+                         "Image noise is above brightest random star --- no random stars added.");
+            } else {
+                // Normalisation, set by the specified stellar density at the specified bright
+                // magnitude
+                float norm = starsDensity * numCols * numRows * 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;
+
+                printf("Adding %ld stars to %f\n", num, faint);
+
+                for (long j = 0; j < num; j++) {
+                    float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star
+
+                    // TO DO: put stars on chips, rather than cells, so they can overlap cell
+                    // boundaries
+
+                    // Position on cell
+                    float x = psRandomUniform(rng) * numCols;
+                    float y = psRandomUniform(rng) * numRows;
+                    star(signal, variance, x, y, flux, roughNoise, seeing, expCorr);
+                }
+            }
+        }
+#endif
+
         pmCell *cell;                   // Cell from chip
         while ((cell = pmFPAviewNextCell(view, fpa, 1))) {
@@ -152 +383 @@
             float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // CCD gain, e/ADU
             if (isnan(gain)) {
-                psWarning("CELL.GAIN is not set; reverting to recipe value.");
-                gain = psMetadataLookupF32(NULL, config->arguments, "GAIN");
-            }
-            float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.RDNOISE"); // CCD read noise, e
+                psWarning("CELL.GAIN is not set; reverting to recipe value GAIN.");
+                gain = psMetadataLookupF32(&mdok, recipe, "GAIN");
+                if (!mdok) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find GAIN in recipe.");
+                    psFree(bounds);
+                    psFree(rng);
+                    return PS_EXIT_CONFIG_ERROR;
+                }
+            }
+            float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE");// CCD read noise, e
             if (isnan(readnoise)) {
-                psWarning("CELL.RDNOISE is not set; reverting to recipe value.");
-                readnoise = psMetadataLookupF32(NULL, config->arguments, "RDNOISE");
+                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.");
+                    psFree(bounds);
+                    psFree(rng);
+                    return PS_EXIT_CONFIG_ERROR;
+                }
+            }
+
+            float saturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
+            if (isnan(saturation)) {
+                psWarning("CELL.SATURATION is not set; reverting to recipe value SATURATION.");
+                readnoise = psMetadataLookupF32(&mdok, recipe, "SATURATION");
+                if (!mdok) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find SATURATION in recipe.");
+                    psFree(bounds);
+                    psFree(rng);
+                    return PS_EXIT_CONFIG_ERROR;
+                }
+            }
+
+            int readdir = psMetadataLookupS32(NULL, cell->concepts, "CELL.READDIR"); // Read direction
+            if (readdir != 1) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true, "CELL.READDIR = 1 is the only supported mode.");
+                psFree(bounds);
+                psFree(rng);
+                return PS_EXIT_CONFIG_ERROR;
+            }
+
+            psList *biassec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC"); // Bias regions
+            int numBias = psListLength(biassec); // Number of bias sections
+            psVector *biasCols;
+            if (numBias > 0) {
+                biasCols = psVectorAlloc(numBias, PS_TYPE_S32);
+                psListIterator *iter = psListIteratorAlloc(biassec, PS_LIST_HEAD, false); // Iterator
+                psRegion *bias;         // Bias region, from iteration
+                int i = 0;              // Counter
+                while ((bias = psListGetAndIncrement(iter))) {
+                    biasCols->data.S32[i++] = bias->x1 = bias->x0;
+                }
+            } else {
+                biasCols = psVectorAlloc(1, PS_TYPE_S32);
+                biasCols->data.S32[0] = psMetadataLookupS32(&mdok, recipe, "OVERSCAN.SIZE");
+                if (!mdok) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find OVERSCAN.SIZE in recipe.");
+                    psFree(biasCols);
+                    psFree(bounds);
+                    psFree(rng);
+                    return PS_EXIT_CONFIG_ERROR;
+                }
+                psRegion *newBias = psRegionAlloc(NAN, NAN, NAN, NAN); // New bias region, to be set later
+                biassec = psListAlloc(newBias);
+                psFree(newBias);
+                psMetadataAdd(cell->concepts, PS_LIST_TAIL, "CELL.BIASSEC", PS_DATA_LIST | PS_META_REPLACE,
+                              "Bias sections", biassec);
+                psFree(biassec);
+                numBias = 1;
             }
 
@@ -172 +465 @@
                 psImage *variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Noise in pixels
 
+                // Polynomial for bias
+                psPolynomial1D *biasPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, biasOrder);
+                for (int j = 0; j < biasOrder + 1; j++) {
+                    biasPoly->coeff[j] = biasRange * psRandomGaussian(rng);
+                }
+                psVector *rowBias = psVectorAlloc(numRows, PS_TYPE_F32); // Bias value, per row
+                int biasOffset = 0.5 * numRows * psRandomUniform(rng); // Offset to prevent common pattern
+
+                psImage *expCorr = NULL; // Exposure correction per pixel, for adding objects
+                if (type == PPSIM_TYPE_OBJECT) {
+                    expCorr = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+                }
+
                 for (int y = 0; y < numRows; y++) {
                     // Adjust bias level for this row
-                    float rowBias = biasLevel += psRandomGaussian(rng) * biasRange;
-                    float yFPA = (y0Chip + yParityChip * (y0Cell + yParityCell * y) - 2 * bounds->y0) /
+                    rowBias->data.F32[y] = psPolynomial1DEval(biasPoly, (float)(y + biasOffset) /
+                                                              (float)numRows - 0.5) + biasLevel;
+                    float yFPA = cell2fpa(y, y0Cell, yParityCell, y0Chip, yParityChip) * 2.0 /
                         (bounds->y1 - bounds->y0); // Relative y position in FPA
 
                     for (int x = 0; x < numCols; x++) {
-                        float xFPA = (x0Chip + xParityChip * (x0Cell + xParityCell * x) - 2 * bounds->x0) /
+                        float xFPA = cell2fpa(x, x0Cell, xParityCell, x0Chip, xParityChip) * 2.0 /
                             (bounds->x1 - bounds->x0); // Relative x position in FPA
 
                         // Bias level
-                        signal->data.F32[y][x] = rowBias;
+                        signal->data.F32[y][x] = rowBias->data.F32[y];
                         variance->data.F32[y][x] = PS_SQR(readnoise);
 
@@ -211 +518 @@
                             continue;
                         }
-
+                        expCorr->data.F32[y][x] = flatValue * realExpTime / expTime;
                         // Sky background
                         float skyFlux = skyRate * realExpTime * flatValue; // Flux from sky
@@ -219 +526 @@
                         // TO DO: Add fringes
 
-                        // TO DO: Add stars
                     }
                 }
 
-                // TO DO: Add overscan
-
-#if 0
-                // Add the noise into the image
-                for (int y = 0; y < numRows; y++) {
-                    for (int x = 0; x < numCols; x++) {
-                        signal->data.F32[y][x] /= gain; // Converting to ADU
-                        signal->data.F32[y][x] += sqrtf(variance->data.F32[y][x] / gain) *
-                            psRandomGaussian(rng);
+                // Rough noise estimate, appropriate for entire cell
+                float roughNoise = sqrtf(PS_SQR(readnoise) + (darkRate + skyRate) * expTime);
+
+
+                // XXX Not sure I've got the scaling right on the peak flux: it's multiplied by the expTime
+                // here, and the exposure up above.
+
+                // Add specified (catalogue) stars
+                if (type == PPSIM_TYPE_OBJECT && ra && dec && mag) {
+                    for (long j = 0; j < ra->n; j++) {
+                        // Position on the cell
+                        float x = fpa2cell(ra->data.F32[j], x0Cell, xParityCell, x0Chip, xParityChip);
+                        float y = fpa2cell(dec->data.F32[j], y0Cell, yParityCell, y0Chip, yParityChip);
+                        star(signal, variance, x, y, mag->data.F32[j] * expTime,
+                             roughNoise, seeing, expCorr);
                     }
                 }
-#endif
-                readout->image = signal;
+
+                // Add random stars
+                if (type == PPSIM_TYPE_OBJECT && starsDensity > 0) {
+                    // Peak fluxes: faintest and brightest levels for random stars
+                    float faint = roughNoise;
+                    float bright = powf(10.0, -0.4 * (starsMag - zp)) / sqrt(2.0*M_PI) / seeing * expTime;
+                    if (bright < faint) {
+                        psLogMsg("ppSim", PS_LOG_INFO,
+                                 "Image noise is above brightest random star --- no random stars added.");
+                    } else {
+                        // Normalisation, set by the specified stellar density at the specified bright
+                        // magnitude
+                        float norm = starsDensity * numCols * numRows * 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;
+
+                        printf("Adding %ld stars to %f\n", num, faint);
+
+                        for (long j = 0; j < num; j++) {
+                            float flux = expf((logf(j + 1) - logf(norm)) / starsLum); // Flux of star
+
+                            // TO DO: put stars on chips, rather than cells, so they can overlap cell
+                            // boundaries
+
+                            // Position on cell
+                            float x = psRandomUniform(rng) * numCols;
+                            float y = psRandomUniform(rng) * numRows;
+                            star(signal, variance, x, y, flux, roughNoise, seeing, expCorr);
+                        }
+                    }
+                }
+
+                readout->image = addNoise(signal, variance, rng, gain);
+                saturate(readout->image, saturation);
                 psFree(variance);
 
+                psFree(expCorr);
+
+                // Add overscan
+                for (int j = 0; j < numBias; j++) {
+                    psImage *signal = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Overscan
+                    for (int y = 0; y < signal->numRows; y++) {
+                        for (int x = 0; x < signal->numCols; x++) {
+                            signal->data.F32[y][x] = rowBias->data.F32[y];
+                        }
+                    }
+                    psImage *variance = psImageAlloc(biasCols->data.S32[j], numRows, PS_TYPE_F32); // Variance
+                    psImageInit(variance, PS_SQR(readnoise));
+                    addNoise(signal, variance, rng, gain);
+                    psFree(variance);
+
+                    psListAdd(readout->bias, PS_LIST_TAIL, signal);
+                    psFree(signal);     // Drop reference
+                }
+
+                psFree(rowBias);
                 psFree(readout);        // Drop reference
             }
+
 
             psMetadataAddF32(cell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE,