Changeset 518 for trunk/src

Timestamp:

02/28/05 21:48:50 (20 years ago)

Author:

kil064

Message:

add Tsys and Tint weightring

Location:

trunk/src

Files:

• SDMath.cc (modified) (13 diffs)
• SDMath.h (modified) (4 diffs)

trunk/src/SDMath.cc

@@ -150 +150 @@
   
   WeightType wtType = NONE;
-  convertWeightString(wtType, weightStr);
+  convertWeightString(wtType, weightStr, True);
 
 // Create output Table by cloning from the first table
@@ -219 +219 @@
 // To get the right shape for the Tsys accumulator we need to
 // access a column from the first table.  The shape of this
-// array must not change
-
-  Array<Float> tSysSum;
+// array must not change.  Note however that since the TSysSqSum
+// array is used in a normalization process, and that I ignore the
+// channel axis replication of values for now, it loses a dimension
+
+  Array<Float> tSysSum, tSysSqSum;
   {
     const Table& tabIn = in[0]->table();
     tSysCol.attach(tabIn,"TSYS");
     tSysSum.resize(tSysCol.shape(0));
+//
+    tSysSqSum.resize(shp2);
   }
   tSysSum =0.0;
+  tSysSqSum = 0.0;
   Array<Float> tSys;
 
@@ -306 +311 @@
 // Normalize data in 'sum' accumulation array according to weighting scheme
 
-           normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
+           normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
 
 // Get ScanContainer for the first row of this averaged Scan
@@ -331 +336 @@
 //
            tSysSum =0.0;
+           tSysSqSum =0.0;
            timeSum = 0.0;
            intSum = 0.0;
@@ -348 +354 @@
 // Accumulate
 
-        accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, 
-                    tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis, 
-                    nAxesSub, useMask, wtType);
+        accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum, tSysSqSum,
+                   tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis, 
+                   nAxesSub, useMask, wtType);
 //
        oldSourceName = sourceName;
@@ -364 +370 @@
 // Normalize data in 'sum' accumulation array according to weighting scheme
 
-  normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
+  normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType, asap::ChanAxis, nAxesSub);
 
 // Create and fill container.  The container we clone will be from
