Changeset 294 for trunk

Timestamp:

01/25/05 14:56:37 (20 years ago)

Author:

kil064

Message:

move some code into function generateVelocityAligners
add arg. doTSys to binaryOperate and unaryOperate

Location:

trunk/src

Files:

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

trunk/src/SDMath.cc

@@ -114 +114 @@
 
    std::vector<std::string> info = in.getCoordInfo();
+
+// Parse unit ("" means channels)
+
    String velUnit(info[0]);
    if (velUnit.length()==0) {
@@ -123 +126 @@
       }
    }
-//
+
+// Parse doppler
+
    String dopplerStr(info[2]);
    String velSystemStr(info[1]);
@@ -130 +135 @@
       throw(AipsError("You have not set a velocity frame different from the initial - use function set_freqframe"));
    }
-//
+
+// Parse frequency system
+
    MFrequency::Types velSystem;
    MFrequency::getType(velSystem, velSystemStr);
@@ -377 +384 @@
 CountedPtr<SDMemTable> SDMath::binaryOperate (const CountedPtr<SDMemTable>& left,
                                               const CountedPtr<SDMemTable>& right,
-                                              const String& op, Bool preserve)  const
+                                              const String& op, Bool preserve,
+                                              Bool doTSys)  const
 {
 
@@ -395 +403 @@
   } else if (op2=="QUOTIENT") {
      what = 4;
+     doTSys = True;
   } else {
     throw( AipsError("Unrecognized operation"));
@@ -416 +425 @@
 // TSys columns
 
-  ROArrayColumn<Float> tSysLeft(tLeft, "TSYS");
-  ROArrayColumn<Float> tSysRight(tRight, "TSYS");
+  ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
+  if (doTSys) {
+     tSysLeftCol.attach(tLeft, "TSYS");
+     tSysRightCol.attach(tRight, "TSYS");
+  }
 
 // First row for right
 
   Array<Float> tSysLeftArr, tSysRightArr;
-  tSysRight.get(0, tSysRightArr);
+  if (doTSys) tSysRightCol.get(0, tSysRightArr);
   MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
   IPosition shpRight = pMRight->shape();
@@ -438 +450 @@
      MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
      IPosition shpLeft = mLeft.shape();
-     tSysLeft.get(i, tSysLeftArr);
+     if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
 //
      if (nRowRight>1) {
@@ -444 +456 @@
         pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
         shpRight = pMRight->shape();
-        tSysRight.get(i, tSysRightArr);
+        if (doTSys) tSysRightCol.get(i, tSysRightArr);
      }
 //
@@ -450 +462 @@
         throw(AipsError("left and right scan tables are not conformant"));
      }
-     if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
-        throw(AipsError("left and right Tsys data are not conformant"));
-     }
-     if (!shpRight.isEqual(tSysRightArr.shape())) {
-        throw(AipsError("left and right scan tables are not conformant"));
+     if (doTSys) {
+        if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
+           throw(AipsError("left and right Tsys data are not conformant"));
+        }
+        if (!shpRight.isEqual(tSysRightArr.shape())) {
+           throw(AipsError("left and right scan tables are not conformant"));
+        }
      }
 
@@ -466 +480 @@
         MaskedArray<Float> tmp = mLeft + *pMRight;
         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        sc.putTsys(tSysLeftArr+tSysRightArr);
+        if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr);
      } else if (what==1) {
         MaskedArray<Float> tmp = mLeft - *pMRight;
         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        sc.putTsys(tSysLeftArr-tSysRightArr);
+        if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr);
      } else if (what==2) {
         MaskedArray<Float> tmp = mLeft * *pMRight;
         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        sc.putTsys(tSysLeftArr*tSysRightArr);
+        if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr);
      } else if (what==3) {
         MaskedArray<Float> tmp = mLeft / *pMRight;
         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        sc.putTsys(tSysLeftArr/tSysRightArr);
+        if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr);
      } else if (what==4) {
         if (preserve) {     
@@ -621 +635 @@
 
 SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
-                                 uInt what) const
+                                 uInt what, Bool doTSys) const
 //
 // what = 0   Multiply
