Index: /branches/pap/psphot/src/psphotFake.c
===================================================================
--- /branches/pap/psphot/src/psphotFake.c	(revision 25308)
+++ /branches/pap/psphot/src/psphotFake.c	(revision 25309)
@@ -3,4 +3,6 @@
 #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
                     PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
+
+#define TESTING
 
 
@@ -180,4 +182,5 @@
 
     float smoothSigma = 0.5*(fwhmMajor + fwhmMinor) / (2.0*sqrtf(2.0*log(2.0))); // Gaussian smoothing sigma
+    int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
 
     psImageMaskType maskVal = psMetadataLookupImageMask(NULL, recipe, "MASK.PSPHOT"); // Value to mask
@@ -186,4 +189,23 @@
     psphotRemoveAllSources(realSources, recipe);
     //    psphotAddNoise(readout, realSources, recipe);
+
+
+
+    {
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                readout->image->data.F32[y][x] = psRandomGaussian(rng);
+                readout->variance->data.F32[y][x] = 1.0;
+                readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+            }
+        }
+        psFree(readout->covariance);
+        readout->covariance = NULL;
+        psFree(rng);
+    }
+
+
+
 
     float magLim;                       // Guess at limiting magnitude
@@ -193,4 +215,10 @@
         return false;
     }
+
+#ifdef TESTING
+    psphotSaveImage(NULL, readout->image, "orig_image.fits");
+    psphotSaveImage(NULL, readout->variance, "orig_variance.fits");
+    psphotSaveImage(NULL, readout->mask, "orig_mask.fits");
+#endif
 
     psImage *xFake = NULL, *yFake = NULL; // Coordinates of sources, each bin in a row
@@ -201,4 +229,10 @@
         return false;
     }
+
+#ifdef TESTING
+    psphotSaveImage(NULL, readout->image, "fake_image.fits");
+    psphotSaveImage(NULL, readout->variance, "fake_variance.fits");
+    psphotSaveImage(NULL, readout->mask, "fake_mask.fits");
+#endif
 
     // XXX Could speed this up significantly by only convolving the central pixels of each fake source
@@ -211,8 +245,11 @@
     }
 
+#ifdef TESTING
+    psphotSaveImage(NULL, readout->mask, "fake_sig.fits");
+#endif
+
     thresh *= thresh;                   // Significance image is actually significance-squared
 
     int numBins = magOffsets->n;                          // Number of bins
-    int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
     psVector *count = psVectorAlloc(numBins, PS_TYPE_S32); // Number of sources found in each bin
     psArray *fakeSources = psArrayAllocEmpty(numSources * numBins); // Fake sources
@@ -221,4 +258,5 @@
     for (int i = 0; i < numBins; i++) {
         int num = 0;               // Number found
+
         for (int j = 0; j < numSources; j++) {
             // Coordinates of interest
@@ -269,13 +307,22 @@
     psVector *magErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin
     psVector *magErrStdev = psVectorAlloc(numBins, PS_TYPE_F32); // Stdev of error in magnitude for each bin
+    psVector *magErrErrMean = psVectorAlloc(numBins, PS_TYPE_F32); // Mean error in magnitude for each bin
     psVector *magErr = psVectorAlloc(numSources, PS_TYPE_F32);   // Magnitude errors
+    psVector *magErrErr = psVectorAlloc(numSources, PS_TYPE_F32); // Error in magnitude error
     psVector *magMask = psVectorAlloc(numSources, PS_TYPE_VECTOR_MASK); // Mask for magnitude errors
-    psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics
+    psStats *magStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    psVector *robustMedian = psVectorAlloc(numBins, PS_TYPE_F32);
+    psVector *robustStdev = psVectorAlloc(numBins, PS_TYPE_F32);
     for (int i = 0, index = 0; i < numBins; i++) {
         psStatsInit(magStats);
         psVectorInit(magMask, 0);
 
+        psString name = NULL;
+        psStringAppend(&name, "fake_%d.dat", i);
+        FILE *file = fopen(name, "w");
+
         float magRef = magLim + magOffsets->data.F32[i]; // Reference magnitude
-        for (int j = 0; index < fakeEnd->data.S32[i]; j++, index++) {
+        int j;                          // Index
+        for (j = 0; index < fakeEnd->data.S32[i]; j++, index++) {
             pmSource *source = fakeSources->data[index]; // Source of interest
             if (!source || !isfinite(source->psfMag)) {
@@ -283,6 +330,12 @@
                 continue;
             }
+
+            fprintf(file, "%f %f %f %f %f\n", source->peak->xf, source->peak->yf,
+                    magRef, source->psfMag, source->errMag);
             magErr->data.F32[j] = source->psfMag - magRef;
-        }
+            magErrErr->data.F32[j] = source->errMag;
+        }
+        magErr->n = j;
+        magMask->n = j;
         if (!psVectorStats(magStats, magErr, NULL, magMask, 0xFF)) {
             // Probably because we don't have enough sources
@@ -291,5 +344,15 @@
         magErrMean->data.F32[i] = magStats->sampleMean;
         magErrStdev->data.F32[i] = magStats->sampleStdev;
-    }
+        robustMedian->data.F32[i] = magStats->robustMedian;
+        robustStdev->data.F32[i] = magStats->robustStdev;
+        if (!psVectorStats(magStats, magErrErr, NULL, magMask, 0xFF)) {
+            // Probably because we don't have enough sources
+            psErrorClear();
+        }
+        magErrErrMean->data.F32[i] = magStats->sampleMean;
+        fclose(file);
+    }
+
+
     psFree(magStats);
     psFree(magMask);
@@ -313,4 +376,10 @@
     psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RMS", PS_META_REPLACE,
                         "Stdev in magnitude differences", magErrStdev);
+    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RM", PS_META_REPLACE,
+                        "Robust median magnitude differences", robustMedian);
+    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.RS", PS_META_REPLACE,
+                        "Robust stdev in magnitude differences", robustStdev);
+    psMetadataAddVector(stats, PS_LIST_TAIL, "FAKE.MEANERR", PS_META_REPLACE,
+                        "Mean error in magnitude differences", magErrErrMean);
     psMetadataConfigWrite(stats, "fake.stats");
     psFree(stats);
@@ -320,4 +389,7 @@
     psFree(magErrStdev);
 
+    psFree(robustMedian);
+    psFree(robustStdev);
+
     return true;
 }