@@ -807 +813 @@
 {
    WeightType wtType = NONE;
-   convertWeightString(wtType, weightStr);
+   convertWeightString(wtType, weightStr, False);
 
    const uInt nRows = in.nRow();
@@ -1232 +1238 @@
 //
 // phase in degrees
-// Applies to all Beams and IFs
-// Might want to optionally select on Beam/IF
+// assumes linear correlations
 //
 {
    if (in.nPol() != 4) {
       throw(AipsError("You must have 4 polarizations to run this function"));
+   }
+//
+   SDHeader sh = in.getSDHeader();
+   Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
+   SDAttr sdAtt;
+   if (sdAtt.feedPolType(inst) != LINEAR) {
+      throw(AipsError("Only linear polarizations are supported"));
    }
 //    
@@ -1284 +1296 @@
 //
 // phase in degrees
-// Applies to all Beams and IFs
-// Might want to optionally select on Beam/IF
+// assumes linear correlations
 //
 {
    if (in.nPol() != 4) {
       throw(AipsError("You must have 4 polarizations to run this function"));
+   }
+//
+   SDHeader sh = in.getSDHeader();
+   Instrument inst = SDAttr::convertInstrument (sh.antennaname, False);
+   SDAttr sdAtt;
+   if (sdAtt.feedPolType(inst) != LINEAR) {
+      throw(AipsError("Only linear polarizations are supported"));
    }
 //    
@@ -1700 +1718 @@
 }
 
+void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
+                        MaskedArray<Float>& sum, Array<Float>& sumSq, 
+                        Array<Float>& nPts, Array<Float>& tSysSum, 
+                        Array<Float>& tSysSqSum,
+                        const Array<Float>& tSys, const Array<Float>& nInc, 
+                        const Vector<Bool>& mask, Double time, Double interval,
+                        const Block<CountedPtr<SDMemTable> >& in,
+                        uInt iTab, uInt iRow, uInt axis, 
+                        uInt nAxesSub, Bool useMask,
+                        WeightType wtType) const
+{
+
+// Get data
+
+   MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
+   Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
+   const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
+//
+   if (wtType==NONE) {
+      const MaskedArray<Float> n(nInc,dataIn.getMask());
+      nPts += n;                               // Only accumulates where mask==T
+   } else if (wtType==TINT) {
+
+// We are weighting the data by integration time.
+
+     valuesIn *= Float(interval);
+
+   } else if (wtType==VAR) {
+
+// We are going to average the data, weighted by the noise for each pol, beam and IF.
+// So therefore we need to iterate through by spectrum (axis 3)
+
+      VectorIterator<Float> itData(valuesIn, axis);
+      ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
+      Float fac = 1.0;
+      IPosition pos(nAxesSub,0);  
+//
+      while (!itData.pastEnd()) {
+
+// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
+
+         if (useMask) {
+            MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
+            fac = 1.0/variance(tmp);
+         } else {
+            MaskedArray<Float> tmp(itData.vector(),itMask.vector());
+            fac = 1.0/variance(tmp);
+         }
+
+// Scale data
+
+         itData.vector() *= fac;     // Writes back into 'dataIn'
+//
+// Accumulate variance per if/pol/beam averaged over spectrum
+// This method to get pos2 from itData.pos() is only valid
+// because the spectral axis is the last one (so we can just
+// copy the first nAXesSub positions out)
+
+         pos = itData.pos().getFirst(nAxesSub);
+         sumSq(pos) += fac;
+//
+         itData.next();
+         itMask.next();
+      }
+   } else if (wtType==TSYS) {
+
+// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
+// So therefore we need to iterate through by spectrum (axis 3).  Although
+// Tsys is stored as a vector of length nChan, the values are replicated.
+// We will take a short cut and just use the value from the first channel
+// for now.
+// 
+      VectorIterator<Float> itData(valuesIn, axis);
+      ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
+      Float fac = 1.0;
+      IPosition pos(nAxesSub,0);  
+//
+      while (!itData.pastEnd()) {
+         Float t = itTSys.vector()[0];
+         fac = 1.0/t/t;
+
+// Scale data
+
+         itData.vector() *= fac;     // Writes back into 'dataIn'
+//
+// Accumulate Tsys  per if/pol/beam averaged over spectrum
+// This method to get pos2 from itData.pos() is only valid
+// because the spectral axis is the last one (so we can just
+// copy the first nAXesSub positions out)
+
+         pos = itData.pos().getFirst(nAxesSub);
+         tSysSqSum(pos) += fac;
+//
+         itData.next();
+         itTSys.next();
+      }
+   }
+
+// Accumulate sum of (possibly scaled) data
+
+   sum += dataIn;
+
+// Accumulate Tsys, time, and interval
+
+   tSysSum += tSys;
+   timeSum += time;
+   intSum += interval;
+   nAccum += 1;
+}
+
+
 void SDMath::normalize(MaskedArray<Float>& sum,
-                        const Array<Float>& sumSq,
-                        const Array<Float>& nPts,
-                        WeightType wtType, Int axis,
-                        Int nAxesSub) const
+                       const Array<Float>& sumSq,
+                       const Array<Float>& tSysSqSum,
+                       const Array<Float>& nPts,
+                       Double intSum,
+                       WeightType wtType, Int axis,
+                       Int nAxesSub) const
 {
    IPosition pos2(nAxesSub,0);
@@ -1716 +1847 @@
       MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
       sum /= t;
+   } else if (wtType==TINT) {
+
+// Average by sum of Tint
+
+      sum /= Float(intSum);
    } else if (wtType==VAR) {
 
@@ -1729 +1865 @@
       }
    } else if (wtType==TSYS) {
-   }
-}
-
-
-void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
-                        MaskedArray<Float>& sum, Array<Float>& sumSq, 
-                        Array<Float>& nPts, Array<Float>& tSysSum, 
-                        const Array<Float>& tSys, const Array<Float>& nInc, 
-                        const Vector<Bool>& mask, Double time, Double interval,
-                        const Block<CountedPtr<SDMemTable> >& in,
-                        uInt iTab, uInt iRow, uInt axis, 
-                        uInt nAxesSub, Bool useMask,
-                        WeightType wtType) const
-{
-
-// Get data
-
-   MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
-   Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
-   const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
-//
-   if (wtType==NONE) {
-      const MaskedArray<Float> n(nInc,dataIn.getMask());
-      nPts += n;                               // Only accumulates where mask==T
-   } else if (wtType==VAR) {
-
-// We are going to average the data, weighted by the noise for each pol, beam and IF.
-// So therefore we need to iterate through by spectrum (axis 3)
-
-      VectorIterator<Float> itData(valuesIn, axis);
-      ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
-      Float fac = 1.0;
-      IPosition pos(nAxesSub,0);  
-//
+   
+// Normalize each spectrum by sum(1/Tsys**2) where the pseudo
+// replication over channel for Tsys has been dropped.
+
+      Array<Float>& data = sum.getRWArray();
+      VectorIterator<Float> itData(data, axis);
       while (!itData.pastEnd()) {
-
-// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
-
-        if (useMask) {
-           MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
-           fac = 1.0/variance(tmp);
-        } else {
-           MaskedArray<Float> tmp(itData.vector(),itMask.vector());
-           fac = 1.0/variance(tmp);
-        }
-
-// Scale data
-
-        itData.vector() *= fac;     // Writes back into 'dataIn'
-//
-// Accumulate variance per if/pol/beam averaged over spectrum
-// This method to get pos2 from itData.pos() is only valid
-// because the spectral axis is the last one (so we can just
-// copy the first nAXesSub positions out)
-
-        pos = itData.pos().getFirst(nAxesSub);
-        sumSq(pos) += fac;
-//
-        itData.next();
-        itMask.next();
+         pos2 = itData.pos().getFirst(nAxesSub);
+         itData.vector() /= tSysSqSum(pos2);
+         itData.next();
       }
-   } else if (wtType==TSYS) {
-   }
-
-// Accumulate sum of (possibly scaled) data
-
-   sum += dataIn;
-
-// Accumulate Tsys, time, and interval
-
-   tSysSum += tSys;
-   timeSum += time;
-   intSum += interval;
-   nAccum += 1;
-}
+   }
+}
+
+
 
 
@@ -1835 +1911 @@
 
 
-void SDMath::convertWeightString(WeightType& wtType, const String& weightStr) const
+void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
+                                 Bool listType) const
 {
   String tStr(weightStr);
   tStr.upcase();
+  String msg;
   if (tStr.contains(String("NONE"))) {
      wtType = NONE;
+     msg = String("Weighting type selected : None");
   } else if (tStr.contains(String("VAR"))) {
      wtType = VAR;
+     msg = String("Weighting type selected : Variance");
+  } else if (tStr.contains(String("TINT"))) {
+     wtType = TINT;
+     msg = String("Weighting type selected : Tnt");
   } else if (tStr.contains(String("TSYS"))) {
      wtType = TSYS;
-     throw(AipsError("T_sys weighting not yet implemented"));
+     msg = String("Weighting type selected : Tsys");
   } else {
-    throw(AipsError("Unrecognized weighting type"));
-  }
+     msg = String("Weighting type selected : None");
+     throw(AipsError("Unrecognized weighting type"));
+  }
+//
+  if (listType) cout << msg << endl;
 }
 