@@ -628 +642 @@
    SDMemTable* pOut = new SDMemTable(in,False);
    const Table& tOut = pOut->table();
-   ArrayColumn<Float> spec(tOut,"SPECTRA");  
+   ArrayColumn<Float> specCol(tOut,"SPECTRA");  
+   ArrayColumn<Float> tSysCol(tOut,"TSYS");  
+   Array<Float> tSysArr;
 //
    if (doAll) {
       for (uInt i=0; i < tOut.nrow(); i++) {
+
+// Modify data
+
          MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
-//
          if (what==0) {
             dataIn  *= val;
@@ -639 +657 @@
             dataIn += val;
          }
-//
-         spec.put(i, dataIn.getArray());
+         specCol.put(i, dataIn.getArray());
+
+// Modify Tsys
+
+         if (doTSys) {
+            tSysCol.get(i, tSysArr);
+            if (what==0) {
+               tSysArr *= val;
+            } else if (what==1) {
+               tSysArr += val;
+            }
+            tSysCol.put(i, tSysArr);
+         }
       }
    } else {
@@ -650 +679 @@
 //
       for (uInt i=0; i < tOut.nrow(); i++) {
+
+// Modify data
+
          MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
          MaskedArray<Float> dataIn2 = dataIn(start,end);    // Reference
-//
          if (what==0) {
             dataIn2 *= val;
@@ -658 +689 @@
             dataIn2 += val;
          }
-//
-         spec.put(i, dataIn.getArray());
+         specCol.put(i, dataIn.getArray());
+
+// Modify Tsys
+
+         if (doTSys) {
+            tSysCol.get(i, tSysArr);
+            Array<Float> tSysArr2 = tSysArr(start,end);     // Reference
+            if (what==0) {
+               tSysArr2 *= val;
+            } else if (what==1) {
+               tSysArr2 += val;
+            }
+            tSysCol.put(i, tSysArr);
+         }
       }
    }
@@ -1118 +1161 @@
 // Get Columns from Table
 
-  ROScalarColumn<Double> mjdCol(tabIn, "TIME");
-  ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
-  ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
-//
-  Vector<Double> times = mjdCol.getColumn();
-  Vector<String> srcNames = srcCol.getColumn();
-  Vector<uInt> freqID;
+   ROScalarColumn<Double> mjdCol(tabIn, "TIME");
+   ROScalarColumn<String> srcCol(tabIn, "SRCNAME");
+   ROArrayColumn<uInt> fqIDCol(tabIn, "FREQID");
+//
+   Vector<Double> times = mjdCol.getColumn();
+   Vector<String> srcNames = srcCol.getColumn();
+   Vector<uInt> freqID;
 
 // Generate Source table
@@ -1134 +1177 @@
    cerr << "Found " << srcTab.nelements() << " sources to align " << endl;
 
-// Set reference Epoch to time of first row or given String
+// Get reference Epoch to time of first row or given String
 
    Unit DAY(String("d"));
@@ -1146 +1189 @@
    cerr << "Aligning at reference Epoch " << formatEpoch(refEpoch) << endl;
 
-// Set Reference Position
+// Get Reference Position
 
    MPosition refPos = in.getAntennaPosition();
@@ -1159 +1202 @@
 
    PtrBlock<VelocityAligner<Float>* > vA(nFreqIDs*nSrcTab);
