Changeset 1757 for branches/alma/src

Timestamp:

06/09/10 19:03:06 (14 years ago)

Author:

Kana Sugimoto

Message:

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.

Location:

branches/alma

Files:

• . (modified) (2 props)
• src/IndexedCompare.h (added)
• src/Makefile (modified) (3 diffs)
• src/MathUtils.cpp (modified) (2 diffs)
• src/MathUtils.h (modified) (3 diffs)
• src/RowAccumulator.cpp (modified) (3 diffs)
• src/RowAccumulator.h (modified) (1 diff)
• src/STAsciiWriter.cpp (modified) (2 diffs)
• src/STAtmosphere.cpp (added)
• src/STAtmosphere.h (added)
• src/STAttr.cpp (modified) (1 diff)
• src/STCoordinate.h (added)
• src/STDefs.h (modified) (2 diffs)
• src/STFITSImageWriter.cpp (modified) (3 diffs)
• src/STFiller.cpp (modified) (6 diffs)
• src/STFocus.cpp (modified) (11 diffs)
• src/STFocus.h (modified) (3 diffs)
• src/STFrequencies.cpp (modified) (5 diffs)
• src/STLineFinder.cpp (modified) (17 diffs)
• src/STLineFinder.h (modified) (2 diffs)
• src/STMath.cpp (modified) (18 diffs)
• src/STMath.h (modified) (5 diffs)
• src/STMathWrapper.h (modified) (3 diffs)
• src/STPol.h (modified) (3 diffs)
• src/STWriter.cpp (modified) (10 diffs)
• src/STWriter.h (modified) (2 diffs)
• src/Scantable.cpp (modified) (16 diffs)
• src/Scantable.h (modified) (11 diffs)
• src/ScantableWrapper.h (modified) (2 diffs)
• src/python_STAtmosphere.cpp (added)
• src/python_STCoordinate.cpp (added)
• src/python_Scantable.cpp (modified) (1 diff)
• src/python_asap.cpp (modified) (1 diff)
• src/python_asap.h (modified) (1 diff)

branches/alma/src/Makefile

@@ -123 +123 @@
              STAsciiWriter.o \
              STFITSImageWriter.o \
+             STAtmosphere.o \
              Scantable.o \
              Templates.o 
@@ -136 +137 @@
              python_LineCatalog.o \
              python_SrcType.o \
+             python_STAtmosphere.o \
+             python_STCoordinate.o \
              python_asap.o
 
@@ -172 +175 @@
              STWriter.h \
              STAsciiWriter.h \
-             STFITSImageWriter.h
+             STFITSImageWriter.h \
+             IndexedCompare.h \
+             STAtmosphere.h \
+             STCoordinate.h
 
 STATICCCLIB := libasap.a

branches/alma/src/MathUtils.cpp

@@ -39 +39 @@
 #include <scimath/Mathematics/MedianSlider.h>
 #include <casa/Exceptions/Error.h>
+
+#include <scimath/Fitting/LinearFit.h>
+#include <scimath/Functionals/Polynomial.h>
+#include <scimath/Mathematics/AutoDiff.h>
+
 
 #include "MathUtils.h"
@@ -183 +188 @@
   MedianSlider ms(hwidth);
   Slice sl(0, fwidth-1);
-  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl), 
+  Float medval = ms.add(const_cast<Vector<Float>& >(in)(sl),
                   const_cast<Vector<Bool>& >(flag)(sl));
   uInt n = in.nelements();
   for (uInt i=hwidth; i<(n-hwidth); ++i) {
     // add data value
-    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]); 
-    outflag[i] = (ms.nval() == 0);    
-  }
-  // replicate edge values from fisrt value with full width of values
+    out[i] = ms.add(in[i+hwidth], flag[i+hwidth]);
+    outflag[i] = (ms.nval() == 0);
+  }
+  // replicate edge values from first value with full width of values
   for (uInt i=0;i<hwidth;++i) {
     out[i] = out[hwidth];
-    outflag[i] = outflag[hwidth];    
+    outflag[i] = outflag[hwidth];
     out[n-1-i] = out[n-1-hwidth];
-    outflag[n-1-i] = outflag[n-1-hwidth];    
-  }
-}
+    outflag[n-1-i] = outflag[n-1-hwidth];
+  }
+}
+
+void mathutil::polyfit(Vector<Float>& out, Vector<Bool>& outmask,
+                       const Vector<Float>& in, const Vector<Bool>& mask,
+                       float width, int order)
+{
+  Int hwidth = Int(width+0.5);
+  Int fwidth = hwidth*2+1;
+  out.resize(in.nelements());
+  outmask.resize(mask.nelements());
+  LinearFit<Float> fitter;
+  Polynomial<Float> poly(order);
+  fitter.setFunction(poly);
+  Vector<Float> sigma(fwidth);
+  sigma = 1.0;
+  Vector<Float> parms;
+  Vector<Float> x(fwidth);
+  indgen(x);
+
+  uInt n = in.nelements();
+
+  for (uInt i=hwidth; i<(n-hwidth); ++i) {
+    // add data value
+    if (mask[i]) {
+      Slice sl(i-hwidth, fwidth);
+      const Vector<Float> &y = const_cast<Vector<Float>& >(in)(sl);
+      const Vector<Bool> &m = const_cast<Vector<Bool>& >(mask)(sl);
+      parms = fitter.fit(x, y, sigma, &m);
+
+      poly.setCoefficients(parms);
+      out[i] = poly(x[hwidth]);//cout << in[i] <<"->"<<out[i]<<endl;
+    } else {
+      out[i] = in[i];
+    }
+    outmask[i] = mask[i];
+  }
+  // replicate edge values from first value with full width of values
+  for (uInt i=0;i<hwidth;++i) {
+    out[i] = out[hwidth];
+    outmask[i] = outmask[hwidth];
+    out[n-1-i] = out[n-1-hwidth];
+    outmask[n-1-i] = outmask[n-1-hwidth];
+  }
+}

branches/alma/src/MathUtils.h

@@ -51 +51 @@
  * @param ignoreOther drop every second channel (NYI)
  */
