From 938af19e562424847a28612f2e5196f12a33e48e Mon Sep 17 00:00:00 2001
From: Ger van Diepen <diepen@astron.nl>
Date: Mon, 27 Jul 2015 09:10:33 +0000
Subject: [PATCH] Task #8157: Merged trunk into branch

---
 CEP/DP3/DPPP/src/StationAdder.cc   | 54 +++++++++++++++++++++++-------
 CEP/DP3/DPPP/test/tStationAdder.cc | 18 ++++++++--
 2 files changed, 57 insertions(+), 15 deletions(-)

diff --git a/CEP/DP3/DPPP/src/StationAdder.cc b/CEP/DP3/DPPP/src/StationAdder.cc
index 7abe640edff..f159a9c555f 100644
--- a/CEP/DP3/DPPP/src/StationAdder.cc
+++ b/CEP/DP3/DPPP/src/StationAdder.cc
@@ -56,7 +56,7 @@ namespace LOFAR {
         itsMinNPoint    (parset.getUint  (prefix+"minpoints", 1)),
         itsMakeAutoCorr (parset.getBool  (prefix+"autocorr", false)),
         itsSumAutoCorr  (parset.getBool  (prefix+"sumauto", true)),
-        itsDoAverage    (parset.getBool  (prefix+"average", false)),
+        itsDoAverage    (parset.getBool  (prefix+"average", true)),
         itsUseWeight    (parset.getBool  (prefix+"useweights", true))
     {
     }
@@ -341,6 +341,10 @@ namespace LOFAR {
         for (uint k=0; k<nrfr; ++k) {
           frfPtr[k] = true;
         }
+        for (uint k=0; k<3; ++k) {
+          uvwPtr[k]  = 0.;
+        }
+        double uvwWghtSum = 0.;
         // Sum the baselines forming the new baselines.
         for (uint j=0; j<itsBufRows[i].size(); ++j) {
           // Get the baseline number to use.
@@ -351,7 +355,7 @@ namespace LOFAR {
             blnr = -blnr;
             useConj = true;
           }
-          blnr--;
+          blnr--;     // decrement because blnr+1 is stored in itsBufRows
           // Get pointers to the input baseline data.
           const Complex* inDataPtr = (itsBuf.getData().data() +
                                       blnr*nrcc);
@@ -361,15 +365,21 @@ namespace LOFAR {
                                       blnr*nrcc);
           const Bool*    inFrfPtr  = (itsBuf.getFullResFlags().data() +
                                       blnr*nrfr);
-          // Add the data and weights if not flagged.
+          const Double*  inUvwPtr  = (itsBuf.getUVW().data() +
+                                      blnr*3);
+          // Add the data, uvw, and weights if not flagged.
           // Write 4 loops to avoid having to test inside the loop.
           if (useConj) {
             if (itsUseWeight) {
               for (uint k=0; k<nrcc; ++k) {
                 if (!inFlagPtr[k]) {
                   npoints[k]++;
-                  dataPtr[k] += conj(inDataPtr[k]);
+                  dataPtr[k] += conj(inDataPtr[k]) * inWghtPtr[k];
                   wghtPtr[k] += inWghtPtr[k];
+                  for (int ui=0; ui<3; ++ui) {
+                    uvwPtr[ui] -= inUvwPtr[ui] * inWghtPtr[k];
+                  }
+                  uvwWghtSum += inWghtPtr[k];
                 }
               }
             } else {
@@ -378,6 +388,10 @@ namespace LOFAR {
                   npoints[k]++;
                   dataPtr[k] += conj(inDataPtr[k]);
                   wghtPtr[k] += 1.;
+                  for (int ui=0; ui<3; ++ui) {
+                    uvwPtr[ui] -= inUvwPtr[ui];
+                  }
+                  uvwWghtSum += 1;
                 }
               }
             }
@@ -386,8 +400,12 @@ namespace LOFAR {
               for (uint k=0; k<nrcc; ++k) {
                 if (!inFlagPtr[k]) {
                   npoints[k]++;
-                  dataPtr[k] += inDataPtr[k];
+                  dataPtr[k] += inDataPtr[k] * inWghtPtr[k];
                   wghtPtr[k] += inWghtPtr[k];
+                  for (int ui=0; ui<3; ++ui) {
+                    uvwPtr[ui] += inUvwPtr[ui] * inWghtPtr[k];
+                  }
+                  uvwWghtSum += inWghtPtr[k];
                 }
               }
             } else {
@@ -396,6 +414,10 @@ namespace LOFAR {
                   npoints[k]++;
                   dataPtr[k] += inDataPtr[k];
                   wghtPtr[k] += 1.;
+                  for (int ui=0; ui<3; ++ui) {
+                    uvwPtr[ui] += inUvwPtr[ui];
+                  }
+                  uvwWghtSum += 1;
                 }
               }
             }
@@ -417,14 +439,20 @@ namespace LOFAR {
             }
           }
         }
