Index: /trunk/psModules/src/detrend/Makefile.am
===================================================================
--- /trunk/psModules/src/detrend/Makefile.am	(revision 6998)
+++ /trunk/psModules/src/detrend/Makefile.am	(revision 6999)
@@ -5,4 +5,5 @@
 libpsmoduledetrend_la_SOURCES  = \
 	pmFlatField.c \
+	pmFlatNormalize.c \
 	pmFringeStats.c \
 	pmMaskBadPixels.c \
@@ -13,4 +14,5 @@
 	pmFlatField.h \
 	pmFlatFieldErrors.h \
+	pmFlatNormalize.h \
 	pmFringeStats.h \
 	pmMaskBadPixelsErrors.h \
Index: /trunk/psModules/src/detrend/pmFlatNormalize.c
===================================================================
--- /trunk/psModules/src/detrend/pmFlatNormalize.c	(revision 6999)
+++ /trunk/psModules/src/detrend/pmFlatNormalize.c	(revision 6999)
@@ -0,0 +1,125 @@
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <pslib.h>
+
+#include "pmFlatNormalize.h"
+
+
+// Estimate the flat-field normalisation
+bool pmFlatNormalize(psVector *sourceFlux, // The source flux in each image; modified for return
+                     psVector *chipGains, // Initial guess of the chip gains; modified for return
+                     psImage *fluxLevels, // Fluxes for each integration (row) and chip (col); modified
+                     unsigned int maxIter, // Maximum number of iterations
+                     double tolerance   // Tolerance level before dying
+                    )
+{
+    int numSources = sourceFlux->n;     // Number of measurements made
+    int numChips = chipGains->n;        // Number of chips with which each measurement is made
+
+    // Sanity checks
+    assert(fluxLevels->numCols == numSources);
+    assert(fluxLevels->numRows == numChips);
+    assert(sourceFlux->type.type == PS_TYPE_F64);
+    assert(chipGains->type.type == PS_TYPE_F64);
+    assert(fluxLevels->type.type == PS_TYPE_F64);
+    assert(maxIter >= 1);
+    assert(tolerance > 0);
+
+    // Take the logarithms
+    psImage *fluxMask = psImageAlloc(numSources, numChips, PS_TYPE_U8); // Mask for bad measurements
+    psImageInit(fluxMask, 0);
+    psVector *gainMask = psVectorAlloc(numChips, PS_TYPE_U8); // Mask for bad gains
+    psVectorInit(gainMask, 0);
+    psVector *sourceMask = psVectorAlloc(numSources, PS_TYPE_U8); // Mask for bad integrations
+    psVectorInit(sourceMask, 0);
+    for (int i = 0; i < numChips; i++) {
+        if (isfinite(chipGains->data.F64[i]) && chipGains->data.F64[i] > 0) {
+            chipGains->data.F64[i] = log(chipGains->data.F64[i]);
+        } else {
+            // Blank out this chip
+            gainMask->data.U8[i] = 1;
+            chipGains->data.F64[i] = NAN;
+        }
+
+        for (int j = 0; j < numSources; j++) {
+            if (isfinite(fluxLevels->data.F64[j][i]) && fluxLevels->data.F64[j][i] > 0) {
+                fluxLevels->data.F64[j][i] = log(fluxLevels->data.F64[j][i]);
+            } else {
+                // Blank out this measurement
+                fluxMask->data.U8[j][i] = 1;
+                fluxLevels->data.F64[j][i] = NAN;
+            }
+        }
+    }
+    // Don't need to initialise sourceFlux, since that is changed immediately
+
+
+    double diff = INFINITY;             // Difference from previous iteration
+    for (int iter = 0; iter < maxIter && diff > tolerance; iter++) {
+        diff = 0.0;
+
+        // Improve on the fluxes
+        for (int i = 0; i < numSources; i++) {
+            if (sourceMask->data.U8[i]) {
+                continue;
+            }
+            double previous = sourceFlux->data.F64[i]; // Value from previous iteration
+            double sum = 0.0;           // Sum of F_ij - G_j
+            int number = 0;             // Number of chips contributing
+            for (int j = 0; j < numChips; j++) {
+                if (!gainMask->data.U8[j] && !fluxMask->data.U8[i][j]) {
+                    sum += fluxLevels->data.F64[i][j] - chipGains->data.F64[j];
+                    number++;
+                }
+            }
+            if (number > 0) {
+                sourceFlux->data.F64[i] = sum / (double)number;
+                diff += abs((sourceFlux->data.F64[i] - previous) / sourceFlux->data.F64[i]);
+            } else {
+                sourceMask->data.U8[i] = 1;
+                sourceFlux->data.F64[i] = NAN;
+            }
+        }
+
+        // Improve on the gains
+        for (int i = 0; i < numChips; i++) {
+            if (gainMask->data.U8[i]) {
+                continue;
+            }
+            double previous = chipGains->data.F64[i]; // Value from previous iteration
+            double sum = 0.0;           // Sum of F_ji - S_j
+            int number = 0;             // Numer of sources contributing
+            for (int j = 0; j < numSources; j++) {
+                if (!fluxMask->data.U8[j][i]) {
+                    sum += fluxLevels->data.F64[j][i] - sourceFlux->data.F64[j];
+                    number++;
+                }
+            }
+            if (number > 0) {
+                chipGains->data.F64[i] = sum / (double)number;
+                diff += abs((chipGains->data.F64[i] - previous) / chipGains->data.F64[i]);
+            } else {
+                gainMask->data.U8[i] = 1;
+                chipGains->data.F64[i] = NAN;
+            }
+        }
+    }
+    psFree(fluxMask);
+
+    // Un-log everything
+    for (int i = 0; i < numChips; i++) {
+        if (!gainMask->data.U8[i]) {
+            chipGains->data.F64[i] = exp(chipGains->data.F64[i]);
+        }
+    }
+    for (int i = 0; i < numSources; i++) {
+        if (!sourceMask->data.U8[i]) {
+            sourceFlux->data.F64[i] = exp(sourceFlux->data.F64[i]);
+        }
+    }
+    psFree(gainMask);
+    psFree(sourceMask);
+
+    return (diff <= tolerance);         // Did we converge?
+}
Index: /trunk/psModules/src/detrend/pmFlatNormalize.h
===================================================================
--- /trunk/psModules/src/detrend/pmFlatNormalize.h	(revision 6999)
+++ /trunk/psModules/src/detrend/pmFlatNormalize.h	(revision 6999)
@@ -0,0 +1,15 @@
+#ifndef PM_FLAT_NORMALIZE_H
+#define PM_FLAT_NORMALIZE_H
+
+
+// Normalise the flat-field measurements (f_ij = g_i s_j where f_ij is the flux recorded for chip i and
+// integration j, g_i is the gain for the i-th chip, s_j is the flux of the source in the j-th integration).
+bool pmFlatNormalize(psVector *sourceFlux, // The source flux in each image; modified for return
+                     psVector *chipGains, // Initial guess of the chip gains; modified for return
+                     psImage *fluxLevels, // Fluxes for each integration (row) and chip (col); modified
+                     unsigned int maxIter, // Maximum number of iterations
+                     double tolerance // Tolerance level before dying
+                    );
+
+
+#endif