Changeset 418 for trunk/src

Timestamp:

02/14/05 12:35:18 (20 years ago)

Author:

kil064

Message:

First go at handling polarimetric data

• Add STOKES (virtual) column
• Add function getStokesSpectrum
• Add function rotateXYPhase

Location:

trunk/src

Files:

• SDMemTable.cc (modified) (14 diffs)
• SDMemTable.h (modified) (5 diffs)

trunk/src/SDMemTable.cc

@@ -41 +41 @@
 #include <casa/Arrays/ArrayAccessor.h>
 #include <casa/Arrays/Vector.h>
+#include <casa/BasicMath/Math.h>
 #include <casa/Quanta/MVAngle.h>
 
@@ -61 +62 @@
 #include "SDContainer.h"
 #include "MathUtils.h"
+#include "SDPol.h"
 
 
@@ -164 +166 @@
   td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
   td.addColumn(ArrayColumnDesc<Float>("TSYS"));
+  td.addColumn(ArrayColumnDesc<Float>("STOKES"));
   td.addColumn(ScalarColumnDesc<Int>("SCANID"));
   td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
@@ -178 +181 @@
   td.addColumn(ArrayColumnDesc<String>("HISTORY"));
 
-  // Now create a new table from the description.
+  // Now create Table SetUp from the description.
 
   SetupNewTable aNewTab("dummy", td, Table::New);
+
+// Bind the Stokes Virtual machine to the STOKES column
+// Because we don't know how many polarizations will
+// be in the data at this point, we must bind the
+// Virtual Engine regardless.  The STOKES column won't
+// be accessed if not appropriate (nPol=4)
+
+   SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));
+   aNewTab.bindColumn ("STOKES", stokesEngine);
+
+// Create Table
+
   table_ = Table(aNewTab, Table::Memory, 0);
 }
+
 
 void SDMemTable::attach()
@@ -191 +207 @@
   flagsCol_.attach(table_, "FLAGTRA");
   tsCol_.attach(table_, "TSYS");
+  stokesCol_.attach(table_, "STOKES");
   scanCol_.attach(table_, "SCANID");
   integrCol_.attach(table_, "INTERVAL");
@@ -318 +335 @@
   return mask;
 }
+
+
+
 std::vector<float> SDMemTable::getSpectrum(Int whichRow) const 
 {
-  std::vector<float> spectrum;
   Array<Float> arr;
   specCol_.get(whichRow, arr);
+//
+  return getFloatSpectrum (arr);
+}
+
+std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, Bool doPol) const 
+//
+// Gets
+//  doPol=False  : I,Q,U,V
+//  doPol=True   : I,P,PA,V   ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
+//
+{
+  if (nPol()!=4) {
+     throw (AipsError("You must have 4 polarizations to get the Stokes parameters"));
+  }
+  Array<Float> arr;
+  stokesCol_.get(whichRow, arr);
+//
+  if (doPol && (polSel_==1 || polSel_==2)) {
+     const IPosition& shape = arr.shape();
+     IPosition start(asap::nAxes,0);
+     IPosition end(shape-1);
+//
+     start(asap::PolAxis) = 1;                       // Q
+     end (asap::PolAxis) = 1;
+     Array<Float> Q = arr(start,end);
+//
+     start(asap::PolAxis) = 2;                       // U
+     end (asap::PolAxis) = 2;
+     Array<Float> U = arr(start,end);
+//
+     Array<Float> out;
+     if (polSel_==1) {                                        // P
+        out = SDPolUtil::polarizedIntensity(Q,U);
+     } else if (polSel_==2) {                                 // P.A.
+        out = SDPolUtil::positionAngle(Q,U);
+     }
+//
+     IPosition vecShape(1,shape(asap::ChanAxis));
+     Vector<Float> outV = out.reform(vecShape);
+     std::vector<float> spectrum(out.nelements());
+     for (uInt i=0; i<out.nelements(); i++) {
+        spectrum[i] = outV[i];
+     }
+     return spectrum;
+  } else {
+     return getFloatSpectrum (arr);
+  }
+}
+
+std::vector<float> SDMemTable::getFloatSpectrum (const Array<Float>& arr) const
+{
+
+// Iterate and extract
+
   ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
   aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
@@ -329 +402 @@
   ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
   aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
+//
+  std::vector<float> spectrum;
   for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
     spectrum.push_back(*i);
@@ -334 +409 @@
   return spectrum;
 }