-   for (uInt fqID=0; fqID<nFreqIDs; fqID++) {
-      SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
-      for (uInt iSrc=0; iSrc<nSrcTab; iSrc++) {
-         MDirection refDir = in.getDirection(firstRow[iSrc]);
-         uInt idx = (iSrc*nFreqIDs) + fqID;
-         vA[idx] = new VelocityAligner<Float>(sC, nChan, refEpoch, refDir, refPos,
-                                              velUnit, doppler, velSystem);
+   generateVelocityAligners (vA, in, nChan, nFreqIDs, nSrcTab, firstRow,
+                             velSystem, velUnit, doppler,  refPos, refEpoch);
+
+// New output Table
+
+   SDMemTable* pTabOut = new SDMemTable(in,True);
+
+// Loop over rows in Table
+
+   const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
+   VelocityAligner<Float>::Method method = VelocityAligner<Float>::LINEAR;
+   Bool extrapolate=False;
+   Bool useCachedAbcissa = False;
+   Bool first = True;
+   Bool ok;
+   Vector<Float> yOut;
+   Vector<Bool> maskOut;
+   uInt ifIdx, vaIdx;
+//
+   for (uInt iRow=0; iRow<nRows; ++iRow) {
+      if (iRow%10==0) {
+         cerr << "Processing row " << iRow << endl;
       }
-   }
-
-// New output Table
-
-   SDMemTable* pTabOut = new SDMemTable(in,True);
-
-// Loop over rows in Table
-
-  const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
-  VelocityAligner<Float>::Method method = VelocityAligner<Float>::LINEAR;
-  Bool extrapolate=False;
-  Bool useCachedAbcissa = False;
-  Bool first = True;
-  Bool ok;
-  Vector<Float> yOut;
-  Vector<Bool> maskOut;
-  uInt ifIdx, vaIdx;
-//
-  for (uInt iRow=0; iRow<nRows; ++iRow) {
-     if (iRow%10==0) {
-        cerr << "Processing row " << iRow << endl;
-     }
 
 // Get EPoch
 
-    Quantum<Double> tQ2(times[iRow],DAY);
-    MVEpoch mv2(tQ2);
-    MEpoch epoch(mv2, epochRef);
+     Quantum<Double> tQ2(times[iRow],DAY);
+     MVEpoch mv2(tQ2);
+     MEpoch epoch(mv2, epochRef);
 
 // Get FreqID vector.  One freqID per IF
 
-    fqIDCol.get(iRow, freqID);
+     fqIDCol.get(iRow, freqID);
 
 // Get copy of data
     
-    const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
-    Array<Float> values = mArrIn.getArray();
-    Array<Bool> mask = mArrIn.getMask(); 
+     const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
+     Array<Float> values = mArrIn.getArray();
+     Array<Bool> mask = mArrIn.getMask(); 
 
 // cerr << "values in = " << values(IPosition(4,0,0,0,0),IPosition(4,0,0,0,9)) << endl;
@@ -1216 +1252 @@
 // all beams get same position and velocoity abcissa.
 
-    ArrayIterator<Float> itValuesPlane(values, polChanAxes);
-    ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
-    while (!itValuesPlane.pastEnd()) {
+     ArrayIterator<Float> itValuesPlane(values, polChanAxes);
+     ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
+     while (!itValuesPlane.pastEnd()) {
 
 // Find the IF index and then the VA PtrBlock index
 
-       const IPosition& pos = itValuesPlane.pos();
-       ifIdx = pos(asap::IFAxis);
-       vaIdx = (srcIdx[iRow]*nFreqIDs) + freqID[ifIdx];
-//
-       VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
-       VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
-//
-       first = True;
-       useCachedAbcissa=False;
-       while (!itValuesVec.pastEnd()) {     
-          ok = vA[vaIdx]->align (yOut, maskOut, itValuesVec.vector(),
-                                 itMaskVec.vector(), epoch, useCachedAbcissa,
-                                 method, extrapolate); 
-          itValuesVec.vector() = yOut;
-          itMaskVec.vector() = maskOut;
-//
-          itValuesVec.next();
-          itMaskVec.next();
-//
-          if (first) {
-             useCachedAbcissa = True;
-             first = False;
-          }
-       }
+        const IPosition& pos = itValuesPlane.pos();
+        ifIdx = pos(asap::IFAxis);
+        vaIdx = (srcIdx[iRow]*nFreqIDs) + freqID[ifIdx];
+//
+        VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1); 
+        VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
+//
+        first = True;
+        useCachedAbcissa=False;
+        while (!itValuesVec.pastEnd()) {     
+           ok = vA[vaIdx]->align (yOut, maskOut, itValuesVec.vector(),
+                                  itMaskVec.vector(), epoch, useCachedAbcissa,
+                                  method, extrapolate); 
+           itValuesVec.vector() = yOut;
+           itMaskVec.vector() = maskOut;
+//
+           itValuesVec.next();
+           itMaskVec.next();
+//
+           if (first) {
+              useCachedAbcissa = True;
+              first = False;
+           }
+        }
 //
        itValuesPlane.next();
        itMaskPlane.next();
-    }
+     }
 
 // cerr << "values out = " << values(IPosition(4,0,0,0,0),IPosition(4,0,0,0,9)) << endl;
@@ -1259 +1295 @@
 //
     pTabOut->putSDContainer(sc);
-  }
+   }
 
 // Clean up PointerBlock
