Changeset 262

Timestamp:

01/22/05 17:37:42 (20 years ago)

Author:

kil064

Message:

add functiom VelocityAlignment

Location:

trunk/src

Files:

• SDMath.cc (modified) (12 diffs)
• SDMath.h (modified) (2 diffs)
• SDMathWrapper.cc (modified) (1 diff)
• SDMathWrapper.h (modified) (1 diff)

trunk/src/SDMath.cc

@@ -44 +44 @@
 #include <casa/BasicMath/Math.h>
 #include <casa/Containers/Block.h>
+#include <casa/Exceptions.h>
+#include <casa/Quanta/Quantum.h>
+#include <casa/Quanta/Unit.h>
+#include <casa/Quanta/MVEpoch.h>
 #include <casa/Quanta/QC.h>
 #include <casa/Utilities/Assert.h>
-#include <casa/Exceptions.h>
+
+#include <coordinates/Coordinates/SpectralCoordinate.h>
+#include <coordinates/Coordinates/CoordinateSystem.h>
+#include <coordinates/Coordinates/CoordinateUtil.h>
+#include <coordinates/Coordinates/VelocityAligner.h>
+
+#include <lattices/Lattices/LatticeUtilities.h>
+#include <lattices/Lattices/RebinLattice.h>
+
+#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MDirection.h>
+#include <measures/Measures/MPosition.h>
 
 #include <scimath/Mathematics/VectorKernel.h>
@@ -58 +73 @@
 #include <tables/Tables/ReadAsciiTable.h>
 
-#include <lattices/Lattices/LatticeUtilities.h>
-#include <lattices/Lattices/RebinLattice.h>
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-#include <coordinates/Coordinates/VelocityAligner.h>
-
 #include "MathUtils.h"
 #include "SDDefs.h"
@@ -96 +104 @@
 SDMath::~SDMath()
 {;}
+
+
+
+SDMemTable* SDMath::velocityAlignment (const SDMemTable& in) const
+{
+
+// Get velocity/frame info
+
+   std::vector<std::string> info = in.getCoordInfo();
+   String velUnit(info[0]);
+   if (velUnit.length()==0) {
+      throw(AipsError("You have not set a velocity abcissa unit - use function set_unit"));
+   } else {
+      Unit velUnitU(velUnit);
+      if (velUnitU!=Unit(String("m/s"))) {
+         throw(AipsError("Specified abcissa unit is not consistent with km/s - use function set_unit"));
+      }
+   }
+//
+   String dopplerStr(info[2]);
+   String velSystemStr(info[1]);
+   String velBaseSystemStr(info[3]);
+   if (velBaseSystemStr==velSystemStr) {
+      throw(AipsError("You have not set a velocity frame different from the initial - use function set_frame"));
+   }
+
+cerr << "unit, frame, doppler, baseframe = " << velUnit << ", " << velSystemStr << ", " << 
+               dopplerStr << ", " << velBaseSystemStr << endl;
+//
+   MFrequency::Types velSystem;
+   MFrequency::getType(velSystem, velSystemStr);
+   MDoppler::Types doppler;
+   MDoppler::getType(doppler, dopplerStr);
+
+// Get Header
+
+   SDHeader sh = in.getSDHeader();
+   const uInt nChan = sh.nchan;
+   const uInt nRows = in.nRow();
+
+// Get Table reference
+
+   const Table& tabIn = in.table();
+
+// 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;
+
+// Generate Source table
+
+   Vector<String> srcTab;
+   Vector<uInt> srcIdx, firstRow;
+   generateSourceTable (srcTab, srcIdx, firstRow, srcNames);
+   const uInt nSrcTab = srcTab.nelements();
+/*
+   cerr << "source table = " << srcTab << endl;
+   cerr << "source idx = " << srcIdx << endl;
+   cerr << "first row = " << firstRow << endl;
+*/
+
+// Set reference Epoch to time of first row
+
+   MEpoch::Ref timeRef = MEpoch::Ref(in.getTimeReference());
+   Unit DAY(String("d"));
+   Quantum<Double> tQ(times[0], DAY);
+   MVEpoch mve(tQ);
+   MEpoch refTime(mve, timeRef);
+
+// Set Reference Position
+
+   Vector<Double> antPos = sh.antennaposition;
+   MVPosition mvpos(antPos[0], antPos[1], antPos[2]);
+   MPosition refPos(mvpos);
+
+// Get Frequency Table
+
+   SDFrequencyTable fTab = in.getSDFreqTable();
+   const uInt nFreqIDs = fTab.length();
+
+// Create VelocityAligner Block. One VA for each possible
+// source/freqID combination
+
+   PtrBlock<VelocityAligner<Float>* > vA(nFreqIDs*nSrcTab);
+   for (uInt fqID=0; fqID<nFreqIDs; fqID++) {
+      SpectralCoordinate sC = in.getCoordinate(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, refTime, refDir, refPos,
+                                              velUnit, doppler, velSystem);
+      }
+   }
+
+// 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;
+  Vector<Float> yOut;
+  Vector<Bool> maskOut;
+//
+  for (uInt iRow=0; iRow<nRows; ++iRow) {
+
+// Get EPoch
+
+   Quantum<Double> tQ2(times[iRow],DAY);
+   MVEpoch mv2(tQ2);
+   MEpoch epoch(mv2, timeRef);
+
+// Get FreqID vector.  One freqID per IF
+
+    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(); 
+
+// cerr << "values in = " << values(IPosition(4,0,0,0,0),IPosition(4,0,0,0,9)) << endl;
+
+// For each row, the Velocity abcissa will be the same regardless
+// of polarization.  For all other axes (IF and BEAM) the abcissa
+// will change.  So we iterate through the data by pol-chan planes
+// to mimimize the work.  At this point, I think the Direction
+// is stored as the same for each beam. DOn't know where the 
+// offsets are or what to do about them right now.  For now
+// all beams get same position and velocoity abcissa.
+
+    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();
+       uInt ifIdx = pos(asap::IFAxis);
+       uInt vaIdx = (srcIdx[iRow]*nFreqIDs) + freqID[ifIdx];
+//cerr << "VA idx = " << vaIdx << endl;
+//
+       VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
+       VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
+//
+       first = True;
+       useCachedAbcissa=False;
+       while (!itValuesVec.pastEnd()) {     
+          Bool 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;
+
+// Create and put back
+
+    SDContainer sc = in.getSDContainer(iRow);
+    putDataInSDC(sc, values, mask);
+//
+    pTabOut->putSDContainer(sc);
+  }
+
+// Clean up PointerBlock
+
+  for (uInt i=0; i<vA.nelements(); i++) delete vA[i];
+//
+  return pTabOut;
+}
 
 
@@ -197 +398 @@
   String sourceName, oldSourceName, sourceNameStart;
   Vector<uInt> freqID, freqIDStart, oldFreqID;
