Changeset 7178

Timestamp:

May 22, 2006, 5:17:19 PM (20 years ago)

Author:

Paul Price

Message:

Normalising source fluxes by the mean

File:

• trunk/psModules/src/detrend/pmFlatNormalize.c (modified) (6 diffs)

trunk/psModules/src/detrend/pmFlatNormalize.c

@@ -26 +26 @@
     // Take the logarithms
     psImage *flux = psImageCopy(NULL, fluxLevels, PS_TYPE_F64); // Copy of the input flux levels matrix
-    psImage *fluxMask = psImageAlloc(numSources, numChips, PS_TYPE_U8); // Mask for bad measurements
+    psImage *fluxMask = psImageAlloc(numChips, numSources, 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);
     for (int i = 0; i < numChips; i++) {
@@ -36 +38 @@
             chipGains->data.F64[i] = log(chipGains->data.F64[i]);
         } else {
+            chipGains->data.F64[i] = 0.0; // Take a wild guess
+            #if 0
             // Blank out this chip
             gainMask->data.U8[i] = 1;
             chipGains->data.F64[i] = NAN;
+            #endif
+
         }
 
@@ -52 +58 @@
     }
 
-
     double diff = INFINITY;             // Difference from previous iteration
     psVector *sourceFlux = psVectorAlloc(numSources, PS_TYPE_F64); // The flux in each integration
-    // Don't need to initialise sourceFlux, since that is changed immediately
+    sourceFlux->n = numSources;
+    psVector *oldSourceFlux = psVectorAlloc(numSources, PS_TYPE_F64); // 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
+        double sumFlux = 0.0;           // Total fluxes
+        long numFluxes = 0;             // Number of fluxes
         for (int i = 0; i < numSources; i++) {
             if (sourceMask->data.U8[i]) {
+                psTrace(__func__, 7, "Flux for exposure %d is masked.\n", i);
                 continue;
             }
-            double previous = sourceFlux->data.F64[i]; // Value from previous iteration
+            numFluxes++;
             double sum = 0.0;           // Sum of F_ij - G_j
             int number = 0;             // Number of chips contributing
@@ -75 +87 @@
             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;
             }
+            sumFlux += exp(sourceFlux->data.F64[i]);
+            psTrace(__func__, 7, "Flux for exposure %d is %f\n", i, exp(sourceFlux->data.F64[i]));
+        }
+        // Normalise the mean to unity
+        sumFlux /= (double)numFluxes;
+        sumFlux = log(sumFlux);
+        for (int i = 0; i < numSources; i++) {
+            if (sourceMask->data.U8[i]) {
+                continue;
+            }
+            sourceFlux->data.F64[i] -= sumFlux;
+            diff += abs((sourceFlux->data.F64[i] - oldSourceFlux->data.F64[i]) / sourceFlux->data.F64[i]);
         }
 
@@ -87 +110 @@
                 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
@@ -98 +120 @@
             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;
             }
+            psTrace(__func__, 7, "Gain for chip %d is %f\n", i, exp(chipGains->data.F64[i]));
         }
+
+        psTrace(__func__, 2, "Iteration %d: difference is %f\n", iter, diff);
+
+        // Switch the old and new
+        psVector *temp = sourceFlux;
+        sourceFlux = oldSourceFlux;
+        oldSourceFlux = temp;
     }
     psFree(flux);
     psFree(fluxMask);
+    psFree(oldSourceFlux);
 
     // Un-log the vectors