@@ -1628 +1664 @@
 }
 
+
+void SDMath::generateVelocityAligners (PtrBlock<VelocityAligner<Float>* >& vA, 
+                                       const SDMemTable& in, uInt nChan,
+                                       uInt nFreqIDs, uInt nSrcTab, 
+                                       const Vector<uInt>& firstRow,
+                                       MFrequency::Types velSystem,
+                                       const String& velUnit,
+                                       MDoppler::Types doppler,
+                                       const MPosition& refPos,
+                                       const MEpoch& refEpoch) const
+{
+   for (uInt fqID=0; fqID<nFreqIDs; fqID++) {
+      SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
+      for (uInt iSrc=0; iSrc<nSrcTab; iSrc++) {
+         MDirection refDir = in.getDirection(firstRow[iSrc]);
+         uInt idx = (iSrc*nFreqIDs) + fqID;
+         vA[idx] = new VelocityAligner<Float>(sC, nChan, refEpoch, refDir, refPos,
+                                              velUnit, doppler, velSystem);
+      }
+   }
+}
+

trunk/src/SDMath.h

@@ -36 +36 @@
 #include <casa/aips.h>
 #include <casa/Utilities/CountedPtr.h>
+#include <coordinates/Coordinates/VelocityAligner.h>
+
 #include "SDDefs.h"
 
 class casa::Table;
 class casa::MEpoch;
+class casa::MPosition;
+template<class T> class casa::PtrBlock;
+//template<class T> class casa::VelocityAligner;
 
 
@@ -65 +70 @@
    casa::CountedPtr<SDMemTable> binaryOperate (const casa::CountedPtr<SDMemTable>& left,
                                                const casa::CountedPtr<SDMemTable>& right,
-                                               const casa::String& op, casa::Bool preserve) const;
+                                               const casa::String& op, casa::Bool preserve,
+                                               casa::Bool tSys) const;
 
 // Average in time
@@ -104 +110 @@
 // Simple unary mathematical operations.  what=0 (mul) or 1 (add)
    SDMemTable* unaryOperate(const SDMemTable& in, casa::Float offset, 
-                            casa::Bool doAll, casa::uInt what) const;
+                            casa::Bool doAll, casa::uInt what, casa::Bool tSys) const;
 
 // Average polarizations
@@ -184 +190 @@
                              const casa::Vector<casa::String>& srcNames) const;
 
+// Generate velocity aligners
+   void generateVelocityAligners (casa::PtrBlock<casa::VelocityAligner<casa::Float>* >& vA,
+                                  const SDMemTable& in, casa::uInt nChan,
+                                  casa::uInt nFreqIDs, casa::uInt nSrcTab,
+                                  const casa::Vector<casa::uInt>& firstRow,
+                                  casa::MFrequency::Types velSystem,
+                                  const casa::String& velUnit,
+                                  casa::MDoppler::Types doppler,
+                                  const casa::MPosition& refPos,
+                                  const casa::MEpoch& refEpoch) const;
+
 // Align in Velocity
    SDMemTable* velocityAlign (const SDMemTable& in,