Changeset 469

Timestamp:

02/18/05 11:11:01 (20 years ago)

Author:

kil064

Message:

new functions for the SDWriter to handle Stokes output
for SDFITS / MS. This is bit of a mess...

Location:

trunk/src

Files:

• SDPol.cc (modified) (2 diffs)
• SDPol.h (modified) (4 diffs)
• SDPol2.cc (modified) (2 diffs)

trunk/src/SDPol.cc

@@ -36 +36 @@
 #include <casa/Arrays/ArrayMath.h>
 #include <casa/Arrays/ArrayLogical.h>
+#include <casa/Arrays/VectorIter.h>
 #include <casa/Containers/Record.h>
 #include <casa/BasicSL/Constants.h>
@@ -359 +360 @@
    return stokes;
 }
+
+Array<casa::uChar> SDPolUtil::andArrays (const Array<casa::uChar>& in1,
+                                         const Array<casa::uChar>& in2)
+{
+   Array<uChar> out(in1.shape());
+//
+   Array<uChar>::const_iterator in1Iter;
+   Array<uChar>::const_iterator in2Iter;
+   Array<uChar>::iterator outIter;
+//
+   for (in1Iter=in1.begin(),in2Iter=in2.begin(),outIter=out.begin();
+        in1Iter!=in1.end(); ++in1Iter,++in2Iter,++outIter) {  
+      *outIter = *in1Iter & *in2Iter;
+   }
+   return out;
+}
+
+
+Array<Float> SDPolUtil::extractStokesForWriter (Array<Float>& in, const IPosition& start, const IPosition& end)
+//
+// start/end must already have applied the cursor selection of beam and IF
+// Extract specified Stokes for beam/IF and flip nChan and nPol for bloody SDwriter
+//
+{
+   IPosition shapeIn = in.shape();
+   uInt nChan = shapeIn(asap::ChanAxis);
+   uInt nPol = shapeIn(asap::PolAxis);
+//
+   IPosition shapeOut(2,nChan,nPol);
+   Array<Float> out(shapeOut);
+//
+   Array<Float> sliceRef = in(start,end);                        // Beam and IF now degenerate axes
+   ReadOnlyVectorIterator<Float> itIn(sliceRef, asap::ChanAxis);
+   VectorIterator<Float> itOut(out,0);
+   while (!itIn.pastEnd()) {
+      itOut.vector() = itIn.vector();
+// 
+      itIn.next();
+      itOut.next();
+   }
+//
+   return out;
+}
+

trunk/src/SDPol.h

@@ -155 +155 @@
                                                  const casa::Array<casa::Float>& U);
 
-// Rotate phase of Complex correlation C3+iC4 by phase (degrees)
+ // Rotate phase of Complex correlation C3+iC4 by phase (degrees)
   static void rotateXYPhase (casa::Array<casa::Float>& C3,
                              casa::Array<casa::Float>& C4,
@@ -163 +163 @@
 // already be setup to access the Array at the current cursor location
 // (beam, IF, Pol; see SDMemTable).   The desired Stokes
-// is specfied in the string from "I", "Q", "U", "V"
+// is specified in the string from "I", "Q", "U", "V"
   static casa::Array<casa::Float> getStokesSlice (casa::Array<casa::Float>& input, const casa::IPosition& start,
                                                   const casa::IPosition& end, const casa::String& stokes);
@@ -173 +173 @@
 
 
-// Compute value for STokes parameters by combining the raw correlation values
-// The output shape may change from the input shape according to
-// XX or YY          -> I
-// XX,YY             -> I
-// XX,YY,R(XY),I(XY) -> I,Q,U,V
-// Bool for mask
-// Float for TSys 
+// Compute value for auxilliary spectra when converting to Stokes parameters.
+// This combines the raw correlation values, according to what was combined
+// to actually convert to Stokes parameters.   The output shape may change from 
+// the input shape according to
+// XX or YY          -> I             (1 -> 1)
+// XX,YY             -> I             (2 -> 1)
+// XX,YY,R(XY),I(XY) -> I,Q,U,V       (4 -> 4)
+// 
+// It is meant for tSys (FLoat) and masks (Bool or uChar)
+// The input array must be of shape [nBeam,nIF,nPol,nChan]
   template <class T>
   static casa::Array<T> stokesData (casa::Array<T>& rawData, casa::Bool doLinear);
-  
+
+// Find the number of STokes parameters given the input number
+// of raw polarizations
+   static casa::uInt numberStokes (casa::uInt nPol) {casa::uInt nOut = nPol; if (nPol==2) nOut = 1; return nOut;};
 
 // Find the Stokes type for the given polarization axis (0,1,2,3)
-// You can ask for STokes or raw correltions (linear or circular)
+// You can ask for STokes or raw correlations (linear or circular)
    static casa::Stokes::StokesTypes convertStokes(casa::Int val, casa::Bool toStokes, 
                                                   casa::Bool linear);
+
+
+// These two functions are explicitly for the SDWriter
+//
+// Compute value for auxilliary spectra (Tsys and flags) when converting to Stokes parameters.
+// This combines the raw correlation values, according to what was combined
+// to actually convert to Stokes parameters.   The output shape may change from 
+// the input shape according to
+// XX or YY          -> I             (1 -> 1)
+// XX,YY             -> I             (2 -> 1)
+// XX,YY,R(XY),I(XY) -> I,Q,U,V       (4 -> 4)
+// 
+// The input array must be of shape [nChan,nPol]  (flags)
+// The input array must be of shape [nPol]        (tSys)
+//
+  static casa::Array<casa::uChar> computeStokesFlagsForWriter (casa::Array<casa::uChar>& rawData, casa::Bool doLinear)
+                                                         {return computeStokesDataForWriter(rawData,doLinear);};
+  static casa::Array<casa::Float> computeStokesTSysForWriter (casa::Array<casa::Float>& rawData, casa::Bool doLinear)
+                                                         {return computeStokesDataForWriter(rawData,doLinear);};
+
+// This function is explcitrly for SDWriter. It extracts the Beam and IF slice (specified in start
+// and end - must be length 1 in beam and IF) and flips the pol and channel axes.
+  static casa::Array<casa::Float> extractStokesForWriter (casa::Array<casa::Float>& in, const casa::IPosition& start, 
+                                                          const casa::IPosition& end);
 
 private:
@@ -197 +227 @@
                                             const casa::Array<casa::Bool>& in2) 
                                             {return in1&&in2;}