-void hanning(casa::Vector<casa::Float>& out, 
+void hanning(casa::Vector<casa::Float>& out,
              casa::Vector<casa::Bool>& outmask,
-             const casa::Vector<casa::Float>& in, 
+             const casa::Vector<casa::Float>& in,
              const casa::Vector<casa::Bool>& mask,
              casa::Bool relaxed=casa::False,
@@ -60 +60 @@
 /**
  * Apply a running median to  a masked vector.
- * Edge solution:  The first and last hwidth channels will be replicated 
+ * Edge solution:  The first and last hwidth channels will be replicated
  * from the first/last value from a full window.
  * @param out the smoothed vector
@@ -68 +68 @@
  * @param hwidth half-width of the smoothing window
  */
- void runningMedian(casa::Vector<casa::Float>& out, 
+ void runningMedian(casa::Vector<casa::Float>& out,
+                   casa::Vector<casa::Bool>& outflag,
+                   const casa::Vector<casa::Float>& in,
+                   const casa::Vector<casa::Bool>& flag,
+                   float hwidth);
+
+ void polyfit(casa::Vector<casa::Float>& out,
                    casa::Vector<casa::Bool>& outmask,
-                   const casa::Vector<casa::Float>& in, 
+                   const casa::Vector<casa::Float>& in,
                    const casa::Vector<casa::Bool>& mask,
-                   float hwidth);
+                   float hwidth, int order);
 
 // Generate specified statistic

branches/alma/src/RowAccumulator.cpp

@@ -67 +67 @@
   Float totalweight = weight;
   MaskedArray<Float> data(v,m);
-  if ( weightType_ == asap::VAR ) {
+  if ( weightType_ == asap::W_VAR ) {
     if (m.nelements() == userMask_.nelements()) {
       Float fac = 1.0/variance(data(userMask_));
@@ -90 +90 @@
   Float w = 1.0;
   tsysSum_ += v[0];
-  if ( weightType_ == asap::TSYS  || weightType_ == asap::TINTSYS ) {
+  if ( weightType_ == asap::W_TSYS  || weightType_ == asap::W_TINTSYS ) {
     w /= (v[0]*v[0]);
   }
@@ -105 +105 @@
   Float w = 1.0;
   intervalSum_ += inter;
-  if ( weightType_ == asap::TINT || weightType_ == asap::TINTSYS ) {
+  if ( weightType_ == asap::W_TINT || weightType_ == asap::W_TINTSYS ) {
     w /= Float(inter);
   }

branches/alma/src/RowAccumulator.h

@@ -33 +33 @@
    * Constructor taking a weight type as defined in @ref STDefs
    */
-  explicit RowAccumulator(WeightType wt = asap::NONE);
+  explicit RowAccumulator(WeightType wt = asap::W_NONE);
 
  ~RowAccumulator();

branches/alma/src/STAsciiWriter.cpp

@@ -89 +89 @@
 
    String rootName(fileName);
-   if (rootName.length()==0) rootName = String("ascii");
 
   Block<String> cols(4);
@@ -104 +103 @@
     String dirtype = stable.getDirectionRefString();
     ostringstream onstr;
-    onstr << "SCAN" << rec.asuInt("SCANNO")
-    << "_CYCLE" << rec.asuInt("CYCLENO")
-    << "_BEAM" << rec.asuInt("BEAMNO")
-    << "_IF" << rec.asuInt("IFNO");
+
+    if (rootName.length()==0) {
+      rootName = String("ascii");
+    }
+    if (tab.nrow() > 1) {
+      if (stable.nscan() > 1) 
+        onstr << "_SCAN" << rec.asuInt("SCANNO");
+      if (stable.ncycle(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_CYCLE" << rec.asuInt("CYCLENO");
+      if (stable.nbeam(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_BEAM" << rec.asuInt("BEAMNO");
+      if (stable.nif(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_IF" << rec.asuInt("IFNO");
+    }
+
     String fName = rootName + String(onstr) + String(".txt");
     ofstream of(fName.chars(), ios::trunc);

branches/alma/src/STAttr.cpp

@@ -331 +331 @@
    TidGainElPoly_(2) = -3.219093e-4;
    
+   // 2009-09-15 - 13mm (22.2GHz) receiver
    ParkesGainElPoly_.resize(3);
-   ParkesGainElPoly_(0) = 0.296759e-1;
-   ParkesGainElPoly_(1) = -0.293124e-3;
-   ParkesGainElPoly_(2) = 0.264295e-6;
+   ParkesGainElPoly_(0) = -0.194031;
+   ParkesGainElPoly_(1) = 0.457724e-1;
+   ParkesGainElPoly_(2) = -0.438659e-3;
 }
 

branches/alma/src/STDefs.h

@@ -34 +34 @@
 
 namespace asap {
-  enum AxisNo { BeamAxis=0,
-                IFAxis,
-                PolAxis,
-                ChanAxis,
-                nAxes};
 
   enum Instrument {UNKNOWNINST=0,
@@ -53 +48 @@
   enum FeedPolType {UNKNOWNFEED, LINEAR, CIRCULAR, N_POL};
 
-  enum WeightType {NONE=0, VAR, TSYS, TINT, TINTSYS, N_WEIGHTTYPES};
-
-  enum TableType {MEMORY=0, PERSISTENT};
-
+  enum WeightType {W_NONE=0, W_VAR, W_TSYS, W_TINT, W_TINTSYS, W_N_WEIGHTTYPES};
 
 

branches/alma/src/STFITSImageWriter.cpp

@@ -121 +121 @@
   TableIterator iter(tab, cols);
   // Open data file
+
   while ( !iter.pastEnd() ) {
     Table t = iter.table();
@@ -127 +128 @@
     String dirtype = stable.getDirectionRefString();
     ostringstream onstr;
-    onstr << "SCAN" << rec.asuInt("SCANNO")
-    << "_CYCLE" << rec.asuInt("CYCLENO")
-    << "_BEAM" << rec.asuInt("BEAMNO")
-    << "_IF" << rec.asuInt("IFNO")
-    << "_POL" << rec.asuInt("POLNO");
+    if (rootName.length()==0) {
+      rootName = "fits";
+    }
+    if (tab.nrow() > 1) {
+      if (stable.nscan() > 1) 
+        onstr << "_SCAN" << rec.asuInt("SCANNO");
+      if (stable.ncycle(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_CYCLE" << rec.asuInt("CYCLENO");
+      if (stable.nbeam(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_BEAM" << rec.asuInt("BEAMNO");
+      if (stable.nif(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_IF" << rec.asuInt("IFNO");
+      if (stable.npol(rec.asuInt("SCANNO")) > 1) 
+        onstr << "_POL" << rec.asuInt("POLNO");
+    }
     String fileName = rootName + String(onstr) + String(".fits");
     int row0 = t.rowNumbers(tab)[0];
@@ -260 +271 @@
     }
     fits_close_file(fptr, &status);
-  
+    ostringstream oss;
+    oss << "Wrote " << fileName;
+    pushLog(String(oss));  
     //pushLog(String(oss));
     ++iter;

branches/alma/src/STFiller.cpp

@@ -264 +264 @@
     freqFrame = "REST";
   }
-  table_->frequencies().setFrame(freqFrame);
-     
+  // set both "FRAME" and "BASEFRAME"
+  table_->frequencies().setFrame(freqFrame, false);
+  table_->frequencies().setFrame(freqFrame,true);
+  //table_->focus().setParallactify(true);
 }
 
@@ -329 +331 @@
     if ( status != 0 ) break;
     n += 1;
-
     Regex filterrx(".*[SL|PA]$");
     Regex obsrx("^AT.+");
@@ -391 +392 @@
     /// @todo this has to change when nchan isn't global anymore
     //id = table_->frequencies().addEntry(Double(header_->nchan/2),
-    //                                        pksrec.refFreq, pksrec.freqInc);
+    //                                    pksrec.refFreq, pksrec.freqInc);
     if ( pksrec.nchan == 1 ) {
       id = table_->frequencies().addEntry(Double(0),
@@ -416 +417 @@
     RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
     *mweatheridCol = id;
+
     RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
-    id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
-                                  pksrec.focusRot);
+    id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi, 
+                                  pksrec.focusTan, pksrec.focusRot);
     *mfocusidCol = id;
     RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
@@ -426 +428 @@
     RecordFieldPtr<Float> elCol(rec, "ELEVATION");
     *elCol = pksrec.elevation;
-
-    RecordFieldPtr<Float> parCol(rec, "PARANGLE");
-    *parCol = pksrec.parAngle;
 
     RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
@@ -542 +541 @@
   }
 
-  // set frame keyword of FREQUENCIES table
+  // set FRAME and BASEFRAME keyword of FREQUENCIES table
   if ( header_->freqref != "TOPO" ) {
     table_->frequencies().setFrame( header_->freqref, false ) ;
+    table_->frequencies().setFrame( header_->freqref, true ) ;
   }
 

branches/alma/src/STFocus.cpp

@@ -34 +34 @@
 }
 
-asap::STFocus::STFocus( casa::Table tab ) : STSubTable(tab, name_)
+STFocus::STFocus( casa::Table tab ) : 
+  STSubTable(tab, name_)
 {
+  parangleCol_.attach(table_,"PARANGLE");
   rotationCol_.attach(table_,"ROTATION");
   axisCol_.attach(table_,"AXIS");
@@ -50 +52 @@
 }
 
-STFocus& asap::STFocus::operator =( const STFocus & other )
+STFocus& STFocus::operator =( const STFocus & other )
 {
   if (this != &other) {
     static_cast<STSubTable&>(*this) = other;
+    parangleCol_.attach(table_,"PARANGLE");
     rotationCol_.attach(table_,"ROTATION");
     axisCol_.attach(table_,"AXIS");
@@ -65 +68 @@
   return *this;
 }
-void asap::STFocus::setup( )
+void STFocus::setup( )
 {
   // add to base class table
+  table_.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
   table_.addColumn(ScalarColumnDesc<Float>("ROTATION"));
   table_.addColumn(ScalarColumnDesc<Float>("AXIS"));
@@ -76 +80 @@
   table_.addColumn(ScalarColumnDesc<Float>("XYPHASE"));
   table_.addColumn(ScalarColumnDesc<Float>("XYPHASEOFFSET"));
+  table_.rwKeywordSet().define("PARALLACTIFY", False);
 
   // new cached columns
+  parangleCol_.attach(table_,"PARANGLE");
   rotationCol_.attach(table_,"ROTATION");
   axisCol_.attach(table_,"AXIS");
@@ -88 +94 @@
 }
 
-uInt STFocus::addEntry( Float fax, Float ftan, Float frot, Float hand,
-                        Float user, Float mount,
-                        Float xyphase, Float xyphaseoffset)
+  uInt STFocus::addEntry( Float pa, Float fax, Float ftan, Float frot, Float hand,
+                          Float user, Float mount,
+                          Float xyphase, Float xyphaseoffset)
 {
-  Table result = table_( near(table_.col("ROTATION"), frot)
-                    && near(table_.col("AXIS"), fax)
-                    && near(table_.col("TAN"), ftan)
-                    && near(table_.col("HAND"), hand)
-                    && near(table_.col("USERPHASE"), user)
-                    && near(table_.col("MOUNT"), mount)
-                    && near(table_.col("XYPHASE"), xyphase)
-                    && near(table_.col("XYPHASEOFFSET"), xyphaseoffset)
-                    );
+  Table result = table_(  near(table_.col("PARANGLE"), pa)
+                          && near(table_.col("ROTATION"), frot)
+                          && near(table_.col("AXIS"), fax)
+                          && near(table_.col("TAN"), ftan)
+                          && near(table_.col("HAND"), hand)
+                          && near(table_.col("USERPHASE"), user)
+                          && near(table_.col("MOUNT"), mount)
+                          && near(table_.col("XYPHASE"), xyphase)
+                          && near(table_.col("XYPHASEOFFSET"), xyphaseoffset)
+                          );
   uInt resultid = 0;
   if ( result.nrow() > 0) {
@@ -113 +120 @@
       resultid++;
     }
+    parangleCol_.put(rno, pa);
     rotationCol_.put(rno, frot);
     axisCol_.put(rno, fax);
@@ -126 +134 @@
 }
 
-void asap::STFocus::getEntry( Float& rotation, Float& angle, Float& ftan,
-                              Float& hand, Float& user, Float& mount,
-                              Float& xyphase, Float& xyphaseoffset,
-                              uInt id) const
+  void STFocus::getEntry( Float& pa, Float& rotation, Float& angle, Float& ftan,
+                                Float& hand, Float& user, Float& mount,
+                                Float& xyphase, Float& xyphaseoffset,
+                                uInt id) const
 {
   Table t = table_(table_.col("ID") == Int(id) );
@@ -138 +146 @@
   // get first row - there should only be one matching id
   const TableRecord& rec = row.get(0);
+  pa = rec.asFloat("PARANGLE");
   rotation = rec.asFloat("ROTATION");
   angle = rec.asFloat("AXIS");
@@ -149 +158 @@
 
 
-casa::Float asap::STFocus::getTotalFeedAngle( casa::uInt id ) const
+casa::Float STFocus::getTotalAngle( casa::uInt id ) const
 {
   Float total = 0.0f;
   Table t = table_(table_.col("ID") == Int(id) );
   if (t.nrow() == 0 ) {
-    throw(AipsError("STFocus::getEntry - id out of range"));
+    throw(AipsError("STFocus::getTotalAngle - id out of range"));
+  }
+  if (table_.keywordSet().asBool("PARALLACTIFY")) {
+    return 0.0f;
   }
   ROTableRow row(t);
   // get first row - there should only be one matching id
   const TableRecord& rec = row.get(0);
+  total += rec.asFloat("PARANGLE");  
   total += rec.asFloat("ROTATION");
   total += rec.asFloat("USERPHASE");
@@ -164 +177 @@
   return total;
 }
-}
 
-casa::Float asap::STFocus::getFeedHand( casa::uInt id ) const
+
+casa::Float STFocus::getFeedHand( casa::uInt id ) const
 {
   Table t = table_(table_.col("ID") == Int(id) );
@@ -177 +190 @@
 }
 
+void STFocus::setParallactify(bool istrue) {
+  table_.rwKeywordSet().define("PARALLACTIFY", Bool(istrue));
+}
+
+}

branches/alma/src/STFocus.h

@@ -37 +37 @@
   STFocus& operator=(const STFocus& other);
 
-  casa::uInt addEntry( casa::Float faxis, casa::Float ftan,
+  casa::uInt addEntry( casa::Float pa, casa::Float faxis, casa::Float ftan,
                        casa::Float frot, casa::Float hand=1.0f,
                        casa::Float mount=0.0f, casa::Float user=0.0f,
-                       casa::Float xyphase=0.0f, casa::Float xyphaseoffset=0.0f);
+                       casa::Float xyphase=0.0f, 
+                       casa::Float xyphaseoffset=0.0f);
 
-  void getEntry( casa::Float& fax, casa::Float& ftan,
+  void getEntry( casa::Float& pa, casa::Float& fax, casa::Float& ftan,
                  casa::Float& frot, casa::Float& hand,
                  casa::Float& mount, casa::Float& user,
@@ -48 +49 @@
                  casa::uInt id) const;
 
-  casa::Float getTotalFeedAngle(casa::uInt id) const;
+  casa::Float getTotalAngle(casa::uInt id) const;
+
+  casa::Float getParAngle(casa::uInt id) const {
+    return parangleCol_(id);
+  }
   casa::Float getFeedHand(casa::uInt id) const;
+
+  void setParallactify(bool istrue=false);
 
   const casa::String& name() const { return name_; }
@@ -57 +64 @@
   static const casa::String name_;
   casa::ScalarColumn<casa::Float> rotationCol_, axisCol_,
-                                  tanCol_,handCol_,
-                                  mountCol_,userCol_,
-                                  xyphCol_,xyphoffCol_;
+    tanCol_,handCol_, parangleCol_,
+    mountCol_,userCol_, xyphCol_,xyphoffCol_,;
 };
 

branches/alma/src/STFrequencies.cpp

@@ -170 +170 @@
 **/
 SpectralCoordinate
-  asap::STFrequencies::getSpectralCoordinate( const MDirection& md,
+  STFrequencies::getSpectralCoordinate( const MDirection& md,
                                               const MPosition& mp,
                                               const MEpoch& me,
@@ -242 +242 @@
 {
   Vector<Float> offset(1,0.0);
-  Vector<Float> factors(1,1.0/width);
+  Vector<Float> factors(1,width);
   Vector<Int> newshape;
   CoordinateSystem csys;
@@ -317 +317 @@
         << rec.asDouble("REFVAL") << setw(7)
         << rec.asDouble("REFPIX")
-        << setw(12)
+        << setw(15)
         << rec.asDouble("INCREMENT");
   }
@@ -324 +324 @@
     int f = outstr.find_first_not_of(' ');
     int l = outstr.find_last_not_of(' ', outstr.size());
-    if (f < 0) { 
-      f = 0; 
-    }
-    if ( l < f  || l < f ) { 
+    if (f < 0) {
+      f = 0;
+    }
+    if ( l < f  || l < f ) {
       l = outstr.size();
     }
@@ -439 +439 @@
 }
 
-void STFrequencies::shiftRefPix(int npix, uInt id) 
+void STFrequencies::shiftRefPix(int npix, uInt id)
 {
   Table t = table_(table_.col("ID") == Int(id) );

branches/alma/src/STLineFinder.cpp

@@ -34 +34 @@
 #include "STLineFinder.h"
 #include "STFitter.h"
+#include "IndexedCompare.h"
 
 // STL
@@ -110 +111 @@
 
 protected:
-   // supplementary function to control running mean calculations.
-   // It adds a specified channel to the running mean box and
+   // supplementary function to control running mean/median calculations.
+   // It adds a specified channel to the running box and
    // removes (ch-maxboxnchan+1)'th channel from there
    // Channels, for which the mask is false or index is beyond the
@@ -152 +153 @@
                                            // last point of the detected line
                                            //
+   bool itsUseMedian;                      // true if median statistics is used
+                                           // to determine the noise level, otherwise
+                                           // it is the mean of the lowest 80% of deviations
+                                           // (default)
+   int itsNoiseSampleSize;                 // sample size used to estimate the noise statistics
+                                           // Negative value means the whole spectrum is used (default)
 public:
 
@@ -157 +164 @@
    LFAboveThreshold(std::list<pair<int,int> > &in_lines,
                     int in_min_nchan = 3,
-                    casa::Float in_threshold = 5) throw();
+                    casa::Float in_threshold = 5,
+                    bool use_median = false,
+                    int noise_sample_size = -1) throw();
    virtual ~LFAboveThreshold() throw();
 
@@ -200 +209 @@
 ///////////////////////////////////////////////////////////////////////////////
 
+///////////////////////////////////////////////////////////////////////////////
+//
+// LFNoiseEstimator  a helper class designed to estimate off-line variance
+//                   using statistics depending on the distribution of 
+//                   values (e.g. like a median)
+//
+//                   Two statistics are supported: median and an average of
+//                   80% of smallest values. 
+//
+
+struct LFNoiseEstimator {
+   // construct an object
+   // size - maximum sample size. After a size number of elements is processed
+   // any new samples would cause the algorithm to drop the oldest samples in the
+   // buffer.
+   explicit LFNoiseEstimator(size_t size);
+
+   // add a new sample 
+   // in - the new value
+   void add(float in);
+
+   // median of the distribution
+   float median() const;
+
+   // mean of lowest 80% of the samples
+   float meanLowest80Percent() const;
+
+   // return true if the buffer is full (i.e. statistics are representative)
+   inline bool filledToCapacity() const { return itsBufferFull;}
+
+protected:
+   // update cache of sorted indices
+   // (it is assumed that itsSampleNumber points to the newly
+   // replaced element)
+   void updateSortedCache() const;
+
+   // build sorted cache from the scratch
+   void buildSortedCache() const;
+
+   // number of samples accumulated so far
+   // (can be less than the buffer size)
+   size_t numberOfSamples() const;
+
+   // this helper method builds the cache if
+   // necessary using one of the methods
+   void fillCacheIfNecessary() const;
+
+private:
+   // buffer with samples (unsorted)
+   std::vector<float> itsVariances;
+   // current sample number (<=itsVariances.size())
+   size_t itsSampleNumber;
+   // true, if the buffer all values in the sample buffer are used
+   bool itsBufferFull;
+   // cached indices into vector of samples
+   mutable std::vector<size_t> itsSortedIndices;
+   // true if any of the statistics have been obtained at least
+   // once. This flag allows to implement a more efficient way of
+   // calculating statistics, if they are needed at once and not 
+   // after each addition of a new element
+   mutable bool itsStatisticsAccessed;
+};
+
+//
+///////////////////////////////////////////////////////////////////////////////
+
+
 } // namespace asap
+
+///////////////////////////////////////////////////////////////////////////////
+//
+// LFNoiseEstimator  a helper class designed to estimate off-line variance
+//                   using statistics depending on the distribution of 
+//                   values (e.g. like a median)
+//
+//                   Two statistics are supported: median and an average of
+//                   80% of smallest values. 
+//
+
+// construct an object
+// size - maximum sample size. After a size number of elements is processed
+// any new samples would cause the algorithm to drop the oldest samples in the
+// buffer.
+LFNoiseEstimator::LFNoiseEstimator(size_t size) : itsVariances(size),
+     itsSampleNumber(0), itsBufferFull(false), itsSortedIndices(size),
+     itsStatisticsAccessed(false)
+{
+   AlwaysAssert(size>0,AipsError);
+} 
+
+
+// add a new sample 
+// in - the new value
+void LFNoiseEstimator::add(float in)
+{
+   if (isnan(in)) {
+       // normally it shouldn't happen
+       return;
+   }
+   itsVariances[itsSampleNumber] = in;
+
+   if (itsStatisticsAccessed) {
+       // only do element by element addition if on-the-fly 
+       // statistics are needed
+       updateSortedCache();
+   }
+
+   // advance itsSampleNumber now
+   ++itsSampleNumber;
+   if (itsSampleNumber == itsVariances.size()) {
+       itsSampleNumber = 0;
+       itsBufferFull = true;
+   } 
+   AlwaysAssert(itsSampleNumber<itsVariances.size(),AipsError);
+}
+
+// number of samples accumulated so far
+// (can be less than the buffer size)
+size_t LFNoiseEstimator::numberOfSamples() const
+{
+  // the number of samples accumulated so far may be less than the
+  // buffer size
+  const size_t nSamples = itsBufferFull ? itsVariances.size(): itsSampleNumber;
+  AlwaysAssert( (nSamples > 0) && (nSamples <= itsVariances.size()), AipsError);
+  return nSamples;
+}
+
+// this helper method builds the cache if
+// necessary using one of the methods
+void LFNoiseEstimator::fillCacheIfNecessary() const
+{
+  if (!itsStatisticsAccessed) {
+      if ((itsSampleNumber!=0) || itsBufferFull) {
+          // build the whole cache efficiently 
+          buildSortedCache();
+      } else {
+          updateSortedCache();
+      }
+      itsStatisticsAccessed = true;
+  } // otherwise, it is updated in 'add' using on-the-fly method
+}
+
+// median of the distribution
+float LFNoiseEstimator::median() const
+{
+  fillCacheIfNecessary();
+  // the number of samples accumulated so far may be less than the
+  // buffer size
+  const size_t nSamples = numberOfSamples();
+  const size_t medSample = nSamples / 2;
+  AlwaysAssert(medSample < itsSortedIndices.size(), AipsError);
+  return itsVariances[itsSortedIndices[medSample]];
+}
+
+// mean of lowest 80% of the samples
+float LFNoiseEstimator::meanLowest80Percent() const
+{
+  fillCacheIfNecessary();
+  // the number of samples accumulated so far may be less than the
+  // buffer size
+  const size_t nSamples = numberOfSamples();
+  float result = 0;
+  size_t numpt=size_t(0.8*nSamples);
+  if (!numpt) {
+      numpt=nSamples; // no much else left,
+                     // although it is very inaccurate
+  }
+  AlwaysAssert( (numpt > 0) && (numpt<itsSortedIndices.size()), AipsError);
+  for (size_t ch=0; ch<numpt; ++ch) {
+       result += itsVariances[itsSortedIndices[ch]];
+  }
+  result /= float(numpt);
+  return result;
+}
+
+// update cache of sorted indices
+// (it is assumed that itsSampleNumber points to the newly
+// replaced element)
+void LFNoiseEstimator::updateSortedCache() const
+{
+  // the number of samples accumulated so far may be less than the
+  // buffer size
+  const size_t nSamples = numberOfSamples();
+
+  if (itsBufferFull) {
+      // first find the index of the element which is being replaced
+      size_t index = nSamples;
+      for (size_t i=0; i<nSamples; ++i) {
+           AlwaysAssert(i < itsSortedIndices.size(), AipsError);
+           if (itsSortedIndices[i] == itsSampleNumber) {
+               index = i;
+               break;
+           }
+      }
+      AlwaysAssert( index < nSamples, AipsError);
+
+      const vector<size_t>::iterator indStart = itsSortedIndices.begin();
+      // merge this element with preceeding block first
+      if (index != 0) {
+          // merge indices on the basis of variances
+          inplace_merge(indStart,indStart+index,indStart+index+1, 
+                        indexedCompare<size_t>(itsVariances.begin()));
+      }
+      // merge with the following block
+      if (index + 1 != nSamples) {
+          // merge indices on the basis of variances
+          inplace_merge(indStart,indStart+index+1,indStart+nSamples, 
+                        indexedCompare<size_t>(itsVariances.begin()));
+      }
+  } else {
+      // itsSampleNumber is the index of the new element
+      AlwaysAssert(itsSampleNumber < itsSortedIndices.size(), AipsError);
+      itsSortedIndices[itsSampleNumber] = itsSampleNumber;
+      if (itsSampleNumber >= 1) {
+          // we have to place this new sample in
+          const vector<size_t>::iterator indStart = itsSortedIndices.begin();
+          // merge indices on the basis of variances
+          inplace_merge(indStart,indStart+itsSampleNumber,indStart+itsSampleNumber+1, 
+                        indexedCompare<size_t>(itsVariances.begin()));
+      }
+  }
+}
+
+// build sorted cache from the scratch
+void LFNoiseEstimator::buildSortedCache() const
+{
+  // the number of samples accumulated so far may be less than the
+  // buffer size
+  const size_t nSamples = numberOfSamples();
+  AlwaysAssert(nSamples <= itsSortedIndices.size(), AipsError);
+  for (size_t i=0; i<nSamples; ++i) {
+       itsSortedIndices[i]=i;
+  }
+
+  // sort indices, but check the array of variances
+  const vector<size_t>::iterator indStart = itsSortedIndices.begin();
+  stable_sort(indStart,indStart+nSamples, indexedCompare<size_t>(itsVariances.begin()));
+}
+
+//
+///////////////////////////////////////////////////////////////////////////////
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -275 +524 @@
 }
 
-// supplementary function to control running mean calculations.
-// It adds a specified channel to the running mean box and
+// supplementary function to control running mean/median calculations.
+// It adds a specified channel to the running box and
 // removes (ch-max_box_nchan+1)'th channel from there
 // Channels, for which the mask is false or index is beyond the
@@ -339 +588 @@
                 (meanch2-square(meanch));
       linmean=coeff*(Float(cur_channel)-meanch)+mean;
-      linvariance=sqrt(sumf2/Float(box_chan_cntr)-square(mean)-
-                    square(coeff)*(meanch2-square(meanch)));
+      linvariance=sumf2/Float(box_chan_cntr)-square(mean)-
+                    square(coeff)*(meanch2-square(meanch));
+      if (linvariance<0.) {
+          // this shouldn't happen normally, but could be due to round-off error
+          linvariance = 0; 
+      } else {
+          linvariance = sqrt(linvariance);
+      }
   }
   need2recalculate=False;
@@ -351 +606 @@
 ///////////////////////////////////////////////////////////////////////////////
 //
-// LFAboveThreshold - a running mean algorithm for line detection
+// LFAboveThreshold - a running mean/median algorithm for line detection
 //
 //
@@ -359 +614 @@
 LFAboveThreshold::LFAboveThreshold(std::list<pair<int,int> > &in_lines,
                                    int in_min_nchan,
-                                   casa::Float in_threshold) throw() :
+                                   casa::Float in_threshold,
+                                   bool use_median,
+                                   int noise_sample_size) throw() :
              min_nchan(in_min_nchan), threshold(in_threshold),
-             lines(in_lines), running_box(NULL) {}
+             lines(in_lines), running_box(NULL), itsUseMedian(use_median),
+             itsNoiseSampleSize(noise_sample_size) {}
 
 LFAboveThreshold::~LFAboveThreshold() throw()
@@ -474 +732 @@
       running_box=new RunningBox(spectrum,mask,edge,max_box_nchan);
 
-
       // determine the off-line variance first
       // an assumption made: lines occupy a small part of the spectrum
 
-      std::vector<float> variances(edge.second-edge.first);
-      DebugAssert(variances.size(),AipsError);
-
-      for (;running_box->haveMore();running_box->next())
-           variances[running_box->getChannel()-edge.first]=
-                                running_box->getLinVariance();
-
-      // in the future we probably should do a proper Chi^2 estimation
-      // now a simple 80% of smaller values will be used.
-      // it may degrade the performance of the algorithm for weak lines
-      // due to a bias of the Chi^2 distribution.
-      stable_sort(variances.begin(),variances.end());
-
-      Float offline_variance=0;
-      uInt offline_cnt=uInt(0.8*variances.size());
-      if (!offline_cnt) offline_cnt=variances.size(); // no much else left,
-                                    // although it is very inaccurate
-      for (uInt n=0;n<offline_cnt;++n)
-           offline_variance+=variances[n];
-      offline_variance/=Float(offline_cnt);
+      const size_t noiseSampleSize = itsNoiseSampleSize<0 ? size_t(edge.second-edge.first) : 
+                      std::min(size_t(itsNoiseSampleSize), size_t(edge.second-edge.first));
+      DebugAssert(noiseSampleSize,AipsError);
+      const bool globalNoise = (size_t(edge.second - edge.first) == noiseSampleSize);
+      LFNoiseEstimator ne(noiseSampleSize);
+
+      for (;running_box->haveMore();running_box->next()) {
+           ne.add(running_box->getLinVariance());
+           if (ne.filledToCapacity()) {
+               break;
+           }
+      }
+
+      Float offline_variance = -1; // just a flag that it is unset
+            
+      if (globalNoise) {
+          offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
+      }
 
       // actual search algorithm
@@ -510 +765 @@
                                  running_box->next()) {
            const int ch=running_box->getChannel();
-           if (running_box->getNumberOfBoxPoints()>=minboxnchan)
+           if (!globalNoise) {
+               // add a next point for a local noise estimate
+               ne.add(running_box->getLinVariance());
+           }    
+           if (running_box->getNumberOfBoxPoints()>=minboxnchan) {
+               if (!globalNoise) {
+                   offline_variance = itsUseMedian ? ne.median() : ne.meanLowest80Percent();
+               }
+               AlwaysAssert(offline_variance>0.,AipsError);
                processChannel(mask[ch] && (fabs(running_box->aboveMean()) >=
                   threshold*offline_variance), mask);
-           else processCurLine(mask); // just finish what was accumulated before
+           } else processCurLine(mask); // just finish what was accumulated before
 
            signs[ch]=getAboveMeanSign();
-           // os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
-           // threshold*offline_variance<<endl;
-
-           const Float buf=running_box->aboveMean();
-           if (buf>0) signs[ch]=1;
-           else if (buf<0) signs[ch]=-1;
-           else if (buf==0) signs[ch]=0;
-           //os<<ch<<" "<<spectrum[ch]<<" "<<running_box->getLinMean()<<" "<<
-           //             threshold*offline_variance<<endl;
+            //os<<ch<<" "<<spectrum[ch]<<" "<<fabs(running_box->aboveMean())<<" "<<
+            //threshold*offline_variance<<endl;
       }
       if (lines.size())
@@ -649 +905 @@
 //              Setting a very large value doesn't usually provide
 //              valid detections.
-// in_box_size  the box size for running mean calculation. Default is
+// in_box_size  the box size for running mean/median calculation. Default is
 //              1./5. of the whole spectrum size
+// in_noise_box the box size for off-line noise estimation (if working with
+//              local noise. Negative value means use global noise estimate
+//              Default is -1 (i.e. estimate using the whole spectrum)
+// in_median    true if median statistics is used as opposed to average of
+//              the lowest 80% of deviations (default)
 void STLineFinder::setOptions(const casa::Float &in_threshold,
                               const casa::Int &in_min_nchan,
                               const casa::Int &in_avg_limit,
-                              const casa::Float &in_box_size) throw()
+                              const casa::Float &in_box_size,
+                              const casa::Float &in_noise_box,
+                              const casa::Bool &in_median) throw()
 {
   threshold=in_threshold;
@@ -660 +923 @@
   avg_limit=in_avg_limit;
   box_size=in_box_size;
+  itsNoiseBox = in_noise_box;
+  itsUseMedian = in_median;
 }
 
@@ -683 +948 @@
                 const casa::uInt &whichRow) throw(casa::AipsError)
 {
-  //const int minboxnchan=4;
   if (scan.null())
       throw AipsError("STLineFinder::findLines - a scan should be set first,"
@@ -700 +964 @@
       throw AipsError("STLineFinder::findLines - in_scan and in_mask have different"
             "number of spectral channels.");
+
+  // taking flagged channels into account
+  vector<bool> flaggedChannels = scan->getMask(whichRow);
+  if (flaggedChannels.size()) {
+      // there is a mask set for this row
+      if (flaggedChannels.size() != mask.nelements()) {
+          throw AipsError("STLineFinder::findLines - internal inconsistency: number of mask elements do not match the number of channels");
+      }
+      for (size_t ch = 0; ch<mask.nelements(); ++ch) {
+           mask[ch] &= flaggedChannels[ch];
+      }
+  }
+
   // number of elements in in_edge
   if (in_edge.size()>2)
@@ -732 +1009 @@
       throw AipsError("STLineFinder::findLines - box_size is too small");
 
+  // number of elements in the sample for noise estimate
+  const int noise_box = itsNoiseBox<0 ? -1 : int(nchan * itsNoiseBox); 
+
+  if ((noise_box!= -1) and (noise_box<2))
+      throw AipsError("STLineFinder::findLines - noise_box is supposed to be at least 2 elements");
+
   spectrum.resize();
   spectrum = Vector<Float>(scan->getSpectrum(whichRow));
@@ -755 +1038 @@
      try {
          // line find algorithm
-         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold);
+         LFAboveThreshold lfalg(new_lines,avg_factor*min_nchan, threshold, itsUseMedian,noise_box);
          lfalg.findLines(spectrum,temp_mask,edge,max_box_nchan);
          signs.resize(lfalg.getSigns().nelements());

branches/alma/src/STLineFinder.h

@@ -150 +150 @@
    // in_box_size  the box size for running mean calculation. Default is
    //              1./5. of the whole spectrum size
+   // in_noise_box the box size for off-line noise estimation (if working with
+   //              local noise. Negative value means use global noise estimate
+   //              Default is -1 (i.e. estimate using the whole spectrum)
+   // in_median    true if median statistics is used as opposed to average of
+   //              the lowest 80% of deviations (default)
    void setOptions(const casa::Float &in_threshold=sqrt(3.),
                    const casa::Int &in_min_nchan=3,
                    const casa::Int &in_avg_limit=8,
-                   const casa::Float &in_box_size=0.2) throw();
+                   const casa::Float &in_box_size=0.2,
+                   const casa::Float &in_noise_box=-1.,
+                   const casa::Bool &in_median = casa::False) throw();
 
    // set the scan to work with (in_scan parameter)
@@ -241 +248 @@
    // a buffer for the spectrum
    mutable casa::Vector<casa::Float>  spectrum;
+
+   // the box size for off-line noise estimation (if working with
+   // local noise. Negative value means use global noise estimate
+   // Default is -1 (i.e. estimate using the whole spectrum)
+   casa::Float itsNoiseBox;
+
+   // true if median statistics is used as opposed to average of
+   // the lowest 80% of deviations (default)
+   casa::Bool itsUseMedian;
 };
 

branches/alma/src/STMath.cpp

@@ -52 +52 @@
 #include "RowAccumulator.h"
 #include "STAttr.h"
+#include "STSelector.h"
+
 #include "STMath.h"
-#include "STSelector.h"
-
 using namespace casa;
 
@@ -740 +740 @@
 }
 
-CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left, 
-                                            const CountedPtr<Scantable>& right, 
+CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
+                                            const CountedPtr<Scantable>& right,
                                             const std::string& mode)
 {
@@ -1379 +1379 @@
 {
 
-  
+
   STSelector sel;
   CountedPtr< Scantable > ws = getScantable(s, false);
@@ -1518 +1518 @@
   // for each row
   // assume that the data order are same between sig and ref
-  RowAccumulator acc( asap::TINTSYS ) ; 
+  RowAccumulator acc( asap::W_TINTSYS ) ; 
   for ( int i = 0 ; i < sig->nrow() ; i++ ) {
     // get values
@@ -1952 +1952 @@
   // initialize the lookup table if necessary
   if ( lookup.empty() ) {
-    lookup["NONE"]   = asap::NONE;
-    lookup["TINT"] = asap::TINT;
-    lookup["TINTSYS"]  = asap::TINTSYS;
-    lookup["TSYS"]  = asap::TSYS;
-    lookup["VAR"]  = asap::VAR;
+    lookup["NONE"]   = asap::W_NONE;
+    lookup["TINT"] = asap::W_TINT;
+    lookup["TINTSYS"]  = asap::W_TINTSYS;
+    lookup["TSYS"]  = asap::W_TSYS;
+    lookup["VAR"]  = asap::W_VAR;
   }
 
@@ -2000 +2000 @@
     String msg;
     if ( nc > 0 ) {
-      ppoly = new Polynomial<Float>(nc);
+      ppoly = new Polynomial<Float>(nc-1);
       coeff = coeffs;
       msg = String("user");
@@ -2006 +2006 @@
       STAttr sdAttr;
       coeff = sdAttr.gainElevationPoly(inst);
-      ppoly = new Polynomial<Float>(3);
+      ppoly = new Polynomial<Float>(coeff.nelements()-1);
       msg = String("built in");
     }
@@ -2144 +2144 @@
     if (d < 0) {
       Instrument inst =
-        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), 
+        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
                                   True);
       STAttr sda;
@@ -2207 +2207 @@
 
 CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
-                                         float tau )
+                                         const std::vector<float>& tau )
 {
   CountedPtr< Scantable > out = getScantable(in, false);
 
-  Table tab = out->table();
-  ROScalarColumn<Float> elev(tab, "ELEVATION");
-  ArrayColumn<Float> specCol(tab, "SPECTRA");
-  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
-  ArrayColumn<Float> tsysCol(tab, "TSYS");
-  for ( uInt i=0; i<tab.nrow(); ++i) {
-    Float zdist = Float(C::pi_2) - elev(i);
-    Float factor = exp(tau/cos(zdist));
-    MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
-    ma *= factor;
-    specCol.put(i, ma.getArray());
-    flagCol.put(i, flagsFromMA(ma));
-    Vector<Float> tsys;
-    tsysCol.get(i, tsys);
-    tsys *= factor;
-    tsysCol.put(i, tsys);
+  Table outtab = out->table();
+
+  const uInt ntau = uInt(tau.size());
+  std::vector<float>::const_iterator tauit = tau.begin();
+  AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
+               AipsError);
+  TableIterator iiter(outtab, "IFNO");
+  while ( !iiter.pastEnd() ) {
+    Table itab = iiter.table();
+    TableIterator piter(outtab, "POLNO");
+    while ( !piter.pastEnd() ) {
+      Table tab = piter.table();
+      ROScalarColumn<Float> elev(tab, "ELEVATION");
+      ArrayColumn<Float> specCol(tab, "SPECTRA");
+      ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+      ArrayColumn<Float> tsysCol(tab, "TSYS");
+      for ( uInt i=0; i<tab.nrow(); ++i) {
+        Float zdist = Float(C::pi_2) - elev(i);
+        Float factor = exp(*tauit/cos(zdist));
+        MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
+        ma *= factor;
+        specCol.put(i, ma.getArray());
+        flagCol.put(i, flagsFromMA(ma));
+        Vector<Float> tsys;
+        tsysCol.get(i, tsys);
+        tsys *= factor;
+        tsysCol.put(i, tsys);
+      }
+      if (ntau == in->nif()*in->npol() ) {
+        tauit++;
+      }
+      piter++;
+    }
+    if (ntau >= in->nif() ) {
+      tauit++;
+    }
+    iiter++;
   }
   return out;
@@ -2233 +2254 @@
 CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
                                              const std::string& kernel,
-                                             float width )
+                                             float width, int order)
 {
   CountedPtr< Scantable > out = getScantable(in, false);
@@ -2254 +2275 @@
       mathutil::runningMedian(specout, maskout, spec , mask, width);
       convertArray(flag, maskout);
+    } else if ( kernel == "poly" ) {
+      mathutil::polyfit(specout, maskout, spec, !mask, width, order);
+      convertArray(flag, !maskout);
     }
     flagCol.put(i, flag);
@@ -2262 +2286 @@
 
 CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
-                                        const std::string& kernel, float width )
-{
-  if (kernel == "rmedian"  || kernel == "hanning") {
-    return smoothOther(in, kernel, width);
+                                        const std::string& kernel, float width,
+                                        int order)
+{
+  if (kernel == "rmedian"  || kernel == "hanning" || kernel == "poly") {
+    return smoothOther(in, kernel, width, order);
   }
   CountedPtr< Scantable > out = getScantable(in, false);
@@ -2351 +2376 @@
           id = out->molecules().addEntry(rf, name, fname);
           molidcol.put(k, id);
-          Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
-          (*it)->focus().getEntry(fax, ftan, frot, fhand,
+          Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
+          (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand,
                                   fmount,fuser, fxy, fxyp,
                                   rec.asuInt("FOCUS_ID"));
-          id = out->focus().addEntry(fax, ftan, frot, fhand,
+          id = out->focus().addEntry(fpa, fax, ftan, frot, fhand,
                                      fmount,fuser, fxy, fxyp);
           focusidcol.put(k, id);
@@ -2399 +2424 @@
   cols[3] = String("CYCLENO");
   TableIterator iter(tout, cols);
-  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_, 
+  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
                                                out->getPolType() );
   while (!iter.pastEnd()) {
@@ -2405 +2430 @@
     ArrayColumn<Float> speccol(t, "SPECTRA");
     ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
-    ScalarColumn<Float> parancol(t, "PARANGLE");
     Matrix<Float> pols(speccol.getColumn());
     try {
       stpol->setSpectra(pols);
-      Float fang,fhand,parang;
-      fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
+      Float fang,fhand;
+      fang = in->focusTable_.getTotalAngle(focidcol(0));
       fhand = in->focusTable_.getFeedHand(focidcol(0));
-      parang = parancol(0);
-      /// @todo re-enable this
-      // disable total feed angle to support paralactifying Caswell style
-      stpol->setPhaseCorrections(parang, -parang, fhand);
+      stpol->setPhaseCorrections(fang, fhand);
       // use a member function pointer in STPol.  This only works on
       // the STPol pointer itself, not the Counted Pointer so
@@ -2681 +2702 @@
       uInt row = tab.rowNumbers()[0];
       stpol->setSpectra(in->getPolMatrix(row));
-      Float fang,fhand,parang;
-      fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
+      Float fang,fhand;
+      fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
       fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
-      parang = in->paraCol_(row);
-      /// @todo re-enable this
-      // disable total feed angle to support paralactifying Caswell style
-      stpol->setPhaseCorrections(parang, -parang, fhand);
+      stpol->setPhaseCorrections(fang, fhand);
       Int npolout = 0;
       for (uInt i=0; i<tab.nrow(); ++i) {
@@ -2737 +2755 @@
 CountedPtr< Scantable >
   asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
-                          double frequency, double width )
+                         double start, double end,
+                         const std::string& mode)
 {
   CountedPtr< Scantable > out = getScantable(in, false);
@@ -2755 +2774 @@
       Vector<Float> spec = specCol(i);
       Vector<uChar> flag = flagCol(i);
-      Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
-      Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
+      int fstart = -1;
+      int fend = -1;
       for (unsigned int k=0; k < flag.nelements(); ++k ) {
         if (flag[k] > 0) {
-          spec[k] = 0.0;
+          fstart = k;
+          while (flag[k] > 0 && k < flag.nelements()) {
+            fend = k;
+            k++;
+          }
         }
+        Float interp = 0.0;
+        if (fstart-1 > 0 ) {
+          interp = spec[fstart-1];
+          if (fend+1 < spec.nelements()) {
+            interp = (interp+spec[fend+1])/2.0;
+          }
+        } else {
+          interp = spec[fend+1];
+        }
+        if (fstart > -1 && fend > -1) {
+          for (int j=fstart;j<=fend;++j) {
+            spec[j] = interp;
+          }
+        }
+        fstart =-1;
+        fend = -1;
       }
       Vector<Complex> lags;
       ffts.fft0(lags, spec);
-      Int start =  max(0, lag0);
-      Int end =  min(Int(lags.nelements()-1), lag1);
-      if (start == end) {
-        lags[start] = Complex(0.0);
+      Int lag0(start+0.5);
+      Int lag1(end+0.5);
+      if (mode == "frequency") {
+        lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
+        lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
+      }
+      Int lstart =  max(0, lag0);
+      Int lend =  min(Int(lags.nelements()-1), lag1);
+      if (lstart == lend) {
+        lags[lstart] = Complex(0.0);
       } else {
-        for (int j=start; j <=end ;++j) {
+        if (lstart > lend) {
+          Int tmp = lend;
+          lend = lstart;
+          lstart = tmp;
+        }
+        for (int j=lstart; j <=lend ;++j) {
           lags[j] = Complex(0.0);
         }

branches/alma/src/STMath.h

@@ -146 +146 @@
 
   casa::CountedPtr<Scantable>
-    binaryOperate( const casa::CountedPtr<Scantable>& left, 
-                   const casa::CountedPtr<Scantable>& right, 
+    binaryOperate( const casa::CountedPtr<Scantable>& left,
+                   const casa::CountedPtr<Scantable>& right,
                    const std::string& mode);
 
@@ -290 +290 @@
   casa::CountedPtr<Scantable>
     smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel,
-                      float width);
+                      float width, int order=2);
 
   casa::CountedPtr<Scantable>
@@ -302 +302 @@
 
   casa::CountedPtr<Scantable> opacity(const casa::CountedPtr<Scantable>& in,
-                                      float tau);
+                                      const std::vector<float>& tau);
 
   casa::CountedPtr<Scantable>
@@ -338 +338 @@
    */
   casa::CountedPtr<Scantable>
-    lagFlag( const casa::CountedPtr<Scantable>& in, double frequency,
-              double width);
+    lagFlag( const casa::CountedPtr<Scantable>& in, double start,
+             double end, const std::string& mode="frequency");
 
   // test for average spectra with different channel/resolution
@@ -374 +374 @@
                               bool tokelvin, float cfac);
 
-  casa::CountedPtr< Scantable > 
+  casa::CountedPtr< Scantable >
     smoothOther( const casa::CountedPtr< Scantable >& in,
                  const std::string& kernel,
-                 float width );
+                 float width, int order=2 );
 
   casa::CountedPtr< Scantable >

branches/alma/src/STMathWrapper.h

@@ -146 +146 @@
 
   ScantableWrapper
-    smooth(const ScantableWrapper& in, const std::string& kernel, float width)
-  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width)); }
+    smooth(const ScantableWrapper& in, const std::string& kernel, float width,
+           int order=2)
+  { return ScantableWrapper(STMath::smooth(in.getCP(), kernel, width, order)); }
 
   ScantableWrapper
@@ -164 +165 @@
 
   ScantableWrapper opacity(const ScantableWrapper& in,
-                                      float tau)
+                          const std::vector<float>& tau)
   { return ScantableWrapper(STMath::opacity(in.getCP(), tau)); }
 
@@ -201 +202 @@
 
   ScantableWrapper lagFlag( const ScantableWrapper& in,
-                            double frequency, double width )
-  { return ScantableWrapper(STMath::lagFlag(in.getCP(), frequency, width)); }
+                            double start, double end,
+                            const std::string& mode="frequency" )
+  { return ScantableWrapper(STMath::lagFlag(in.getCP(), start, end,
+                                            mode)); }
 
   // test for average spectra with different channel/resolution

branches/alma/src/STPol.h

@@ -35 +35 @@
 
   typedef  void (STPol::*polOperation)( casa::Float phase );
-  STPol(): totalfeed_(0.0),parangle_(0.0),feedhand_(1.0) {}
+  STPol(): totalangle_(0.0),feedhand_(1.0) {}
   virtual ~STPol() {}
 
@@ -83 +83 @@
 
 
-  void setPhaseCorrections(casa::Float parangle=0.0, casa::Float totalfeed=0.0,
-                           casa::Float feedhand=1.0)
-    { totalfeed_=totalfeed;parangle_=parangle;feedhand_=feedhand;}
+  void setPhaseCorrections(casa::Float totalang=0.0, casa::Float feedhand=1.0)
+    { totalangle_=totalang;feedhand_=feedhand;}
 
-  casa::Float getTotalPhase() const { return totalfeed_+parangle_; }
+  casa::Float getTotalPhase() const { return totalangle_; }
   casa::Float getFeedHand() const { return feedhand_; }
 
@@ -99 +98 @@
   static std::map<std::string, std::map<int, std::string> > labelmap_;
 
-  casa::Float totalfeed_,parangle_,feedhand_;
+  casa::Float totalangle_,feedhand_;
   std::string mode_;
   casa::Matrix<casa::Float> basespectra_;

branches/alma/src/STWriter.cpp

@@ -42 +42 @@
 #include <atnf/PKSIO/PKSMS2writer.h>
 #include <atnf/PKSIO/PKSSDwriter.h>
+#include <atnf/PKSIO/SrcType.h>
 
 #include <tables/Tables/Table.h>
@@ -51 +52 @@
 #include "STAsciiWriter.h"
 #include "STHeader.h"
+#include "STMath.h"
+
 
 #include "STWriter.h"
@@ -104 +107 @@
                     const std::string &filename)
 {
+  // If we write out foreign formats we have to convert the frequency system
+  // into the output frame, as we do everything related to SPectarlCoordinates
+  // in asap on-the-fly.
+
+  CountedPtr<Scantable> inst = in;
+  if (in->frequencies().getFrame(true) != in->frequencies().getFrame(false)) {
+    STMath stm(false);
+    inst = stm.frequencyAlign(in);
+  }
 
   if (format_=="ASCII") {
     STAsciiWriter iw;
+    // ASCII calls SpectralCoordinate::toWorld so no freqAlign use 'in'
     if (iw.write(*in, filename)) {
       return 0;
@@ -117 +130 @@
       iw.setClass(True);
     }
-    if (iw.write(*in, filename)) {
+    if (iw.write(*inst, filename)) {
       return 0;
     }
@@ -127 +140 @@
   // this is a little different from what I have done
   // before. Need to check with the Offline User Test data
-  STHeader hdr = in->getHeader();
+  STHeader hdr = inst->getHeader();
   //const Int nPol  = hdr.npol;
   //const Int nChan = hdr.nchan;
-  std::vector<uint> ifs = in->getIFNos();
-  int nIF = in->nif();//ifs.size();
+  std::vector<uint> ifs = inst->getIFNos();
+  int nIF = inst->nif();//ifs.size();
   Vector<uInt> nPol(nIF),nChan(nIF);
   Vector<Bool> havexpol(nIF);
   String fluxUnit = hdr.fluxunit;
-
+  fluxUnit.upcase();
   nPol = 0;nChan = 0; havexpol = False;
   for (uint i=0;i<ifs.size();++i) {
-    nPol(ifs[i]) = in->npol();
-    nChan(ifs[i]) = in->nchan(ifs[i]);
+    nPol(ifs[i]) = inst->npol();
+    nChan(ifs[i]) = inst->nchan(ifs[i]);
     havexpol(ifs[i]) = nPol(ifs[i]) > 2;
   }
-
-  const Table table = in->table();
+//   Vector<String> obstypes(2);
+//   obstypes(0) = "TR";//on
+//   obstypes(1) = "RF TR";//off
+  const Table table = inst->table();
+
 
   // Create the output file and write static data.
@@ -155 +171 @@
     throw(AipsError("Failed to create output file"));
   }
-
 
   Int count = 0;
@@ -206 +221 @@
           String tcalt;
           Vector<String> stmp0, stmp1;
-          in->frequencies().getEntry(crpix,crval, pksrec.freqInc,
+          inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
                                      rec.asuInt("FREQ_ID"));
-          in->focus().getEntry(pksrec.focusAxi, pksrec.focusTan,
-                               pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
-                               rec.asuInt("FOCUS_ID"));
-          in->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
+          inst->focus().getEntry(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
+                                 pksrec.focusRot, tmp0,tmp1,tmp2,tmp3,tmp4,
+                                 rec.asuInt("FOCUS_ID"));
+          inst->molecules().getEntry(pksrec.restFreq,stmp0,stmp1,
                                    rec.asuInt("MOLECULE_ID"));
-          in->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
+          inst->tcal().getEntry(pksrec.tcalTime, pksrec.tcal,
                               rec.asuInt("TCAL_ID"));
-          in->weather().getEntry(pksrec.temperature, pksrec.pressure,
+          inst->weather().getEntry(pksrec.temperature, pksrec.pressure,
                                  pksrec.humidity, pksrec.windSpeed,
                                  pksrec.windAz, rec.asuInt("WEATHER_ID"));
@@ -230 +245 @@
           pksrec.fieldName = rec.asString("FIELDNAME");
           pksrec.srcName   = rec.asString("SRCNAME");
-          pksrec.obsType   = hdr.obstype;
+          //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
+          pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
           pksrec.bandwidth = nchan * abs(pksrec.freqInc);
           pksrec.azimuth   = rec.asFloat("AZIMUTH");
           pksrec.elevation = rec.asFloat("ELEVATION");
-          pksrec.parAngle  = rec.asFloat("PARANGLE");
           pksrec.refBeam   = rec.asInt("REFBEAMNO") + 1;
           pksrec.sigma.resize(npol);
@@ -293 +308 @@
 {
   String poltype = tab.keywordSet().asString("POLTYPE");
-  if ( poltype == "stokes") {
+// Full stokes is not supported. Just allow stokes I
+  if ( poltype == "stokes" && tab.nrow() != 1) {
     String msg = "poltype = " + poltype + " not yet supported in output.";
     throw(AipsError(msg));
@@ -340 +356 @@
 }
 
-}
+// get obsType string from SRCTYPE value
+String STWriter::getObsTypes( Int srctype )
+{
+  String obsType ;
+  switch( srctype ) {
+  case Int(SrcType::PSON):
+    obsType = "PSON" ;
+    break ;
+  case Int(SrcType::PSOFF):
+    obsType = "PSOFF" ;
+    break ;
+  case Int(SrcType::NOD):
+    obsType = "NOD" ;
+    break ;
+  case Int(SrcType::FSON):
+    obsType = "FSON" ;
+    break ;
+  case Int(SrcType::FSOFF):
+    obsType = "FSOFF" ;
+    break ;
+  case Int(SrcType::SKY):
+    obsType = "SKY" ;
+    break ;
+  case Int(SrcType::HOT):
+    obsType = "HOT" ;
+    break ;
+  case Int(SrcType::WARM):
+    obsType = "WARM" ;
+    break ;
+  case Int(SrcType::COLD):
+    obsType = "COLD" ;
+    break ;
+  case Int(SrcType::PONCAL):
+    obsType = "PSON:CALON" ;
+    break ;
+  case Int(SrcType::POFFCAL):
+    obsType = "PSOFF:CALON" ;
+    break ;
+  case Int(SrcType::NODCAL):
+    obsType = "NOD:CALON" ;
+    break ;
+  case Int(SrcType::FONCAL):
+    obsType = "FSON:CALON" ;
+    break ;
+  case Int(SrcType::FOFFCAL):
+    obsType = "FSOFF:CALOFF" ;
+    break ;
+  case Int(SrcType::FSLO):
+    obsType = "FSLO" ;
+    break ;
+  case Int(SrcType::FLOOFF):
+    obsType = "FS:LOWER:OFF" ;
+    break ;
+  case Int(SrcType::FLOSKY):
+    obsType = "FS:LOWER:SKY" ;
+    break ;
+  case Int(SrcType::FLOHOT):
+    obsType = "FS:LOWER:HOT" ;
+    break ;
+  case Int(SrcType::FLOWARM):
+    obsType = "FS:LOWER:WARM" ;
+    break ;
+  case Int(SrcType::FLOCOLD):
+    obsType = "FS:LOWER:COLD" ;
+    break ;
+  case Int(SrcType::FSHI):
+    obsType = "FSHI" ;
+    break ;
+  case Int(SrcType::FHIOFF):
+    obsType = "FS:HIGHER:OFF" ;
+    break ;
+  case Int(SrcType::FHISKY):
+    obsType = "FS:HIGHER:SKY" ;
+    break ;
+  case Int(SrcType::FHIHOT):
+    obsType = "FS:HIGHER:HOT" ;
+    break ;
+  case Int(SrcType::FHIWARM):
+    obsType = "FS:HIGHER:WARM" ;
+    break ;
+  case Int(SrcType::FHICOLD):
+    obsType = "FS:HIGHER:COLD" ;
+    break ;
+  default:
+    obsType = "NOTYPE" ;
+  }
+
+  return obsType ;
+}
+
+}

branches/alma/src/STWriter.h

@@ -36 +36 @@
 #include <casa/aips.h>
 #include <casa/Utilities/CountedPtr.h>
+#include <casa/BasicSL/String.h>
 
 #include "Logger.h"
@@ -84 +85 @@
   void replacePtTab(const casa::Table& tab, const std::string& fname);
 
+  casa::String getObsTypes( casa::Int srctype ) ;
+
   std::string     format_;
   PKSwriter* writer_;

branches/alma/src/Scantable.cpp

@@ -261 +261 @@
   td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
   td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
-  td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
   td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
 
@@ -303 +302 @@
   elCol_.attach(table_, "ELEVATION");
   dirCol_.attach(table_, "DIRECTION");
-  paraCol_.attach(table_, "PARANGLE");
   fldnCol_.attach(table_, "FIELDNAME");
   rbeamCol_.attach(table_, "REFBEAMNO");
 
+  mweatheridCol_.attach(table_,"WEATHER_ID");
   mfitidCol_.attach(table_,"FIT_ID");
   mfreqidCol_.attach(table_, "FREQ_ID");
@@ -501 +500 @@
 }
 
-MPosition Scantable::getAntennaPosition () const
+MPosition Scantable::getAntennaPosition() const
 {
   Vector<Double> antpos;
@@ -515 +514 @@
   /// @todo reindex SCANNO, recompute nbeam, nif, npol
   inname = path.expandedName();
-  table_.deepCopy(inname, Table::New);
+  // WORKAROUND !!! for Table bug
+  // Remove when fixed in casacore
+  if ( table_.tableType() == Table::Memory  && !selector_.empty() ) {
+    Table tab = table_.copyToMemoryTable(generateName());
+    tab.deepCopy(inname, Table::New);
+    tab.markForDelete();
+
+  } else {
+    table_.deepCopy(inname, Table::New);
+  }
 }
 
@@ -646 +654 @@
 }
 
-std::vector<uint> Scantable::getNumbers(ScalarColumn<uInt>& col)
+std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
 {
   Vector<uInt> nos(col.getColumn());
@@ -840 +848 @@
     specCol_.get(whichrow, arr);
   } else {
-    CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_, basetype));
+    CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
+                                               basetype));
     uInt row = uInt(whichrow);
     stpol->setSpectra(getPolMatrix(row));
     Float fang,fhand,parang;
-    fang = focusTable_.getTotalFeedAngle(mfocusidCol_(row));
+    fang = focusTable_.getTotalAngle(mfocusidCol_(row));
     fhand = focusTable_.getFeedHand(mfocusidCol_(row));
-    parang = paraCol_(row);
-    /// @todo re-enable this
-    // disable total feed angle to support paralactifying Caswell style
-    stpol->setPhaseCorrections(parang, -parang, fhand);
+    stpol->setPhaseCorrections(fang, fhand);
     arr = stpol->getSpectrum(requestedpol, ptype);
   }
@@ -921 +927 @@
       << setw(15) << "Polarisations:" << setw(4) << npol()
       << "(" << getPolType() << ")" << endl
-      << setw(15) << "Channels:"  << setw(4) << nchan() << endl;
-  oss << endl;
+      << setw(15) << "Channels:" << nchan() << endl;
   String tmp;
   oss << setw(15) << "Observer:"
@@ -967 +972 @@
       << setw(10) << "Time" << setw(18) << "Integration" << endl;
   oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
-  oss << setw(10) << "" << setw(3) << "IF" << setw(6) << ""
+  oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
       << setw(8) << "Frame" << setw(16)
-      << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
+      << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
+      << setw(7) << "Channels"
+      << endl;
   oss << asap::SEPERATOR << endl;
   TableIterator iter(table_, "SCANNO");
@@ -1004 +1011 @@
         ROTableRow irow(isubt);
         const TableRecord& irec = irow.get(0);
-        oss << setw(10) << "";
+        oss << setw(9) << "";
         oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
-            << setw(2) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
+            << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
+            << setw(3) << "" << nchan(irec.asuInt("IFNO"))
             << endl;
 
@@ -1040 +1048 @@
     Double tm;
     table_.keywordSet().get("UTC",tm);
-    return MEpoch(MVEpoch(tm));  
+    return MEpoch(MVEpoch(tm));
   }
 }
