Changeset 1066

Timestamp:

07/04/06 11:22:42 (18 years ago)

Author:

mar637

Message:

Enhancement Ticket #50; arbitrary quotients.

Location:

trunk/src

Files:

• STMath.cpp (modified) (6 diffs)
• STMath.h (modified) (2 diffs)
• STMathWrapper.h (modified) (1 diff)
• python_STMath.cpp (modified) (1 diff)

trunk/src/STMath.cpp

@@ -21 +21 @@
 #include <casa/Arrays/ArrayLogical.h>
 #include <casa/Arrays/ArrayMath.h>
+#include <casa/Arrays/Slice.h>
+#include <casa/Arrays/Slicer.h>
 #include <casa/Containers/RecordField.h>
 #include <tables/Tables/TableRow.h>
@@ -66 +68 @@
 {
   if ( avmode == "SCAN" && in.size() != 1 )
-    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
+    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
+                    "Use merge first."));
   WeightType wtype = stringToWeight(weight);
 
@@ -186 +189 @@
 }
 
+
 CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
                                              bool droprows)
@@ -248 +252 @@
 }
 
-CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in,
-                                          const std::string & mode,
-                                          bool preserve )
+CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
+                                              const std::string & mode,
+                                              bool preserve )
 {
   /// @todo make other modes available
@@ -320 +324 @@
   return out;
 }
+
+
+CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
+                                          const CountedPtr< Scantable > & off,
+                                          bool preserve )
+{
+  bool insitu = insitu_;
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(on, false);
+  setInsitu(insitu);
+  Table& tout = out->table();
+  const Table& toff = off->table();
+  TableIterator sit(tout, "SCANNO");
+  TableIterator s2it(toff, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table ton = sit.table();
+    TableRow row(ton);
+    Table t = s2it.table();
+    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
+    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
+    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
+    for (uInt i=0; i < ton.nrow(); ++i) {
+      const TableRecord& rec = row.get(i);
+      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
+                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
+                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
+      TableRow offrow(offsel);
+      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
+      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
+      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
+      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
+      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
+      Vector<Float> specon, tsyson;
+      outtsysCol.get(i, tsyson);
+      outspecCol.get(i, specon);
+      Vector<uChar> flagon;
+      outflagCol.get(i, flagon);
+      MaskedArray<Float> mon = maskedArray(specon, flagon);
+      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
+      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
+      if (preserve) {
+        quot -= tsysoffscalar;
+      } else {
+        quot -= tsyson[0];
+      }
+      outspecCol.put(i, quot.getArray());
+      outflagCol.put(i, flagsFromMA(quot));
+    }
+    ++sit;
+    ++s2it;
+    // take the first off for each on scan which doesn't have a
+    // matching off scan
+    // non <= noff:  matching pairs, non > noff matching pairs then first off
+    if ( s2it.pastEnd() ) s2it.reset();
+  }
+  return out;
+}
+
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
@@ -1196 +1258 @@
   return out;
 }
+

trunk/src/STMath.h

@@ -58 +58 @@
              const std::string& weight = "NONE",
              const std::string& avmode = "SCAN");
+
   casa::CountedPtr< Scantable >
     averagePolarisations( const casa::CountedPtr< Scantable > & in,
@@ -67 +68 @@
                   const std::string& mode, bool tsys=false );
 
-  casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& in,
-                                        const std::string& mode = "NEAREST",
+  casa::CountedPtr<Scantable> autoQuotient(const casa::CountedPtr<Scantable>& in,
+                                           const std::string& mode = "NEAREST",
+                                           bool preserve = true);
+
+  casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& on,
+                                        const casa::CountedPtr<Scantable>& off,
                                         bool preserve = true );
 

trunk/src/STMathWrapper.h

@@ -57 +57 @@
   { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); }
 
-  ScantableWrapper quotient( const ScantableWrapper& in,
-                             const std::string& mode = "NEAREST",
+  ScantableWrapper autoQuotient( const ScantableWrapper& in,
+                                 const std::string& mode = "NEAREST",
+                                 bool preserve = true )
+  { return ScantableWrapper(STMath::autoQuotient(in.getCP(), mode, preserve)); }
+
+  ScantableWrapper quotient( const ScantableWrapper& on,
+                             const ScantableWrapper& off,
                              bool preserve = true )
-  { return ScantableWrapper(STMath::quotient(in.getCP(), mode, preserve)); }
+  { return ScantableWrapper( STMath::quotient( on.getCP(), off.getCP(),
+                                               preserve ) ); }
 
   ScantableWrapper

trunk/src/python_STMath.cpp

@@ -47 +47 @@
         .def("_averagepol", &STMathWrapper::averagePolarisations)
         .def("_unaryop", &STMathWrapper::unaryOperate)
+        .def("_auto_quotient", &STMathWrapper::autoQuotient)
         .def("_quotient", &STMathWrapper::quotient)
         .def("_stats", &STMathWrapper::statistic)