Index: trunk/psModules/src/camera/pmReadoutFake.c
===================================================================
--- trunk/psModules/src/camera/pmReadoutFake.c	(revision 15756)
+++ trunk/psModules/src/camera/pmReadoutFake.c	(revision 15813)
@@ -22,10 +22,39 @@
 #include "pmModelUtils.h"
 
-#define MODEL_TYPE "PS_MODEL_RGAUSS"     // Type of model to use
-
-
-bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
-                              const psVector *xOffset, const psVector *yOffset, pmPSF *psf, float minFlux,
-                              int radius)
+#include "pmReadoutFake.h"
+
+#define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
+#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
+
+
+#ifdef PAP_WORK
+// Given an object model, circularise it by setting the axes to be identical
+static bool circulariseModel(pmModel *model // Model to circularise
+    )
+{
+    assert(model);
+
+    psF32 *params = model->params->data.F32; // Model parameters
+    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
+    // Curiously, the minor axis can be larger than the major axis, so need to check.
+    if (axes.major >= axes.minor) {
+        axes.minor = axes.major;
+    } else {
+        axes.major = axes.minor;
+    }
+    return pmPSF_AxesToModel(params, axes);
+}
+#endif
+
+bool pmReadoutFakeFromSources(
+#ifdef PAP_WORK
+    pmReadout *readout, int numCols, int numRows, const psArray *sources,
+    const psVector *xOffset, const psVector *yOffset, pmPSF *psf,
+    float minFlux, int radius, bool circularise
+#else
+    pmReadout *readout, int numCols, int numRows, const psArray *sources,
+    float fwhm, float minFlux
+#endif
+    )
 {
     PS_ASSERT_PTR_NON_NULL(readout, false);
@@ -33,4 +62,6 @@
     PS_ASSERT_INT_LARGER_THAN(numRows, 0, false);
     PS_ASSERT_ARRAY_NON_NULL(sources, false);
+
+#ifdef PAP_WORK
     if (xOffset || yOffset) {
         PS_ASSERT_VECTOR_NON_NULL(xOffset, false);
@@ -51,5 +82,5 @@
         return false;
     }
-
+#endif
 
     readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
@@ -58,5 +89,5 @@
     int numSources = sources->n;          // Number of stars
 
-#if 0
+#ifndef PAP_WORK
     pmModelType modelType = pmModelClassGetType(MODEL_TYPE); // Type of PSF model
     assert(modelType >= 0);
@@ -85,9 +116,19 @@
         psAbort("Unsupported model type: %d", modelType);
     }
-#endif
-
+#else
     pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0,
                                              1.0); // Fake model, with central intensity of 1.0
+#endif
+
     float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0
+
+#ifdef PAP_WORK
+    if (circularise && !circulariseModel(fakeModel)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
+        psFree(fakeModel);
+        return false;
+    }
+    psFree(fakeModel);
+#endif
 
     for (int i = 0; i < numSources; i++) {
@@ -105,7 +146,20 @@
         }
 
+#ifdef PAP_WORK
+        pmModel *fakeModel = pmModelFromPSFforXY(psf, x, y, powf(10.0, -0.4 * source->psfMag) / flux0);
+        if (!fakeModel) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate model for source %d (%.1f,%.1f)", i, x, y);
+            return false;
+        }
+        if (circularise && !circulariseModel(fakeModel)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
+            psFree(fakeModel);
+            return false;
+        }
+#else
+        fakeModel->params->data.F32[PM_PAR_I0] = powf(10.0, -0.4 * source->psfMag) / flux0;
         fakeModel->params->data.F32[PM_PAR_XPOS] = x;
         fakeModel->params->data.F32[PM_PAR_YPOS] = y;
-        fakeModel->params->data.F32[PM_PAR_I0] = powf(10.0, -0.4 * source->psfMag) / flux0;
+#endif
 
         psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
@@ -115,6 +169,11 @@
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
         fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
-        float fakeRadius = radius > 0 ? radius : fakeModel->modelRadius(fakeModel->params, minFlux); // Radius
-
+        float fakeRadius =
+#ifdef PAP_WORK
+            radius > 0 ? radius :
+#endif
+            fakeModel->modelRadius(fakeModel->params, minFlux); // Radius
+
+#ifdef PAP_WORK
         if (xOffset) {
             if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i],
@@ -132,5 +191,7 @@
                 return false;
             }
-        } else {
+        } else
+#endif
+        {
             if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
@@ -147,7 +208,8 @@
         }
         psFree(fakeSource);
-    }
-
-    psFree(fakeModel);
+#ifdef PAP_WORK
+        psFree(fakeModel);
+#endif
+    }
 
     return true;