@@ -1047 +1055 @@
 {
   return formatDirection(getDirection(uInt(whichrow)));
+}
+
+
+SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
+  const MPosition& mp = getAntennaPosition();
+  const MDirection& md = getDirection(whichrow);
+  const MEpoch& me = timeCol_(whichrow);
+  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  return freqTable_.getSpectralCoordinate(md, mp, me, rf,
+                                          mfreqidCol_(whichrow));
 }
 
@@ -1061 +1080 @@
     return stlout;
   }
-
-  const MPosition& mp = getAntennaPosition();
-  const MDirection& md = getDirection(whichrow);
-  const MEpoch& me = timeCol_(whichrow);
-  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
-  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
-  SpectralCoordinate spc =
-    freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
+  SpectralCoordinate spc = getSpectralCoordinate(whichrow);
   Vector<Double> pixel(nchan);
   Vector<Double> world;
@@ -1225 +1237 @@
            Sort::QuickSort|Sort::NoDuplicates );
   for (uInt i=0; i<fids.nelements(); ++i) {
-    frequencies().shiftRefPix(npix, i);
-  }
-}
-
-std::string asap::Scantable::getAntennaName() const
+    frequencies().shiftRefPix(npix, fids[i]);
+  }
+}
+
+std::string Scantable::getAntennaName() const
 {
   String out;
@@ -1236 +1248 @@
 }
 