-        // Calculate the UVW coordinate of the new station.
-        uint blnr = nrOldBL + i;
-        Vector<Double> uvws = itsUVWCalc.getUVW (getInfo().getAnt1()[blnr],
-                                                 getInfo().getAnt2()[blnr],
-                                                 buf.getTime());
-        uvwPtr[0] = uvws[0];
-        uvwPtr[1] = uvws[1];
-        uvwPtr[2] = uvws[2];
+        // Average or calculate the UVW coordinate of the new station.
+        if (itsDoAverage) {
+          for (int ui=0; ui<3; ++ui) {
+            uvwPtr[ui] /= uvwWghtSum;
+          }
+        } else {
+          uint blnr = nrOldBL + i;
+          Vector<Double> uvws = itsUVWCalc.getUVW (getInfo().getAnt1()[blnr],
+                                                   getInfo().getAnt2()[blnr],
+                                                   buf.getTime());
+          uvwPtr[0] = uvws[0];
+          uvwPtr[1] = uvws[1];
+          uvwPtr[2] = uvws[2];
+        }
         dataPtr += nrcc;
         flagPtr += nrcc;
         wghtPtr += nrcc;
diff --git a/CEP/DP3/DPPP/test/tStationAdder.cc b/CEP/DP3/DPPP/test/tStationAdder.cc
index 69e4a1f1531..54ca01d482a 100644
--- a/CEP/DP3/DPPP/test/tStationAdder.cc
+++ b/CEP/DP3/DPPP/test/tStationAdder.cc
@@ -276,13 +276,25 @@ private:
   void addData (Cube<Complex>& to, const Cube<Complex>& from,
                 Cube<Float>& tow, const Cube<Float>& weights, int bl)
   {
-    to += from(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl));
+    Cube<Complex> tmp=from.copy();
+    Cube<Complex>::iterator tmpit=tmp.begin();
+    Cube<Float>::const_iterator weightit=weights.begin();
+    for (; tmpit!=to.end() && weightit!=weights.end(); tmpit++, weightit++) {
+      *tmpit *= *weightit;
+    }
+    to += tmp(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl));
     tow += weights(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl));
   }
   void addConjData (Cube<Complex>& to, const Cube<Complex>& from,
                     Cube<Float>& tow, const Cube<Float>& weights, int bl)
   {
-    to += conj(from(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl)));
+    Cube<Complex> tmp=from.copy();
+    Cube<Complex>::iterator tmpit=tmp.begin();
+    Cube<Float>::const_iterator weightit=weights.begin();
+    for (; tmpit!=to.end() && weightit!=weights.end(); tmpit++, weightit++) {
+      *tmpit *= *weightit;
+    }
+    to += conj(tmp(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl)));
     tow += weights(IPosition(3,0,0,bl), IPosition(3,to.nrow()-1,to.ncolumn()-1,bl));
   }
   virtual bool process (const DPBuffer& buf)
@@ -472,6 +484,7 @@ void test2(int ntime, int nbl, int nchan, int ncorr)
   parset.add ("stations",
               "{ns1:[rs01.s01, rs02.s01], ns2:[cs01.s02, cs01.s01]}");
   parset.add ("autocorr", "false");
+  parset.add ("average", "false");
   DPStep::ShPtr step2(new StationAdder(in, parset, ""));
   DPStep::ShPtr step3(new TestOutput2(ntime, nbl, nchan, ncorr));
   step1->setNextStep (step2);
@@ -488,6 +501,7 @@ void test3 (const string& stations)
   ParameterSet parset;
   parset.add ("stations", stations);
   parset.add ("autocorr", "true");
+  parset.add ("average", "false");
   DPStep::ShPtr step2(new StationAdder(in, parset, ""));
   DPStep::ShPtr step3(new ThrowStep());
   step1->setNextStep (step2);
-- 
GitLab