-
-// Velocity Aligner. We need an aligner for each Direction and FreqID
-// combination.  I don't think there is anyway to know how many
-// directions there are.
-// For now, assume all Tables have the same Frequency Table
-
-/*
-  {
-     MEpoch::Ref timeRef(MEpoch::UTC);              // Should be in header
-     MDirection::Types dirRef(MDirection::J2000);   // Should be in header
-//
-     SDHeader sh = in[0].getSDHeader();
-     const uInt nChan = sh.nchan;
-//
-     const SDFrequencyTable freqTab = in[0]->getSDFreqTable();
-     const uInt nFreqID = freqTab.length();
-     PtrBlock<const VelocityAligner<Float>* > vA(nFreqID);
-
-// Get first time from first table
-
-     const Table& tabIn0 = in[0]->table();
-     mjdCol.attach(tabIn0, "TIME");
-     Double dTmp;
-     mjdCol.get(0, dTmp);
-     MVEpoch tmp2(Quantum<Double>(dTmp, Unit(String("d"))));
-     MEpoch epoch(tmp2, timeRef);
-//
-     for (uInt freqID=0; freqID<nFreqID; freqID++) {
-        SpectralCoordinate sC = in[0]->getCoordinate(freqID);
-        vA[freqID] = new VelocityAligner<Float>(sC, nChan, epoch, const MDirection& dir, 
-                                                const MPosition& pos, const String& velUnit, 
-                                                MDoppler::Types velType, MFrequency::Types velFreqSystem)
-     }
-  }
-*/
 
 // Loop over tables
@@ -686 +852 @@
 
    const uInt nRows = in.nRow();
-   const uInt polAxis = asap::PolAxis;                     // Polarization axis
-   const uInt chanAxis = asap::ChanAxis;                    // Spectrum axis
 
 // Create output Table and reshape number of polarizations
@@ -701 +865 @@
   const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
   IPosition shapeOut(shapeIn);
-  shapeOut(polAxis) = 1;                          // Average all polarizations
-//
-  const uInt nChan = shapeIn(chanAxis);
+  shapeOut(asap::PolAxis) = 1;                          // Average all polarizations
+//
+  const uInt nChan = shapeIn(asap::ChanAxis);
   const IPosition vecShapeOut(4,1,1,1,nChan);     // A multi-dim form of a Vector shape
   IPosition start(4), end(4);
