Index: /branches/pap/psphot/src/psphotFake.c
===================================================================
--- /branches/pap/psphot/src/psphotFake.c	(revision 25024)
+++ /branches/pap/psphot/src/psphotFake.c	(revision 25025)
@@ -7,4 +7,61 @@
 #define MIN_FLUX 0.1                    // Minimum flux for faint sources
 #endif
+
+
+// Calculate the limiting flux for an image
+//
+// We limit ourselves to calculating the peak flux in the smoothed image, which should be close (modulo
+// non-Gaussian PSF) to the limiting flux in the un-smoothed original image.
+static float fakeLimit(const pmReadout *ro,     // Readout of interest
+                       const psMetadata *recipe // psphot recipe
+                       );
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, NAN);
+    PM_ASSERT_READOUT_VARIANCE(ro, NAN);
+    PS_ASSERT_METADATA_NON_NULL(recipe, NAN);
+
+    float xFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_X");
+    float yFWHM = psMetadataLookupF32(NULL, recipe, "FWHM_Y");
+    if (!isfinite(xFWHM) || !isfinite(yFWHM)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find FWHM_X and FWHM_Y");
+        return NAN;
+    }
+
+    float smoothSigma  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrtf(2.0*log(2.0)));
+    float smoothNsigma = psMetadataLookupF32(NULL, recipe, "PEAKS_SMOOTH_NSIGMA");
+    if (!isfinite(smoothSigma) || !isfinite(smoothNsigma)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to set smoothing parameters");
+        return NAN;
+    }
+
+    psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing
+    psKernel *newCovar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
+    psFree(kernel);
+    float factor = psImageCovarianceFactor(newCovar); // Variance factor
+    psFree(newCovar);
+
+    psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);        // Random number generator
+    if (!psImageBackground(stats, NULL, ro->variance, ro->mask, maskVal, rng) ||
+        !isfinite(stats->robustMedian)) {
+        psError(PSPHOT_ERR_DATA, false, "Unable to determine mean variance");
+        psFree(stats);
+        psFree(rng);
+        return NAN;
+    }
+    psFree(rng);
+    float meanVar = stats->robustMedian; // Mean variance
+    psFree(stats);
+
+    float thresh = psMetadataLookupF32(NULL, recipe, "PEAKS_NSIGMA_LIMIT_2");
+    if (!isfinite(thresh)) {
+        psError(PSPHOT_ERR_CONFIG, false, "Unable to find threshold");
+        return NAN;
+    }
+
+    return thresh * sqrtf(meanVar * factor); // Limiting peak flux in smoothed image
+}