-int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
+int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
 {
   String tbpath;
@@ -1321 +1333 @@
 }
 
-std::vector<double>  asap::Scantable::getDirectionVector(int whichrow) const
+std::vector<double> Scantable::getDirectionVector(int whichrow) const
 {
   Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
@@ -1681 +1693 @@
 }
 
+std::vector<float> Scantable::getWeather(int whichrow) const
+{
+  std::vector<float> out(5);
+  //Float temperature, pressure, humidity, windspeed, windaz;
+  weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
+                         mweatheridCol_(uInt(whichrow)));
+
+
+  return out;
+}
+
 }
 //namespace asap

branches/alma/src/Scantable.h

@@ -22 +22 @@
 #include <casa/Utilities/CountedPtr.h>
 
-#include <casa/Exceptions/Error.h>
-
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-
 #include <tables/Tables/Table.h>
 #include <tables/Tables/ArrayColumn.h>
@@ -32 +28 @@
 #include <measures/TableMeasures/ScalarMeasColumn.h>
 
+#include <coordinates/Coordinates/SpectralCoordinate.h>
+
 #include <casa/Arrays/Vector.h>
 #include <casa/Quanta/Quantum.h>
+
+#include <casa/Exceptions/Error.h>
 
 #include "Logger.h"
@@ -173 +173 @@
    */
   casa::MDirection getDirection( int whichrow ) const;