+
+
 std::vector<string> SDMemTable::getCoordInfo() const
 {
@@ -537 +614 @@
   Array<Float> arr;
   specCol_.get(whichRow, arr);
+
+// Iterate and extract
+
   spectrum.resize(arr.shape()(3));
   ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
@@ -544 +624 @@
   ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
   aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
-
+//
   ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
   for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
@@ -551 +631 @@
   }
 }
+
+
 /*
 void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
@@ -591 +673 @@
   specCol_.get(whichRow, arr);
   flagsCol_.get(whichRow, farr);
-  Array<Bool> barr(farr.shape());convertArray(barr, farr);
+  Array<Bool> barr(farr.shape());
+  convertArray(barr, farr);
   MaskedArray<Float> marr;
   if (useSelection) {
@@ -634 +717 @@
   tsCol_.get(whichRow, arr);
   Float out;
+//
   IPosition ip(arr.shape());
-  ip(0) = beamSel_;ip(1) = IFSel_;ip(2) = polSel_;ip(3)=0;
+  ip(asap::BeamAxis) = beamSel_;
+  ip(asap::IFAxis) = IFSel_;
+  ip(asap::PolAxis) = polSel_;
+  ip(asap::ChanAxis)=0;               // First channel only
+//
   out = arr(ip);
   return out;
@@ -1476 +1564 @@
 }
 
+
+void SDMemTable::rotateXYPhase (Float value) 
+//
+// phase in degrees
+//
+{
+   if (nPol() != 4) {
+      throw(AipsError("You must have 4 polarizations to run this function"));
+   }
+//
+   IPosition start(asap::nAxes,0);
+   IPosition end(asap::nAxes);
+//
+   uInt nRow = specCol_.nrow();
+   Array<Float> data;
+   for (uInt i=0; i<nRow;++i) {
+      specCol_.get(i,data);
+      end = data.shape()-1;
+
+// Get polarization slice references
+
+      start(asap::PolAxis) = 2;
+      end(asap::PolAxis) = 2;
+      Array<Float> C3 = data(start,end);
+//
+      start(asap::PolAxis) = 3;
+      end(asap::PolAxis) = 3;
+      Array<Float> C4 = data(start,end);
+
+// Rotate
+
+      SDPolUtil::rotateXYPhase(C3, C4, value);
+
+// Put
+
+      specCol_.put(i,data);
+   }
+}

trunk/src/SDMemTable.h

@@ -94 +94 @@
   virtual std::vector<float> getSpectrum(casa::Int whichRow=0) const;
   virtual std::vector<bool> getMask(casa::Int whichRow=0) const;
+  virtual std::vector<float> getStokesSpectrum(casa::Int whichRow=0, 
+                                               casa::Bool doPol=casa::False) const;
 
   virtual casa::Float getTsys(casa::Int whichRow=0) const;
@@ -99 +101 @@
   virtual void getSpectrum(casa::Vector<casa::Float>& spectrum, 
                            casa::Int whichRow=0) const;
-
   //virtual void getMask(Vector<Bool>& mask,Int whichRow=0) const;
+
   std::vector<double> getRestFreqs() const;
   
@@ -218 +220 @@
   casa::MPosition getAntennaPosition() const;
 
+// Rotate phase of XY correlation by specified value (degrees)
+  void rotateXYPhase (casa::Float angle);
+
 // Helper function to check instrument (antenna) name and give enum
   static Instrument convertInstrument(const casa::String& instrument,
@@ -226 +231 @@
   casa::String formatSec(casa::Double x) const;
   casa::String formatDirection(const casa::MDirection& md) const;
+  std::vector<float> getFloatSpectrum (const casa::Array<casa::Float>& arr) const;
   void setup();
   void attach();
@@ -245 +251 @@
   casa::ArrayColumn<casa::uInt> freqidCol_, restfreqidCol_;
   casa::ArrayColumn<casa::String> histCol_;
+  casa::ArrayColumn<casa::Float> stokesCol_;
 };