Index: /branches/pap/psphot/src/psphotFake.c
===================================================================
--- /branches/pap/psphot/src/psphotFake.c	(revision 25277)
+++ /branches/pap/psphot/src/psphotFake.c	(revision 25278)
@@ -1,10 +1,6 @@
 # include "psphotInternal.h"
 
-#define TESTING
-
-
-#ifdef TESTING
-#define MIN_FLUX 0.1                    // Minimum flux for faint sources
-#endif
+#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
+                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
 
 
@@ -14,10 +10,11 @@
 // non-Gaussian PSF) to the limiting flux in the un-smoothed original image.
 static bool fakeLimit(float *magLim,           // Limiting magntiude, to return
-                      float *minFlux,          // Minimum flux, to return
+                      int *radius,             // Radius for fake sources
                       const pmReadout *ro,     // Readout of interest
+                      const pmPSF *psf,        // Point-spread function
                       float thresh,            // Threshold for source identification
-                      float fwhmMajor, float fwhmMinor, // Size of PSF
-                      float smoothNsigma,       // Smoothing limit
-                      psImageMaskType maskVal   // Value to mask
+                      float smoothSigma,       // Gaussian smoothing sigma
+                      float smoothNsigma,      // Smoothing limit
+                      psImageMaskType maskVal  // Value to mask
                       )
 {
@@ -25,5 +22,4 @@
     PM_ASSERT_READOUT_VARIANCE(ro, false);
 
-    float smoothSigma  = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0)));
     psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
     psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix
@@ -45,9 +41,34 @@
     psFree(stats);
 
+    // Need to normalise out difference between Gaussian and real PSF
+    pmModel *normModel = pmModelFromPSFforXY(psf, 0.5 * ro->variance->numCols,
+                                             0.5 * ro->variance->numRows, 1.0); // Model for normalisation
+    if (!normModel || (normModel->flags & MODEL_MASK)) {
+        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate PSF model.");
+        psFree(normModel);
+        return false;
+    }
+    float norm = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0
+    psFree(normModel);
+
     if (magLim) {
-        *magLim = -2.5 * log10(thresh * sqrtf(meanVar * factor));
-    }
-    if (minFlux) {
-        *minFlux = sqrtf(meanVar);
+        // The signal-to-noise of the smoothed image is: S/N ~ I_smooth / sqrt(variance * factor)
+
+        // since the variance factor tells us the variance in the smoothed image.  Now, the trick is working
+        // out what the intensity in the smoothed image is, and how it is related to the flux.  We are
+        // convolving Io G_* with G_PSF/2pi.w^2, where G_* is the approximately-Gaussian PSF of the star,
+        // G_PSF is the pretty-close-matching Gaussian of the convolution kernel, Io is the peak flux in the
+        // unsmoothed image, and w is the width of the Gaussian.  Now, a normalised 2D Gaussian convolved
+        // with itself has a normalisation of 1/2pi(w^2+w^2) = 1/4pi.w^2.  Therefore:
+        // I_smooth = Flux / 2*norm.
+        float peakLim = thresh * sqrtf(meanVar * factor); // Limiting peak value in smoothed image
+        float fluxLim = 2.0 * norm * peakLim; // Limiting flux in original
+        *magLim = -2.5 * log10f(fluxLim);
+        psTrace("psphot.fake", 1, "Limiting peak: %f\n", peakLim);
+        psTrace("psphot.fake", 1, "Limiting flux: %f\n", fluxLim);
+        psTrace("psphot.fake", 1, "Limiting mag: %f\n", *magLim);
+    }
+    if (radius) {
+        *radius = smoothSigma * smoothNsigma;
     }
 
@@ -62,5 +83,5 @@
                          int numSources,                 // Number of fake sources for each bin
                          float refMag,                   // Reference magnitude
-                         float minFlux                   // Minimum flux level for fake image
+                         int radius                      // Radius for fake sources
                          )
 {
@@ -99,5 +120,5 @@
     pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
     if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll,
-                                  NULL, NULL, psf, minFlux, 0, false, false)) {
+                                  NULL, NULL, psf, NAN, radius, false, true)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image");
         psFree(fakeRO);
@@ -130,5 +151,4 @@
 
     // Collect recipe information
-    bool mdok;                                                 // Status of MD lookup
     float fwhmMajor = psMetadataLookupF32(NULL, recipe, "FWHM_MAJ"); // PSF size in x
     float fwhmMinor = psMetadataLookupF32(NULL, recipe, "FWHM_MIN"); // PSF size in y
@@ -157,8 +177,6 @@
         return NULL;
     }
-    float minFluxFudge = psMetadataLookupF32(&mdok, recipe, "FAKE.MINFLUX"); // Fudge factor for minimum flux
-    if (!mdok) {
-        minFluxFudge = 1.0;
-    }
+
+    float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma
 
     psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
@@ -166,16 +184,15 @@
     // remove all sources, adding noise for subtracted sources
     psphotRemoveAllSources(realSources, recipe);
-    psphotAddNoise(readout, realSources, recipe);
+    //    psphotAddNoise(readout, realSources, recipe);
 
     float magLim;                       // Guess at limiting magnitude
-    float minFlux;                      // Minimum flux for fake image
-    if (!fakeLimit(&magLim, &minFlux, readout, thresh, fwhmMajor, fwhmMinor, smoothNsigma, maskVal)) {
+    int radius;                         // Radius for fake sources
+    if (!fakeLimit(&magLim, &radius, readout, psf, thresh, smoothSigma, smoothNsigma, maskVal)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to determine limits for image");
         return false;
     }
-    minFlux *= minFluxFudge;
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
-    if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, minFlux)) {
+    if (!fakeGenerate(&xFake, &yFake, readout, psf, magOffsets, numSources, magLim, radius)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources");
         psFree(xFake);
@@ -195,10 +212,13 @@
 
     int numBins = magOffsets->n;                          // Number of bins
+    int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
     psVector *frac = psVectorAlloc(numBins, PS_TYPE_F32); // Fraction of sources in each bin
     for (int i = 0; i < numBins; i++) {
         int numFound = 0;               // Number found
         for (int j = 0; j < numSources; j++) {
-            float x = xFake->data.F32[i][j], y = yFake->data.F32[i][j]; // Coordinates of interest
-            int xPix = x, yPix = y;                                     // Pixel coordinates
+            // Coordinates of interest
+            float x = PS_MIN(xFake->data.F32[i][j] + 0.5, numCols - 1);
+            float y = PS_MIN(yFake->data.F32[i][j] + 0.5, numRows - 1);
+            int xPix = x, yPix = y;     // Pixel coordinates
             if (significance->data.F32[yPix][xPix] > thresh) {
                 numFound++;