-  
+
   /**
    * get the direction type as a string, e.g. "J2000"
@@ -191 +191 @@
    * @return a string describing the direction reference
    */
-  std::string getDirectionRefString() const;    
+  std::string getDirectionRefString() const;
 
   /**
@@ -298 +298 @@
 
   int getBeam(int whichrow) const;
-  std::vector<uint> getBeamNos() { return getNumbers(beamCol_); }
+  std::vector<uint> getBeamNos() const { return getNumbers(beamCol_); }
 
   int getIF(int whichrow) const;
-  std::vector<uint> getIFNos() { return getNumbers(ifCol_); }
+  std::vector<uint> getIFNos() const { return getNumbers(ifCol_); }
 
   int getPol(int whichrow) const;
-  std::vector<uint> getPolNos() { return getNumbers(polCol_); }
-
-  std::vector<uint> getScanNos() { return getNumbers(scanCol_); }
+  std::vector<uint> getPolNos() const { return getNumbers(polCol_); }
+
+  std::vector<uint> getScanNos() const { return getNumbers(scanCol_); }
   int getScan(int whichrow) const { return scanCol_(whichrow); }
 
@@ -335 +335 @@
     { return azCol_(whichrow); }
   float getParAngle(int whichrow) const
-    { return paraCol_(whichrow); }
+    { return focus().getParAngle(mfocusidCol_(whichrow)); }
   int getTcalId(int whichrow) const
     { return mtcalidCol_(whichrow); }
@@ -384 +384 @@
 
   std::vector<double> getAbcissa(int whichrow) const;
+
+  std::vector<float> getWeather(int whichrow) const;
 
   std::string getAbcissaLabel(int whichrow) const;
@@ -408 +410 @@
 
   void shift(int npix);
+
+  casa::SpectralCoordinate getSpectralCoordinate(int whichrow) const;
 
   void convertDirection(const std::string& newframe);
@@ -443 +447 @@
    * against the information found in GBT_GO table for
    * scan number orders to get correct pairs.
-   * @param[in] scan list 
-   * @return status 
+   * @param[in] scan list
+   * @return status
    */
   int checkScanInfo(const std::vector<int>& scanlist) const;
 
   /**
-   * Get the direction as a vector, for a specific row  
+   * Get the direction as a vector, for a specific row
    * @param[in] whichrow the row numbyyer
-   * @return the direction in a vector 
+   * @return the direction in a vector
    */
   std::vector<double> getDirectionVector(int whichrow) const;