+  static casa::Array<casa::uChar> andArrays (const casa::Array<casa::uChar>& in1,
+                                            const casa::Array<casa::uChar>& in2);
+  
+  template <class T>
+  static casa::Array<T> computeStokesDataForWriter (casa::Array<T>& rawData, casa::Bool doLinear);
 };
 

trunk/src/SDPol2.cc

@@ -54 +54 @@
 //
 // Designed for use 
-//   Bool for mask
+//   Bool or uChar for mask
 //   Float for TSys
+//
+// Input arrays should be of shape
+//    [nBeam,nIF,nPol,nChan]
 //
 {
    T* t;
    DataType type = whatType(t);
-   AlwaysAssert(type==TpFloat || type==TpBool, AipsError);
+   AlwaysAssert(type==TpFloat || type==TpBool || type==TpUChar, AipsError);
 //
    IPosition shapeIn = rawData.shape();
+   const uInt nDim = shapeIn.nelements();
+   uInt polAxis = 0;
+   AlwaysAssert(nDim==asap::nAxes, AipsError);
+//
    uInt nPol = shapeIn(asap::PolAxis);
-   const uInt nDim = shapeIn.nelements();
    Array<T> stokesData;
 //
@@ -135 +141 @@
 }
 
+
+
+template <class T>
+Array<T> SDPolUtil::computeStokesDataForWriter (Array<T>& rawData, Bool doLinear)
+//
+// Generate data for each Stokes parameter from the
+// raw flags.  This is a lot of computational work and may
+// not be worth the effort.  This function is specifically
+// for the SDWriter 
+//
+// Designed for use 
+//   Bool or uChar for mask
+//   Float for TSys
+//
+// Input arrays should be of shape
+//    [nChan,nPol]   Bool/uChar
+//    [nPol]         Float      
+//
+{
+   T* t;
+   DataType type = whatType(t);
+   AlwaysAssert(type==TpFloat || type==TpBool || type==TpUChar, AipsError);
+//
+   IPosition shapeIn = rawData.shape();
+   const uInt nDim = shapeIn.nelements();
+   uInt polAxis = 0;
+   if (type==TpFloat) {
+      AlwaysAssert(nDim==1,AipsError);
+   } else {
+      AlwaysAssert(nDim==2,AipsError);
+      polAxis = 1;
+   }
+//
+   uInt nPol = shapeIn(polAxis);
+   Array<T> stokesData;
+//
+   IPosition start(nDim,0);
+   IPosition end(shapeIn-1);
+   IPosition shapeOut = shapeIn;
+//
+   if (doLinear) {
+      if (nPol==1) {
+         stokesData.resize(shapeOut);
+         stokesData = rawData;
+      } else if (nPol==2 || nPol==4) {
+
+// Set shape of output array
+
+         if (nPol==2) {
+            shapeOut(polAxis) = 1;
+         } else {
+            shapeOut(polAxis) = 4;
+         }
+         stokesData.resize(shapeOut);
+
+// Get reference slices and assign/compute
+
+         start(polAxis) = 0;
+         end(polAxis) = 0;
+         Array<T> M1In = rawData(start,end);
+//
+         start(polAxis) = 1;
+         end(polAxis) = 1;
+         Array<T> M2In = rawData(start,end);
+//
+         start(polAxis) = 0;
+         end(polAxis) = 0;
+         Array<T> M1Out = stokesData(start,end);         
+         M1Out = SDPolUtil::andArrays (M1In, M2In);         // I
+//
+         if (nPol==4) {   
+            start(polAxis) = 2;
+            end(polAxis) = 2;
+            Array<T> M3In = rawData(start,end);
+//
+            start(polAxis) = 3;
+            end(polAxis) = 3;
+            Array<T> M4In = rawData(start,end);
+//
+            start(polAxis) = 1;
+            end(polAxis) = 1;
+            Array<T> M2Out = stokesData(start,end);
+            M2Out = M1Out;                                  // Q
+//
+            start(polAxis) = 2;
+            end(polAxis) = 2;
+            Array<T> M3Out = stokesData(start,end);
+            M3Out = M3In;                                   // U
+//
+            start(polAxis) = 3;
+            end(polAxis) = 3;
+            Array<T> M4Out = stokesData(start,end);
+            M4Out = M4In;                                   // V
+         }
+      } else {
+         throw(AipsError("Can only handle 1,2 or 4 polarizations"));
+      }
+   } else {
+      throw (AipsError("Only implemented for Linear polarizations"));
+   }
+//
+   return stokesData;
+}
+
+