Index: trunk/psModules/src/detrend/pmFlatNormalize.c
===================================================================
--- trunk/psModules/src/detrend/pmFlatNormalize.c	(revision 8246)
+++ trunk/psModules/src/detrend/pmFlatNormalize.c	(revision 8397)
@@ -6,46 +6,76 @@
 #include "pmFlatNormalize.h"
 
+// I'm not sure that many many iterations are required, but rather suspect that the system converges within a
+// few with absolutely no trouble (it *is* over-constrained).  For this reason, I'm putting the maximum number
+// of iterations and tolerance as preset values.
+#define MAXITER 10
+#define TOLERANCE 1e-3
 
 // Estimate the flat-field normalisation; return the source flux in each integration
-psVector *pmFlatNormalize(bool *converge, // Did we converge?
-                          psVector *chipGains, // Initial guess of the chip gains; modified for return
-                          const psImage *fluxLevels, // Fluxes for each integration (row) and chip (col)
-                          unsigned int maxIter, // Maximum number of iterations
-                          float tolerance   // Tolerance level before dying
-                         )
+bool pmFlatNormalize(psVector **expFluxesPtr, // Flux in each exposure; modified for return
+                     psVector **chipGainsPtr, // Initial guess of the chip gains; modified for return
+                     const psImage *bgMatrix
+                    )
 {
-    PS_ASSERT_PTR_NON_NULL(chipGains, NULL);
-    PS_ASSERT_PTR_NON_NULL(fluxLevels, NULL);
+    PS_ASSERT_PTR_NON_NULL(bgMatrix, false);
+    PS_ASSERT_IMAGE_NON_NULL(bgMatrix, false);
 
-    int numSources = fluxLevels->numRows; // Number of integrations
-    int numChips = fluxLevels->numCols; // Number of chips with which each integration is made
+    int numExps = bgMatrix->numRows; // Number of exposures
+    int numChips = bgMatrix->numCols; // Number of chips with which each exposure is made
 
-    // Sanity checks
-    assert(chipGains->n == numChips);
-    assert(chipGains->type.type == PS_TYPE_F32);
-    assert(maxIter >= 1);
-    assert(tolerance > 0);
+    psVector *expFluxes;                // Dereferenced version of expFluxesPtr
+    if (expFluxesPtr) {
+        if (*expFluxesPtr) {
+            PS_ASSERT_VECTOR_TYPE(*expFluxesPtr, PS_TYPE_F32, false);
+            PS_ASSERT_VECTOR_SIZE(*expFluxesPtr, numExps, false);
+        } else {
+            *expFluxesPtr = psVectorAlloc(numExps, PS_TYPE_F32);
+            (*expFluxesPtr)->n = numExps;
+        }
+        expFluxes = psMemIncrRefCounter(*expFluxesPtr);
+    } else {
+        expFluxes = psVectorAlloc(numExps, PS_TYPE_F32);
+        expFluxes->n = numExps;
+    }
+
+    psVector *chipGains;                // Dereferenced version of chipGainsPtr
+    if (chipGainsPtr) {
+        if (*chipGainsPtr) {
+            PS_ASSERT_VECTOR_TYPE(*chipGainsPtr, PS_TYPE_F32, false);
+            PS_ASSERT_VECTOR_SIZE(*chipGainsPtr, numChips, false);
+        } else {
+            *chipGainsPtr = psVectorAlloc(numChips, PS_TYPE_F32);
+            (*chipGainsPtr)->n = numChips;
+            psVectorInit(*chipGainsPtr, 1.0);
+        }
+        chipGains = psMemIncrRefCounter(*chipGainsPtr);
+    } else {
+        chipGains = psVectorAlloc(numChips, PS_TYPE_F32);
+        chipGains->n = numChips;
+        psVectorInit(chipGains, 1.0);
+    }
 
     // Take the logarithms
-    psImage *flux = psImageCopy(NULL, fluxLevels, PS_TYPE_F32); // Copy of the input flux levels matrix
-    psImage *fluxMask = psImageAlloc(numChips, numSources, PS_TYPE_U8); // Mask for bad measurements
+    psImage *flux = psImageCopy(NULL, bgMatrix, PS_TYPE_F32); // Copy of the input flux levels matrix
+    psImage *fluxMask = psImageAlloc(numChips, numExps, PS_TYPE_U8); // Mask for bad measurements
     psImageInit(fluxMask, 0);
     psVector *gainMask = psVectorAlloc(numChips, PS_TYPE_U8); // Mask for bad gains
     gainMask->n = numChips;
     psVectorInit(gainMask, 0);
-    psVector *sourceMask = psVectorAlloc(numSources, PS_TYPE_U8); // Mask for bad integrations
-    sourceMask->n = numSources;
-    psVectorInit(sourceMask, 0);
+    psVector *expMask = psVectorAlloc(numExps, PS_TYPE_U8); // Mask for bad exposures
+    expMask->n = numExps;
+    psVectorInit(expMask, 0);
     for (int i = 0; i < numChips; i++) {
         // Note: the input gains are in e/ADU; we want to work with ADU/e (bg [ADU] = g [ADU/e] * f [e])
+        // Hence the minus sign
         if (isfinite(chipGains->data.F32[i]) && chipGains->data.F32[i] > 0) {
-            chipGains->data.F32[i] = -log(chipGains->data.F32[i]);
+            chipGains->data.F32[i] = -logf(chipGains->data.F32[i]);
         } else {
-            chipGains->data.F32[i] = 0.0; // Take a wild guess
+            chipGains->data.F32[i] = 0.0; // Take a wild guess, gain ~ 1 e/ADU
         }
 
-        for (int j = 0; j < numSources; j++) {
+        for (int j = 0; j < numExps; j++) {
             if (isfinite(flux->data.F32[j][i]) && flux->data.F32[j][i] > 0) {
-                flux->data.F32[j][i] = log(flux->data.F32[j][i]);
+                flux->data.F32[j][i] = logf(flux->data.F32[j][i]);
             } else {
                 // Blank out this measurement
@@ -58,22 +88,17 @@
     // Not really sure that we need to iterate, but here we go anyway...
 
-    float diff = INFINITY;             // Difference from previous iteration
-    psVector *sourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The flux in each integration
-    sourceFlux->n = numSources;
-    psVector *oldSourceFlux = psVectorAlloc(numSources, PS_TYPE_F32); // The fluxes in the previous iteration
-    oldSourceFlux->n = numSources;
-    psVectorInit(oldSourceFlux, 0.0);
-    for (int iter = 0; iter < maxIter && diff > tolerance; iter++) {
-        diff = 0.0;
-
-        // Improve on the fluxes
-        long numFluxes = 0;             // Number of fluxes
-        for (int i = 0; i < numSources; i++) {
-            if (sourceMask->data.U8[i]) {
+    float diff = INFINITY;              // Difference from previous iteration
+    psVector *oldExpFluxes = NULL;      // The fluxes in the previous iteration
+    psVector *oldChipGains = NULL;      // Chip gains in the previous iteration
+    for (int iter = 0; iter < MAXITER && diff > TOLERANCE; iter++) {
+        // Improve on the exposure fluxes
+        int numFluxes = 0;              // Number of fluxes
+        for (int i = 0; i < numExps; i++) {
+            if (expMask->data.U8[i]) {
                 psTrace("psModules.detrend", 7, "Flux for exposure %d is masked.\n", i);
                 continue;
             }
             numFluxes++;
-            float sum = 0.0;           // Sum of F_ij - G_j
+            float sum = 0.0;            // Sum of F_ij - G_j
             int number = 0;             // Number of chips contributing
             for (int j = 0; j < numChips; j++) {
@@ -84,19 +109,15 @@
             }
             if (number > 0) {
-                sourceFlux->data.F32[i] = sum / (float)number;
+                expFluxes->data.F32[i] = sum / (float)number;
             } else {
-                sourceMask->data.U8[i] = 1;
-                sourceFlux->data.F32[i] = NAN;
+                expMask->data.U8[i] = 1;
+                expFluxes->data.F32[i] = NAN;
             }
-            psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, exp(sourceFlux->data.F32[i]));
-        }
-
-        for (int i = 0; i < numSources; i++) {
-            if (!sourceMask->data.U8[i]) {
-                diff += abs((sourceFlux->data.F32[i] - oldSourceFlux->data.F32[i]) / sourceFlux->data.F32[i]);
-            }
+            psTrace("psModules.detrend", 7, "Flux for exposure %d is %lf\n", i, expf(expFluxes->data.F32[i]));
         }
 
         // Improve on the gains
+        float meanGain = 0.0;           // Mean gain
+        int numGains = 0;               // Number of gains
         for (int i = 0; i < numChips; i++) {
             if (gainMask->data.U8[i]) {
@@ -105,7 +126,7 @@
             float sum = 0.0;           // Sum of F_ji - S_j
             int number = 0;             // Numer of sources contributing
-            for (int j = 0; j < numSources; j++) {
+            for (int j = 0; j < numExps; j++) {
                 if (!fluxMask->data.U8[j][i]) {
-                    sum += flux->data.F32[j][i] - sourceFlux->data.F32[j];
+                    sum += flux->data.F32[j][i] - expFluxes->data.F32[j];
                     number++;
                 }
@@ -117,36 +138,57 @@
                 chipGains->data.F32[i] = NAN;
             }
-            psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, exp(-chipGains->data.F32[i]));
+            psTrace("psModules.detrend", 7, "Gain for chip %d is %lf\n", i, expf(-chipGains->data.F32[i]));
+            meanGain += expf(chipGains->data.F32[i]);
+            numGains++;
+        }
+
+        // Normalise the mean gain to unity, and measure the difference
+        meanGain /= (float)numGains;
+        meanGain = logf(meanGain);
+        if (iter > 0) {
+            diff = 0.0;
+            for (int i = 0; i < numChips; i++) {
+                if (gainMask->data.U8[i]) {
+                    continue;
+                }
+                chipGains->data.F32[i] -= meanGain;
+                diff += abs((chipGains->data.F32[i] - oldChipGains->data.F32[i]) / chipGains->data.F32[i]);
+            }
+            for (int i = 0; i < numExps; i++) {
+                if (expMask->data.U8[i]) {
+                    continue;
+                }
+                diff += abs((expFluxes->data.F32[i] - oldExpFluxes->data.F32[i]) / expFluxes->data.F32[i]);
+            }
         }
 
         psTrace("psModules.detrend", 2, "Iteration %d: difference is %e\n", iter, diff);
 
-        // Switch the old and new
-        psVector *temp = sourceFlux;
-        sourceFlux = oldSourceFlux;
-        oldSourceFlux = temp;
+        // Copy the new to the old
+        oldChipGains = psVectorCopy(oldChipGains, chipGains, PS_TYPE_F32);
+        oldExpFluxes = psVectorCopy(oldExpFluxes, expFluxes, PS_TYPE_F32);
     }
     psFree(flux);
     psFree(fluxMask);
-    psFree(oldSourceFlux);
+    psFree(oldChipGains);
+    psFree(oldExpFluxes);
 
     // Un-log the vectors
     for (int i = 0; i < numChips; i++) {
         if (!gainMask->data.U8[i]) {
-            chipGains->data.F32[i] = exp(chipGains->data.F32[i]);
+            chipGains->data.F32[i] = expf(chipGains->data.F32[i]);
         }
     }
-    for (int i = 0; i < numSources; i++) {
-        if (!sourceMask->data.U8[i]) {
-            sourceFlux->data.F32[i] = exp(sourceFlux->data.F32[i]);
+    for (int i = 0; i < numExps; i++) {
+        if (!expMask->data.U8[i]) {
+            expFluxes->data.F32[i] = expf(expFluxes->data.F32[i]);
         }
     }
     psFree(gainMask);
-    psFree(sourceMask);
+    psFree(expMask);
 
-    if (converge) {
-        *converge = (diff < tolerance); // Did we converge?
-    }
+    psFree(chipGains);
+    psFree(expFluxes);
 
-    return sourceFlux;
+    return (diff < TOLERANCE); // Did we converge?
 }