+
+  /**
+   * Set a flag indicating whether the data was parallactified
+   * @param[in] flag true or false
+   */
+  void parallactify(bool flag)
+    { focus().setParallactify(flag); }
 
   /**
@@ -523 +534 @@
   int rowToScanIndex(int therow);
 
-  std::vector<uint> getNumbers(casa::ScalarColumn<casa::uInt>& col);
-
-  static const casa::uInt version_ = 2;
+  std::vector<uint> getNumbers(const casa::ScalarColumn<casa::uInt>& col) const;
+
+  static const casa::uInt version_ = 3;
 
   STSelector selector_;
@@ -549 +560 @@
   casa::ScalarColumn<casa::Float> azCol_;
   casa::ScalarColumn<casa::Float> elCol_;
-  casa::ScalarColumn<casa::Float> paraCol_;
   casa::ScalarColumn<casa::String> srcnCol_, fldnCol_;
   casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_;

branches/alma/src/ScantableWrapper.h

@@ -20 +20 @@
 #include "STFit.h"
 #include "Scantable.h"
+#include "STCoordinate.h"
 
 namespace asap {
@@ -233 +234 @@
     { return table_->checkScanInfo(scanlist); }
  
-  std::vector<double>  getDirectionVector(int whichrow) const
+  std::vector<double> getDirectionVector(int whichrow) const
     { return table_->getDirectionVector(whichrow); }
+
+  void parallactify(bool flag)
+    { table_->parallactify(flag); }
+
+  STCoordinate getCoordinate(int row) {
+    return STCoordinate(table_->getSpectralCoordinate(row));
+  }
+
+  std::vector<float> getWeather(int whichrow) const
+    { return table_->getWeather(whichrow); }
 
   void reshapeSpectrum( int nmin, int nmax )

branches/alma/src/python_Scantable.cpp

@@ -133 +133 @@
     .def("_setsourcetype", &ScantableWrapper::setSourceType)
     .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
+    .def("_parallactify", &ScantableWrapper::parallactify)
+    .def("get_coordinate", &ScantableWrapper::getCoordinate)
+    .def("_get_weather", &ScantableWrapper::getWeather)
     .def("_reshape", &ScantableWrapper::reshapeSpectrum, 
          (boost::python::arg("nmin")=-1, 

branches/alma/src/python_asap.cpp

@@ -75 +75 @@
   asap::python::python_LineCatalog();
   asap::python::python_Logger();
+  asap::python::python_STCoordinate();
+  asap::python::python_STAtmosphere();
   asap::python::python_SrcType();
 

branches/alma/src/python_asap.h

@@ -47 +47 @@
     void python_LineCatalog();
     void python_Logger();
+    void python_STCoordinate();
+    void python_STAtmosphere();
     void python_SrcType();