@@ -711 +875 @@
   Array<Float> outData(shapeOut, 0.0);
   Array<Bool> outMask(shapeOut, True);
-  const IPosition axes(2, 2, 3);              // pol-channel plane
+  const IPosition axes(2, asap::PolAxis, asap::ChanAxis);              // pol-channel plane
 // 
-  const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
+  const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
 
 // Loop over rows
@@ -787 +951 @@
         start = pos;
         end = pos; 
-        end(chanAxis) = nChan-1;
+        end(asap::ChanAxis) = nChan-1;
         outData(start,end) = vecSum.getArray().reform(vecShapeOut);
         outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
@@ -854 +1018 @@
     
     const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(ri));
-    AlwaysAssert(dataIn.shape()(chanAxis)==nChan, AipsError);
+    AlwaysAssert(dataIn.shape()(asap::ChanAxis)==nChan, AipsError);
 //
     Array<Float> valuesIn = dataIn.getArray();
@@ -886 +1050 @@
 // Set multi-dim Vector shape
 
-       shapeOut(chanAxis) = valuesIn.shape()(chanAxis);
+       shapeOut(asap::ChanAxis) = valuesIn.shape()(chanAxis);
 
 // Stuff about with shapes so that we don't have conformance run-time errors
@@ -916 +1080 @@
   return pTabOut;
 }
+
 
 
@@ -1438 +1603 @@
 
 
+void SDMath::generateSourceTable (Vector<String>& srcTab,
+                                  Vector<uInt>& srcIdx,
+                                  Vector<uInt>& firstRow,
+                                  const Vector<String>& srcNames) const
+//
+// This algorithm assumes that if there are multiple beams
+// that the source names are diffent.  Oterwise we would need
+// to look atthe direction for each beam...
+//
+{
+   const uInt nRow = srcNames.nelements();
+   srcTab.resize(0);
+   srcIdx.resize(nRow);
+   firstRow.resize(0);
+//
+   uInt nSrc = 0;
+   for (uInt i=0; i<nRow; i++) {
+      String srcName = srcNames[i];
+      
+// Do we have this source already ?
+
+      Int idx = -1;
+      if (nSrc>0) {
+         for (uInt j=0; j<nSrc; j++) {
+           if (srcName==srcTab[j]) {
+              idx = j;
+              break;
+           }
+         }
+      }
+
+// Add new entry if not found
+
+      if (idx==-1) {
+         nSrc++;
+         srcTab.resize(nSrc,True);
+         srcTab(nSrc-1) = srcName;
+         idx = nSrc-1;
+//
+         firstRow.resize(nSrc,True);
+         firstRow(nSrc-1) = i;       // First row for which this source occurs
+      }
+
+// Set index for this row
+
+      srcIdx[i] = idx;
+   }
+}

trunk/src/SDMath.h

@@ -94 +94 @@
                               const casa::String& method, casa::Bool doAll) const;
 
+// Velocity Alignment
+   SDMemTable* velocityAlignment (const SDMemTable& in) const;
+
 // Opacity correction
    SDMemTable* opacity (const SDMemTable& in, casa::Float tau, casa::Bool doAll) const;
@@ -169 +172 @@
                            casa::Bool doAll, const casa::Vector<casa::Float>& factor) const;
 
-
-
 // Read ascii file into a Table
 
    casa::Table readAsciiFile (const casa::String& fileName) const;
+
+// Generate source table
+   void generateSourceTable (casa::Vector<casa::String>& srcTab,
+                             casa::Vector<casa::uInt>& srcIdx,
+                             casa::Vector<casa::uInt>& firstRow,
+                             const casa::Vector<casa::String>& srcNames) const;
+
 };
 

trunk/src/SDMathWrapper.cc

@@ -228 +228 @@
 }
 
+SDMemTableWrapper SDMathWrapper::velocityAlignment (const SDMemTableWrapper& in)
+{
+  const CountedPtr<SDMemTable>& pIn = in.getCP();
+  SDMath sdm;
+  return CountedPtr<SDMemTable>(sdm.velocityAlignment(*pIn));
+}
+
 void SDMathWrapper::opacityInSitu(SDMemTableWrapper& in, float tau, bool doAll)
 {

trunk/src/SDMathWrapper.h

@@ -84 +84 @@
                                   const std::string& method, bool doAll);
 
+// Apply velocity alignment
+  SDMemTableWrapper velocityAlignment(const SDMemTableWrapper& in);
+
 // Apply opacity correction
   void opacityInSitu (SDMemTableWrapper& in, float tau, bool doAll);