trunk/src/SDMath.h

@@ -136 +136 @@
 // Weighting type for time averaging
 
-  enum WeightType {NONE,VAR,TSYS};
+  enum WeightType {NONE,VAR,TSYS,TINT};
 
 // Function to use accumulate data during time averaging
@@ -143 +143 @@
                    casa::MaskedArray<casa::Float>& sum, casa::Array<casa::Float>& sumSq,
                    casa::Array<casa::Float>& nPts, casa::Array<casa::Float>& tSysSum,
+                   casa::Array<casa::Float>& tSysSqSum,
                    const casa::Array<casa::Float>& tSys,  const casa::Array<casa::Float>& nInc,
                    const casa::Vector<casa::Bool>& mask, casa::Double time, casa::Double interval,
@@ -155 +156 @@
 // Convert weight string to enum value
 
-   void convertWeightString (WeightType& wt, const casa::String& weightStr) const;
+   void convertWeightString (WeightType& wt, const casa::String& weightStr,
+                             casa::Bool listType) const;
 
 // Convert interpolation type string
@@ -233 +235 @@
   void normalize (casa::MaskedArray<casa::Float>& data,
                   const casa::Array<casa::Float>& sumSq,
+                  const casa::Array<casa::Float>& tSysSumSq,
                   const casa::Array<casa::Float>& nPts,
+                  casa::Double intSum,
                   WeightType wtType, casa::Int axis, casa::Int nAxes) const;