Changeset 805

Timestamp:

02/16/06 12:02:18 (18 years ago)

Author:

mar637

Message:

Code replacemnts after the rename

Location:

trunk/src

Files:

• STFiller.cpp (modified) (12 diffs)
• STFiller.h (modified) (3 diffs)
• STMath.cpp (modified) (3 diffs)
• STMath.h (modified) (1 diff)
• Scantable.cpp (modified) (17 diffs)

trunk/src/STFiller.cpp

@@ -1 @@
-//#---------------------------------------------------------------------------
-//# SDReader.cc: A class to read single dish spectra from SDFITS, RPFITS
-//#---------------------------------------------------------------------------
-//# Copyright (C) 2004
-//# ATNF
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the Free
-//# Software Foundation; either version 2 of the License, or (at your option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
-//# WITHOUT ANY WARRANTY; without even the implied warranty of
-//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-//# Public License for more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning this software should be addressed as follows:
-//#        Internet email: Malte.Marquarding@csiro.au
-//#        Postal address: Malte Marquarding,
-//#                        Australia Telescope National Facility,
-//#                        P.O. Box 76,
-//#                        Epping, NSW, 2121,
-//#                        AUSTRALIA
-//#
-//# $Id:
-//#---------------------------------------------------------------------------
+//
+// C++ Implementation: STFiller
+//
+// Description:
+//
+//
+// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
 #include <casa/iostream.h>
 #include <casa/iomanip.h>
@@ -36 +17 @@
 #include <casa/OS/File.h>
 #include <casa/Quanta/Unit.h>
+#include <casa/Arrays/ArrayMath.h>
+#include <casa/Utilities/Regex.h>
+
+#include <casa/Containers/RecordField.h>
+
+#include <tables/Tables/TableRow.h>
+
 #include <atnf/PKSIO/PKSreader.h>
 
-#include "SDReader.h"
 #include "SDDefs.h"
 #include "SDAttr.h"
 
+#include "STFiller.h"
+
 using namespace casa;
-using namespace asap;
-
-SDReader::SDReader() :
+
+namespace asap {
+
+STFiller::STFiller() :
   reader_(0),
   header_(0),
-  frequencies_(0),
-  table_(new SDMemTable()),
-  haveXPol_(False)
-{
-  cursor_ = 0;
-}
-SDReader::SDReader(const std::string& filename,
-                   int whichIF, int whichBeam) :
+  table_(0)
+{
+}
+
+STFiller::STFiller( CountedPtr< Scantable > stbl ) :
   reader_(0),
   header_(0),
-  frequencies_(0),
-  table_(new SDMemTable()),
-  haveXPol_(False)
-{
-  cursor_ = 0;
-  open(filename, whichIF, whichBeam);
-}
-
-SDReader::SDReader(CountedPtr<SDMemTable> tbl) :
+  table_(stbl)
+{
+}
+
+STFiller::STFiller(const std::string & filename,
+                                 int whichIF, int whichBeam ) :
   reader_(0),
   header_(0),
-  table_(tbl),
-  haveXPol_(False)
-{
-  cursor_ = 0;
-}
-
-SDReader::~SDReader() {
-  if (reader_) delete reader_;
-  if (header_) delete header_;
-  if (frequencies_) delete frequencies_;
-}
-
-void SDReader::reset() {
-  cursor_ = 0;
-  table_ = new SDMemTable();
-  open(filename_,ifOffset_, beamOffset_);
-}
-
-
-void SDReader::close() {
-  //cerr << "disabled" << endl;
-}
-
-void SDReader::open(const std::string& filename,
-                    int whichIF, int whichBeam) {
-  if (reader_) delete reader_; reader_ = 0;
-  Bool   haveBase, haveSpectra;
+  table_(0)
+{
+  open(filename, whichIF, whichBeam);
+}
+
+STFiller::~STFiller()
+{
+  close();
+}
+
+void STFiller::open( const std::string & filename,
+                            int whichIF, int whichBeam )
+{
+  if (table_.null()) table_ = new Scantable();
+  if (reader_)  { delete reader_; reader_ = 0; }
+  Bool haveBase, haveSpectra;
 
   String inName(filename);
@@ -102 +75 @@
 
   File file(inName);
-  if (!file.exists()) {
+  if ( !file.exists() ) {
      throw(AipsError("File does not exist"));
   }
@@ -110 +83 @@
   String format;
   Vector<Bool> beams;
-  if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
+  if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, nIF_,
                               nChan_, nPol_, haveBase, haveSpectra,
-                              haveXPol_)) == 0)  {
+                              haveXPol_)) == 0 )  {
     throw(AipsError("Creation of PKSreader failed"));
   }
@@ -124 +97 @@
   nBeam_ = beams.nelements();
   // Get basic parameters.
-  if (haveXPol_) {
+  if ( haveXPol_ ) {
     pushLog("Cross polarization present");
     nPol_ += 2;                          // Convert Complex -> 2 Floats
@@ -145 +118 @@
     delete reader_;
     reader_ = 0;
+    delete header_;
+    header_ = 0;
     throw(AipsError("Failed to get header."));
-    return;
   }
   if ((header_->obstype).matches("*SW*")) {
@@ -153 +127 @@
                "setting # of IFs = 1 ");
     nIF_ = 1;
+    header_->obstype = String("fswitch");
   }
 
@@ -166 +141 @@
   header_->nif = nIF_;
   header_->epoch = "UTC";
-
+  // *** header_->frequnit = "Hz"
   // Apply selection criteria.
 
@@ -176 +151 @@
   if (whichIF>=0) {
     if (whichIF>=0 && whichIF<nIF_) {
-       IFsel = False;
-       IFsel(whichIF) = True;
-       header_->nif = 1;
-       nIF_ = 1;
-       ifOffset_ = whichIF;
+      IFsel = False;
+      IFsel(whichIF) = True;
+      header_->nif = 1;
+      nIF_ = 1;
+      ifOffset_ = whichIF;
     } else {
-       throw(AipsError("Illegal IF selection"));
+      delete reader_;
+      reader_ = 0;
+      delete header_;
+      header_ = 0;
+      throw(AipsError("Illegal IF selection"));
     }
   }
@@ -188 +167 @@
   beamOffset_ = 0;
   if (whichBeam>=0) {
-     if (whichBeam>=0 && whichBeam<nBeam_) {
-        beamSel = False;
-        beamSel(whichBeam) = True;
-        header_->nbeam = 1;
-        nBeam_ = 1;
-        beamOffset_ = whichBeam;
-     } else {
-       throw(AipsError("Illegal Beam selection"));
-     }
+    if (whichBeam>=0 && whichBeam<nBeam_) {
+      beamSel = False;
+      beamSel(whichBeam) = True;
+      header_->nbeam = 1;
+      nBeam_ = 1;
+      beamOffset_ = whichBeam;
+    } else {
+      delete reader_;
+      reader_ = 0;
+      delete header_;
+      header_ = 0;
+      throw(AipsError("Illegal Beam selection"));
+    }
   }
   Vector<Int> start(nIF_, 1);
@@ -202 +185 @@
   reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol_);
   table_->putSDHeader(*header_);
-
-  if (frequencies_) delete frequencies_;
-  frequencies_ = new SDFrequencyTable();
-  frequencies_->setRefFrame(header_->freqref);
-  frequencies_->setBaseRefFrame(header_->freqref);
-  frequencies_->setRestFrequencyUnit("Hz");
-  frequencies_->setEquinox(header_->equinox);
-}
-
-int SDReader::read(const std::vector<int>& seq) {
+}
+
+void STFiller::close( )
+{
+  delete reader_;reader_=0;
+  delete header_;header_=0;
+  table_ = 0;
+}
+
+int asap::STFiller::read( )
+{
   int status = 0;
 
@@ -226 +210 @@
   Complex         xCalFctr;
   Vector<Complex> xPol;
-  uInt n = seq.size();
-
-
-  uInt stepsize = header_->nif*header_->nbeam;
-  uInt seqi = 0;
-  Bool getAll = False;
-
-  if (seq[n-1] == -1) getAll = True;
-  while ( (cursor_ <= seq[n-1]) || getAll) {
-    // only needs to be create if in seq
-    SDContainer sc(header_->nbeam,header_->nif,header_->npol,header_->nchan);
-    // iterate over one correlator integration cycle = nBeam*nIF
-    for (uInt row=0; row < stepsize; row++) {
-      // stepsize as well
-      // spectra(nChan,nPol)!!!
-      status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
-                             srcName, srcDir, srcPM, srcVel, IFno, refFreq,
-                             bandwidth, freqInc, restFreq, tcal, tcalTime,
-                             azimuth, elevation, parAngle, focusAxi,
-                             focusTan, focusRot, temperature, pressure,
-                             humidity, windSpeed, windAz, refBeam,
-                             beamNo, direction, scanRate,
-                             tsys, sigma, calFctr, baseLin, baseSub,
-                             spectra, flagtra, xCalFctr, xPol);
-
-      // Make sure beam/IF numbers are 0-relative - dealing with
-      // possible IF or Beam selection
-      beamNo = beamNo - beamOffset_ - 1;
-      IFno = IFno - ifOffset_ - 1;
-
-      if (status) {
-        if (status == -1) {
-          // EOF.
-          if (row > 0 && row < stepsize-1)
-            pushLog("incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan are present.");
-
-          // flush frequency table
-          table_->putSDFreqTable(*frequencies_);
-          return status;
-        }
-      }
-      // if in the given list
-      if (cursor_ == seq[seqi] || getAll) {
-        // add integration cycle
-        if (row==0) {
-          //add invariant info: scanNo, mjd, interval, fieldName,
-          //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
-          //focusRot, temperature, pressure, humidity, windSpeed,
-          //windAz  srcDir, srcPM, srcVel
-          sc.timestamp = mjd;
-          sc.interval = interval;
-          sc.sourcename = srcName;
-          sc.fieldname = fieldName;
-          sc.azimuth = azimuth;
-          sc.elevation = elevation;
-        }
-        // add specific info
-        // refPix = nChan/2+1 in  1-rel Integer arith.!
-        Int refPix = header_->nchan/2;       // 0-rel
-        uInt freqID = frequencies_->addFrequency(refPix, refFreq, freqInc);
-        uInt restFreqID = frequencies_->addRestFrequency(restFreq);
-
-        sc.setFrequencyMap(freqID, IFno);
-        sc.setRestFrequencyMap(restFreqID, IFno);
-
-        sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1];
-        sc.tcaltime = tcalTime;
-        sc.parangle = parAngle;
-        sc.refbeam = -1; //nbeams == 1
-        if (nBeam_ > 1) // circumvent a bug "asap0000" in read which
-                        // returns a random refbema number on multiple
-                        // reads
-          sc.refbeam = refBeam-1;//make it 0-based;
-        sc.scanid = scanNo-1;//make it 0-based
-        if (haveXPol_) {
-           sc.setSpectrum(spectra, xPol, beamNo, IFno);
-           sc.setFlags(flagtra,  beamNo, IFno, True);
-        } else {
-           sc.setSpectrum(spectra, beamNo, IFno);
-           sc.setFlags(flagtra,  beamNo, IFno, False);
-        }
-        sc.setTsys(tsys, beamNo, IFno, haveXPol_);
-        sc.setDirection(direction, beamNo);
-      }
-    }
-
-    if (cursor_ == seq[seqi] || getAll) {
-      // insert container into table/list
-      table_->putSDContainer(sc);
-      seqi++;// next in list
-    }
-    cursor_++;// increment position in file
-
-  }
-  table_->putSDFreqTable(*frequencies_);
+  while ( status == 0 ) {
+    status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
+                          srcName, srcDir, srcPM, srcVel, IFno, refFreq,
+                          bandwidth, freqInc, restFreq, tcal, tcalTime,
+                          azimuth, elevation, parAngle, focusAxi,
+                          focusTan, focusRot, temperature, pressure,
+                          humidity, windSpeed, windAz, refBeam,
+                          beamNo, direction, scanRate,
+                          tsys, sigma, calFctr, baseLin, baseSub,
+                          spectra, flagtra, xCalFctr, xPol);
+    if ( status != 0 ) break;
+    TableRow row(table_->table());
+    TableRecord& rec = row.record();
+    RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
+    *scanoCol = scanNo-1;
+    RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
+    *cyclenoCol = cycleNo-1;
+    RecordFieldPtr<Double> mjdCol(rec, "TIME");
+    *mjdCol = mjd;
+    RecordFieldPtr<Double> intCol(rec, "INTERVAL");
+    *intCol = interval;
+    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
+    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
+    // try to auto-identify if it is on or off.
+    Regex rx(".*_[e|w|R]$");
+    Regex rx2("_[e|w|R|S]$");
+    Int match = srcName.matches(rx);
+    if (match) {
+      *srcnCol = srcName;
+    } else {
+      *srcnCol = srcName.before(rx2);
+    }
+    //*srcnCol = srcName;//.before(rx2);
+    *srctCol = match;
+    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
+    *beamCol = beamNo-beamOffset_-1;
+    RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
+    Int rb = -1;
+    if (nBeam_ > 1 ) rb = refBeam-1;
+    *rbCol = rb;
+    RecordFieldPtr<uInt> ifCol(rec, "IFNO");
+    *ifCol = IFno-ifOffset_- 1;
+    uInt id;
+    /// @todo this has to change when nchan isn't global anymore
+    id = table_->frequencies().addEntry(Double(header_->nchan/2),
+                                            refFreq, freqInc);
+    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
+    *mfreqidCol = id;
+
+    id = table_->molecules().addEntry(restFreq);
+    RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
+    *molidCol = id;
+
+    id = table_->tcal().addEntry(tcalTime, tcal);
+    RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
+    *mcalidCol = id;
+    id = table_->weather().addEntry(temperature, pressure, humidity,
+                                    windSpeed, windAz);
+    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
+    *mweatheridCol = id;
+    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
+    id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
+    *mfocusidCol = id;
+    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
+    *dirCol = direction;
+    RecordFieldPtr<Double> azCol(rec, "AZIMUTH");
+    *azCol = azimuth;
+    RecordFieldPtr<Double> elCol(rec, "ELEVATION");
+    *elCol = elevation;
+
+    RecordFieldPtr<Float> parCol(rec, "PARANGLE");
+    *parCol = parAngle;
+
+    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
+    RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
+    RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
+
+    RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
+    // Turn the (nchan,npol) matrix and possible complex xPol vector
+    // into 2-4 rows in the scantable
+    Vector<Float> tsysvec(1);
+    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
+    uInt npol = (spectra.ncolumn()==1 ? 1: 2);
+    for ( uInt i=0; i< npol; ++i ) {
+      tsysvec = tsys(i);
+      *tsysCol = tsysvec;
+      *polnoCol = i;
+
+      *specCol = spectra.column(i);
+      *flagCol = flagtra.column(i);
+      table_->table().addRow();
+      row.put(table_->table().nrow()-1, rec);
+    }
+    if ( haveXPol_ ) {
+      // no tsys given for xpol, so emulate it
+      tsysvec = sqrt(tsys[0]*tsys[1]);
+      *tsysCol = tsysvec;
+      // add real part of cross pol
+      *polnoCol = 2;
+      Vector<Float> r(real(xPol));
+      *specCol = r;
+      // make up flags from linears
+      /// @fixme this has to be a bitwise or of both pols
+      *flagCol = flagtra.column(0);// | flagtra.column(1);
+      table_->table().addRow();
+      row.put(table_->table().nrow()-1, rec);
+      // ad imaginary part of cross pol
+      *polnoCol = 3;
+      Vector<Float> im(imag(xPol));
+      *specCol = im;
+      table_->table().addRow();
+      row.put(table_->table().nrow()-1, rec);
+    }
+  }
+  if (status > 0) {
+    close();
+    throw(AipsError("Reading error occured, data possibly corrupted."));
+  }
+
   return status;
 }
+
+}//namespace asap

trunk/src/STFiller.h

@@ -1 @@
-//#---------------------------------------------------------------------------
-//# SDReader.h: A class to read single dish spectra from SDFITS, RPFITS
-//#---------------------------------------------------------------------------
-//# Copyright (C) 2004
-//# ATNF
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the Free
-//# Software Foundation; either version 2 of the License, or (at your option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
-//# WITHOUT ANY WARRANTY; without even the implied warranty of
-//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-//# Public License for more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning this software should be addressed as follows:
-//#        Internet email: Malte.Marquarding@csiro.au
-//#        Postal address: Malte Marquarding,
-//#                        Australia Telescope National Facility,
-//#                        P.O. Box 76,
-//#                        Epping, NSW, 2121,
-//#                        AUSTRALIA
-//#
-//# $Id:
-//#---------------------------------------------------------------------------
-#ifndef SDREADER_H
-#define SDREADER_H
+//
+// C++ Interface: STFiller
+//
+// Description:
+//
+//
+// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef STFILLER_H
+#define STFILLER_H
 
 #include <vector>
@@ -41 +22 @@
 #include <casa/Arrays/Vector.h>
 
-#include "SDMemTable.h"
+#include "Scantable.h"
 #include "SDContainer.h"
 #include "SDLog.h"
@@ -49 +30 @@
 namespace asap {
 
-class SDReader : public SDLog {
+/**
+This class fills a Scantable from external data formats using the PKSReader class.
+
+@author   Malte Marquarding
+@date     2006/01/16
+@version  2.0a
+*/
+class STFiller : public SDLog {
 public:
-  SDReader();
-  SDReader(const std::string& filename, 
-           int whichIF=-1, int whichBeam=-1);
-  SDReader(casa::CountedPtr<SDMemTable> tbl);
-  virtual ~SDReader();
 
-  void open(const std::string& filename,
-            int whichIF=-1,
-            int whichBeam=-1);
-  void close();
-  int read(const std::vector<int>& seq);
+  STFiller();
 
-  casa::CountedPtr<SDMemTable> getTable() const { return table_;}
-
-  void reset();
-
-  std::vector<int> pseudoHeader() const {
-    std::vector<int> v;
-    v.push_back(nBeam_);v.push_back(nIF_);
-    v.push_back(nPol_);v.push_back(nChan_);
-    return v;
-  }
+  STFiller(casa::CountedPtr< Scantable > stbl);
 
 
-protected:
-  
+  /**
+    * A constructor for a filler with associated input file
+    * @param filename the input file (rpf,sdfite or ms)
+    * @param whichIF read a specific IF only (default -1 means all IFs)
+    * @param whichBeam read a specific beam only (default -1 means all beams)
+    */
+  STFiller( const std::string& filename, int whichIF=-1,
+                  int whichBeam=-1 );
+
+  ~STFiller();
+
+  /**
+   * associate the Filler with a file on disk
+   * @param filename the input file (rpf,sdfite or ms)
+   * @param whichIF read a specific IF only (default -1 means all IFs)
+   * @param whichBeam read a specific beam only (default -1 means all beams)
+   * @exception AipsError Creation of PKSreader failed
+   */
+  void open( const std::string& filename, int whichIF=-1, int whichBeam=-1 );
+
+  /**
+   * detatch from file
+   */
+  void close( );
+
+  /**
+   * Read in "rows" from the source file attached with open()
+   * @return a status flag passed on by PKSreader
+   */
+  int read( );
+
+  casa::CountedPtr<Scantable> getTable() const { return table_;}
+
 private:
-  casa::Int nBeam_,nIF_,nPol_,nChan_;
-  PKSreader* reader_;  
+
+  PKSreader* reader_;
   SDHeader* header_;
-  SDFrequencyTable* frequencies_;
-  casa::CountedPtr<SDMemTable> table_;
   casa::String filename_;
-  casa::uInt cursor_;
-  casa::Double timestamp_;
-  casa::uInt beamOffset_, ifOffset_;
+  casa::CountedPtr< Scantable > table_;
+  casa::Int nIF_, nBeam_, nPol_, nChan_;
+  casa::uInt ifOffset_, beamOffset_;
   casa::Bool haveXPol_;
 };
 
-}// namespace
+} // namespace
+
 #endif

trunk/src/STMath.cpp

@@ -1 @@
-//#---------------------------------------------------------------------------
-//# SDMath.cc: A collection of single dish mathematical operations
-//#---------------------------------------------------------------------------
-//# Copyright (C) 2004
-//# ATNF
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the Free
-//# Software Foundation; either version 2 of the License, or (at your option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
-//# WITHOUT ANY WARRANTY; without even the implied warranty of
-//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-//# Public License for more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning this software should be addressed as follows:
-//#        Internet email: Malte.Marquarding@csiro.au
-//#        Postal address: Malte Marquarding,
-//#                        Australia Telescope National Facility,
-//#                        P.O. Box 76,
-//#                        Epping, NSW, 2121,
-//#                        AUSTRALIA
-//#
-//# $Id:
-//#---------------------------------------------------------------------------
-#include <vector>
-
-
-#include <casa/aips.h>
-#include <casa/iostream.h>
-#include <casa/sstream.h>
+//
+// C++ Implementation: STMath
+//
+// Description:
+//
+//
+// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
 #include <casa/iomanip.h>
+#include <casa/Exceptions/Error.h>
+#include <casa/Containers/Block.h>
 #include <casa/BasicSL/String.h>
-#include <casa/Arrays/IPosition.h>
-#include <casa/Arrays/Array.h>
-#include <casa/Arrays/ArrayIter.h>
-#include <casa/Arrays/VectorIter.h>
+#include <tables/Tables/TableIter.h>
+#include <tables/Tables/TableCopy.h>
+#include <casa/Arrays/MaskArrLogi.h>
+#include <casa/Arrays/MaskArrMath.h>
+#include <casa/Arrays/ArrayLogical.h>
 #include <casa/Arrays/ArrayMath.h>
-#include <casa/Arrays/ArrayLogical.h>
-#include <casa/Arrays/MaskedArray.h>
-#include <casa/Arrays/MaskArrMath.h>
-#include <casa/Arrays/MaskArrLogi.h>
-#include <casa/Arrays/Matrix.h>
-#include <casa/BasicMath/Math.h>
-#include <casa/Exceptions.h>
-#include <casa/Quanta/Quantum.h>
-#include <casa/Quanta/Unit.h>
-#include <casa/Quanta/MVEpoch.h>
-#include <casa/Quanta/MVTime.h>
-#include <casa/Utilities/Assert.h>
-
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-#include <coordinates/Coordinates/FrequencyAligner.h>
+#include <casa/Containers/RecordField.h>
+#include <tables/Tables/TableRow.h>
+#include <tables/Tables/TableVector.h>
+#include <tables/Tables/ExprNode.h>
+#include <tables/Tables/TableRecord.h>
+#include <tables/Tables/ReadAsciiTable.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>
 #include <scimath/Mathematics/Convolver.h>
-#include <scimath/Mathematics/InterpolateArray1D.h>
 #include <scimath/Functionals/Polynomial.h>
 
-#include <tables/Tables/Table.h>
-#include <tables/Tables/ScalarColumn.h>
-#include <tables/Tables/ArrayColumn.h>
-#include <tables/Tables/ReadAsciiTable.h>
-
 #include "MathUtils.h"
-#include "SDDefs.h"
+#include "RowAccumulator.h"
 #include "SDAttr.h"
-#include "SDContainer.h"
-#include "SDMemTable.h"
-
-#include "SDMath.h"
-#include "SDPol.h"
+#include "STMath.h"
 
 using namespace casa;
+
 using namespace asap;
 
-
-SDMath::SDMath()
-{
-}
-
-SDMath::SDMath(const SDMath& other)
-{
-
-// No state
-
-}
-
-SDMath& SDMath::operator=(const SDMath& other)
-{
-  if (this != &other) {
-// No state
-  }
-  return *this;
-}
-
-SDMath::~SDMath()
-{;}
-
-
-SDMemTable* SDMath::frequencyAlignment(const SDMemTable& in,
-                                       const String& refTime,
-                                       const String& method,
-                                       Bool perFreqID)
-{
-  // Get frame info from Table
-   std::vector<std::string> info = in.getCoordInfo();
-
-   // Parse frequency system
-   String systemStr(info[1]);
-   String baseSystemStr(info[3]);
-   if (baseSystemStr==systemStr) {
-     throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
-   }
-
-   MFrequency::Types freqSystem;
-   MFrequency::getType(freqSystem, systemStr);
-
-   return frequencyAlign(in, freqSystem, refTime, method, perFreqID);
-}
-
-
-
-CountedPtr<SDMemTable>
-SDMath::average(const std::vector<CountedPtr<SDMemTable> >& in,
-                const Vector<Bool>& mask, Bool scanAv,
-                const String& weightStr, Bool alignFreq)
-// Weighted averaging of spectra from one or more Tables.
-{
-  // Convert weight type
-  WeightType wtType = NONE;
-  convertWeightString(wtType, weightStr, True);
-
-  // Create output Table by cloning from the first table
-  SDMemTable* pTabOut = new SDMemTable(*in[0],True);
-  if (in.size() > 1) {
-    for (uInt i=1; i < in.size(); ++i) {
-      pTabOut->appendToHistoryTable(in[i]->getHistoryTable());
-    }
-  }
-  // Setup
-  IPosition shp = in[0]->rowAsMaskedArray(0).shape();      // Must not change
-  Array<Float> arr(shp);
-  Array<Bool> barr(shp);
-  const Bool useMask = (mask.nelements() == shp(asap::ChanAxis));
-
-  // Columns from Tables
-  ROArrayColumn<Float> tSysCol;
-  ROScalarColumn<Double> mjdCol;
-  ROScalarColumn<String> srcNameCol;
-  ROScalarColumn<Double> intCol;
-  ROArrayColumn<uInt> fqIDCol;
+STMath::STMath(bool insitu) :
+  insitu_(insitu)
+{
+}
+
+
+STMath::~STMath()
+{
+}
+
+CountedPtr<Scantable>
+STMath::average( const std::vector<CountedPtr<Scantable> >& in,
+                 const Vector<Bool>& mask,
+                 const std::string& weight,
+                 const std::string& avmode,
+                 bool alignfreq)
+{
+  if ( avmode == "SCAN" && in.size() != 1 )
+    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
+  WeightType wtype = stringToWeight(weight);
+  // output
+  // clone as this is non insitu
+  bool insitu = insitu_;
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(in[0], true);
+  setInsitu(insitu);
+
+  Table& tout = out->table();
+
+  /// @todo check if all scantables are conformant
+
+  ArrayColumn<Float> specColOut(tout,"SPECTRA");
+  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
+  ArrayColumn<Float> tsysColOut(tout,"TSYS");
+  ScalarColumn<Double> mjdColOut(tout,"TIME");
+  ScalarColumn<Double> intColOut(tout,"INTERVAL");
+
+  // set up the output table rows. These are based on the structure of the
+  // FIRST scantabel in the vector
+  const Table& baset = in[0]->table();
+
+  Block<String> cols(3);
+  cols[0] = String("BEAMNO");
+  cols[1] = String("IFNO");
+  cols[2] = String("POLNO");
+  if ( avmode == "SOURCE" ) {
+    cols.resize(4);
+    cols[3] = String("SRCNAME");
+  }
+  if ( avmode == "SCAN"  && in.size() == 1) {
+    cols.resize(4);
+    cols[3] = String("SCANNO");
+  }
+  uInt outrowCount = 0;
+  TableIterator iter(baset, cols);
+  while (!iter.pastEnd()) {
+    Table subt = iter.table();
+    // copy the first row of this selection into the new table
+    tout.addRow();
+    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
+    ++outrowCount;
+    ++iter;
+  }
+  RowAccumulator acc(wtype);
+  acc.setUserMask(mask);
+  ROTableRow row(tout);
+  ROArrayColumn<Float> specCol, tsysCol;
+  ROArrayColumn<uChar> flagCol;
+  ROScalarColumn<Double> mjdCol, intCol;
   ROScalarColumn<Int> scanIDCol;
 
-  // Create accumulation MaskedArray. We accumulate for each
-  // channel,if,pol,beam Note that the mask of the accumulation array
-  // will ALWAYS remain ALL True.  The MA is only used so that when
-  // data which is masked Bad is added to it, that data does not
-  // contribute.
-
-  Array<Float> zero(shp);
-  zero=0.0;
-  Array<Bool> good(shp);
-  good = True;
-  MaskedArray<Float> sum(zero,good);
-
-  // Counter arrays
-  Array<Float> nPts(shp);             // Number of points
-  nPts = 0.0;
-  Array<Float> nInc(shp);             // Increment
-  nInc = 1.0;
-
-  // Create accumulation Array for variance. We accumulate for each
-  // if,pol,beam, but average over channel.  So we need a shape with
-  // one less axis dropping channels.
-  const uInt nAxesSub = shp.nelements() - 1;
-  IPosition shp2(nAxesSub);
-  for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
-     if (i!=asap::ChanAxis) {
-       shp2(j) = shp(i);
-       j++;
-     }
-  }
-  Array<Float> sumSq(shp2);
-  sumSq = 0.0;
-  IPosition pos2(nAxesSub,0);                        // For indexing
-
-  // Time-related accumulators
-  Double time;
-  Double timeSum = 0.0;
-  Double intSum = 0.0;
-  Double interval = 0.0;
-
-  // To get the right shape for the Tsys accumulator we need to access
-  // a column from the first table.  The shape of this array must not
-  // change.  Note however that since the TSysSqSum array is used in a
-  // normalization process, and that I ignore the channel axis
-  // replication of values for now, it loses a dimension
-
-  Array<Float> tSysSum, tSysSqSum;
-  {
-    const Table& tabIn = in[0]->table();
-    tSysCol.attach(tabIn,"TSYS");
-    tSysSum.resize(tSysCol.shape(0));
-    tSysSqSum.resize(shp2);
-  }
-  tSysSum = 0.0;
-  tSysSqSum = 0.0;
-  Array<Float> tSys;
-
-  // Scan and row tracking
-  Int oldScanID = 0;
-  Int outScanID = 0;
-  Int scanID = 0;
-  Int rowStart = 0;
-  Int nAccum = 0;
-  Int tableStart = 0;
-
-  // Source and FreqID
-  String sourceName, oldSourceName, sourceNameStart;
-  Vector<uInt> freqID, freqIDStart, oldFreqID;
-
-  // Loop over tables
-  Float fac = 1.0;
-  const uInt nTables = in.size();
-  for (uInt iTab=0; iTab<nTables; iTab++) {
-
-    // Should check that the frequency tables don't change if doing
-    // FreqAlignment
-
-    // Attach columns to Table
-     const Table& tabIn = in[iTab]->table();
-     tSysCol.attach(tabIn, "TSYS");
-     mjdCol.attach(tabIn, "TIME");
-     srcNameCol.attach(tabIn, "SRCNAME");
-     intCol.attach(tabIn, "INTERVAL");
-     fqIDCol.attach(tabIn, "FREQID");
-     scanIDCol.attach(tabIn, "SCANID");
-
-     // Loop over rows in Table
-     const uInt nRows = in[iTab]->nRow();
-     for (uInt iRow=0; iRow<nRows; iRow++) {
-       // Check conformance
-        IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
-        if (!shp.isEqual(shp2)) {
-          delete pTabOut;
-           throw (AipsError("Shapes for all rows must be the same"));
+  for (uInt i=0; i < tout.nrow(); ++i) {
+    for ( int j=0; j < in.size(); ++j ) {
+      const Table& tin = in[j]->table();
+      const TableRecord& rec = row.get(i);
+      ROScalarColumn<Double> tmp(tin, "TIME");
+      Double td;tmp.get(0,td);
+      Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
+                       && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
+                       && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
+      Table subt;
+      if ( avmode == "SOURCE") {
+        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
+      } else if (avmode == "SCAN") {
+        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
+      } else {
+        subt = basesubt;
+      }
+      specCol.attach(subt,"SPECTRA");
+      flagCol.attach(subt,"FLAGTRA");
+      tsysCol.attach(subt,"TSYS");
+      intCol.attach(subt,"INTERVAL");
+      mjdCol.attach(subt,"TIME");
+      Vector<Float> spec,tsys;
+      Vector<uChar> flag;
+      Double inter,time;
+      for (uInt k = 0; k < subt.nrow(); ++k ) {
+        flagCol.get(k, flag);
+        Vector<Bool> bflag(flag.shape());
+        convertArray(bflag, flag);
+        if ( allEQ(bflag, True) ) {
+          continue;//don't accumulate
         }
-
-        // If we are not doing scan averages, make checks for source
-        // and frequency setup and warn if averaging across them
-        scanIDCol.getScalar(iRow, scanID);
-
-        // Get quantities from columns
-        srcNameCol.getScalar(iRow, sourceName);
-        mjdCol.get(iRow, time);
-        tSysCol.get(iRow, tSys);
-        intCol.get(iRow, interval);
-        fqIDCol.get(iRow, freqID);
-
-        // Initialize first source and freqID
-        if (iRow==0 && iTab==0) {
-          sourceNameStart = sourceName;
-          freqIDStart = freqID;
-        }
-
-        // If we are doing scan averages, see if we are at the end of
-        // an accumulation period (scan).  We must check soutce names
-        // too, since we might have two tables with one scan each but
-        // different source names; we shouldn't average different
-        // sources together
-        if (scanAv && ( (scanID != oldScanID)  ||
-                        (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
-
-          // Normalize data in 'sum' accumulation array according to
-          // weighting scheme
-           normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType,
-                     asap::ChanAxis, nAxesSub);
-
-           // Get ScanContainer for the first row of this averaged Scan
-           SDContainer scOut = in[iTab]->getSDContainer(rowStart);
-
-           // Fill scan container. The source and freqID come from the
-           // first row of the first table that went into this average
-           // ( should be the same for all rows in the scan average)
-
-           Float nR(nAccum);
-           fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
-                   timeSum/nR, intSum, sourceNameStart, freqIDStart);
-
-           // Write container out to Table
-           pTabOut->putSDContainer(scOut);
-
-           // Reset accumulators
-           sum = 0.0;
-           sumSq = 0.0;
-           nAccum = 0;
-
-           tSysSum =0.0;
-           tSysSqSum =0.0;
-           timeSum = 0.0;
-           intSum = 0.0;
-           nPts = 0.0;
-
-           // Increment
-           rowStart = iRow;              // First row for next accumulation
-           tableStart = iTab;            // First table for next accumulation
-           sourceNameStart = sourceName; // First source name for next
-                                         // accumulation
-           freqIDStart = freqID;         // First FreqID for next accumulation
-
-           oldScanID = scanID;
-           outScanID += 1;               // Scan ID for next
-                                         // accumulation period
-        }
-
-        // Accumulate
-        accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts,
-                   tSysSum, tSysSqSum, tSys,
-                   nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
-                   nAxesSub, useMask, wtType);
-        oldSourceName = sourceName;
-        oldFreqID = freqID;
-     }
-  }
-
-  // OK at this point we have accumulation data which is either
-  //   - accumulated from all tables into one row
-  // or
-  //   - accumulated from the last scan average
-  //
-  // Normalize data in 'sum' accumulation array according to weighting
-  // scheme
-
-  normalize(sum, sumSq, tSysSqSum, nPts, intSum, wtType,
-            asap::ChanAxis, nAxesSub);
-
-  // Create and fill container.  The container we clone will be from
-  // the last Table and the first row that went into the current
-  // accumulation.  It probably doesn't matter that much really...
-  Float nR(nAccum);
-  SDContainer scOut = in[tableStart]->getSDContainer(rowStart);
-  fillSDC(scOut, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
-          timeSum/nR, intSum, sourceNameStart, freqIDStart);
-  pTabOut->putSDContainer(scOut);
-  pTabOut->resetCursor();
-
-  return CountedPtr<SDMemTable>(pTabOut);
-}
-
-
-
-CountedPtr<SDMemTable>
-SDMath::binaryOperate(const CountedPtr<SDMemTable>& left,
-                      const CountedPtr<SDMemTable>& right,
-                      const String& op, Bool preserve, Bool doTSys)
-{
-
-  // Check operator
-  String op2(op);
-  op2.upcase();
-  uInt what = 0;
-  if (op2=="ADD") {
-     what = 0;
-  } else if (op2=="SUB") {
-     what = 1;
-  } else if (op2=="MUL") {
-     what = 2;
-  } else if (op2=="DIV") {
-     what = 3;
-  } else if (op2=="QUOTIENT") {
-     what = 4;
-     doTSys = True;
-  } else {
-    throw( AipsError("Unrecognized operation"));
-  }
-
-  // Check rows
-  const uInt nRowLeft = left->nRow();
-  const uInt nRowRight = right->nRow();
-  Bool ok = (nRowRight==1 && nRowLeft>0) ||
-            (nRowRight>=1 && nRowLeft==nRowRight);
-  if (!ok) {
-     throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
-  }
-
-  // Input Tables
-  const Table& tLeft = left->table();
-  const Table& tRight = right->table();
-
-  // TSys columns
-  ROArrayColumn<Float> tSysLeftCol, tSysRightCol;
-  if (doTSys) {
-    tSysLeftCol.attach(tLeft, "TSYS");
-    tSysRightCol.attach(tRight, "TSYS");
-  }
-
-  // First row for right
-  Array<Float> tSysLeftArr, tSysRightArr;
-  if (doTSys) tSysRightCol.get(0, tSysRightArr);
-  MaskedArray<Float>* pMRight =
-    new MaskedArray<Float>(right->rowAsMaskedArray(0));
-
-  IPosition shpRight = pMRight->shape();
-
-  // Output Table cloned from left
-  SDMemTable* pTabOut = new SDMemTable(*left, True);
-  pTabOut->appendToHistoryTable(right->getHistoryTable());
-
-  // Loop over rows
-  for (uInt i=0; i<nRowLeft; i++) {
-
-    // Get data
-    MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
-    IPosition shpLeft = mLeft.shape();
-    if (doTSys) tSysLeftCol.get(i, tSysLeftArr);
-
-    if (nRowRight>1) {
-      delete pMRight;
-      pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
-      shpRight = pMRight->shape();
-      if (doTSys) tSysRightCol.get(i, tSysRightArr);
-    }
-
-    if (!shpRight.isEqual(shpLeft)) {
-      delete pTabOut;
-      delete pMRight;
-      throw(AipsError("left and right scan tables are not conformant"));
-    }
-    if (doTSys) {
-      if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
-        delete pTabOut;
-        delete pMRight;
-        throw(AipsError("left and right Tsys data are not conformant"));
+        specCol.get(k, spec);
+        tsysCol.get(k, tsys);
+        intCol.get(k, inter);
+        mjdCol.get(k, time);
+        // spectrum has to be added last to enable weighting by the other values
+        acc.add(spec, !bflag, tsys, inter, time);
       }
-      if (!shpRight.isEqual(tSysRightArr.shape())) {
-        delete pTabOut;
-        delete pMRight;
-        throw(AipsError("left and right scan tables are not conformant"));
+    }
+    //write out
+    specColOut.put(i, acc.getSpectrum());
+    const Vector<Bool>& msk = acc.getMask();
+    Vector<uChar> flg(msk.shape());
+    convertArray(flg, !msk);
+    flagColOut.put(i, flg);
+    tsysColOut.put(i, acc.getTsys());
+    intColOut.put(i, acc.getInterval());
+    mjdColOut.put(i, acc.getTime());
+    acc.reset();
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
+                                             bool droprows)
+{
+  if (insitu_) return in;
+  else {
+    // clone
+    Scantable* tabp = new Scantable(*in, Bool(droprows));
+    return CountedPtr<Scantable>(tabp);
+  }
+}
+
+CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
+                                              float val,
+                                              const std::string& mode,
+                                              bool tsys )
+{
+  // modes are "ADD" and "MUL"
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tab = out->table();
+  ArrayColumn<Float> specCol(tab,"SPECTRA");
+  ArrayColumn<Float> tsysCol(tab,"TSYS");
+  for (uInt i=0; i<tab.nrow(); ++i) {
+    Vector<Float> spec;
+    Vector<Float> ts;
+    specCol.get(i, spec);
+    tsysCol.get(i, ts);
+    if (mode == "MUL") {
+      spec *= val;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts *= val;
+        tsysCol.put(i, ts);
       }
-    }
-
-    // Make container
-     SDContainer sc = left->getSDContainer(i);
-
-     // Operate on data and TSys
-     if (what==0) {
-        MaskedArray<Float> tmp = mLeft + *pMRight;
-        putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        if (doTSys) sc.putTsys(tSysLeftArr+tSysRightArr);
-     } else if (what==1) {
-        MaskedArray<Float> tmp = mLeft - *pMRight;
-        putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        if (doTSys) sc.putTsys(tSysLeftArr-tSysRightArr);
-     } else if (what==2) {
-        MaskedArray<Float> tmp = mLeft * *pMRight;
-        putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        if (doTSys) sc.putTsys(tSysLeftArr*tSysRightArr);
-     } else if (what==3) {
-        MaskedArray<Float> tmp = mLeft / *pMRight;
-        putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-        if (doTSys) sc.putTsys(tSysLeftArr/tSysRightArr);
-     } else if (what==4) {
-       if (preserve) {
-         MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
-           tSysRightArr;
-         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-       } else {
-         MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) -
-           tSysLeftArr;
-         putDataInSDC(sc, tmp.getArray(), tmp.getMask());
-       }
-       sc.putTsys(tSysRightArr);
-     }
-
-     // Put new row in output Table
-     pTabOut->putSDContainer(sc);
-  }
-  if (pMRight) delete pMRight;
-  pTabOut->resetCursor();
-
-  return CountedPtr<SDMemTable>(pTabOut);
-}
-
-
-std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
-                                     const Vector<Bool>& mask,
-                                     const String& which, Int row) const
-//
-// Perhaps iteration over pol/beam/if should be in here
-// and inside the nrow iteration ?
-//
-{
-  const uInt nRow = in->nRow();
-
-// Specify cursor location
-
-  IPosition start, end;
-  Bool doAll = False;
-  setCursorSlice(start, end, doAll, *in);
-
-// Loop over rows
-
-  const uInt nEl = mask.nelements();
-  uInt iStart = 0;
-  uInt iEnd = in->nRow()-1;
-//
-  if (row>=0) {
-     iStart = row;
-     iEnd = row;
-  }
-//
-  std::vector<float> result(iEnd-iStart+1);
-  for (uInt ii=iStart; ii <= iEnd; ++ii) {
-
-// Get row and deconstruct
-
-     MaskedArray<Float> dataIn = (in->rowAsMaskedArray(ii))(start,end);
-     Array<Float> v = dataIn.getArray().nonDegenerate();
-     Array<Bool>  m = dataIn.getMask().nonDegenerate();
-
-// Access desired piece of data
-
-//     Array<Float> v((arr(start,end)).nonDegenerate());
-//     Array<Bool> m((barr(start,end)).nonDegenerate());
-
-// Apply OTF mask
-
-     MaskedArray<Float> tmp;
-     if (m.nelements()==nEl) {
-       tmp.setData(v,m&&mask);
-     } else {
-       tmp.setData(v,m);
-     }
-
-// Get statistic
-
-     result[ii-iStart] = mathutil::statistics(which, tmp);
-  }
-//
-  return result;
-}
-
-
-SDMemTable* SDMath::bin(const SDMemTable& in, Int width)
-{
-  SDHeader sh = in.getSDHeader();
-  SDMemTable* pTabOut = new SDMemTable(in, True);
-
-// Bin up SpectralCoordinates
-
-  IPosition factors(1);
-  factors(0) = width;
-  for (uInt j=0; j<in.nCoordinates(); ++j) {
-    CoordinateSystem cSys;
-    cSys.addCoordinate(in.getSpectralCoordinate(j));
-    CoordinateSystem cSysBin =
-      CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
-//
-    SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
-    pTabOut->setCoordinate(sCBin, j);
-  }
-
-// Use RebinLattice to find shape
-
-  IPosition shapeIn(1,sh.nchan);
-  IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
-  sh.nchan = shapeOut(0);
-  pTabOut->putSDHeader(sh);
-
-// Loop over rows and bin along channel axis
-
-  for (uInt i=0; i < in.nRow(); ++i) {
-    SDContainer sc = in.getSDContainer(i);
-//
-    Array<Float> tSys(sc.getTsys());                           // Get it out before sc changes shape
-
-// Bin up spectrum
-
-    MaskedArray<Float> marr(in.rowAsMaskedArray(i));
-    MaskedArray<Float> marrout;
-    LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width);
-
-// Put back the binned data and flags
-
-    IPosition ip2 = marrout.shape();
-    sc.resize(ip2);
-//
-    putDataInSDC(sc, marrout.getArray(), marrout.getMask());
-
-// Bin up Tsys.
-
-    Array<Bool> allGood(tSys.shape(),True);
-    MaskedArray<Float> tSysIn(tSys, allGood, True);
-//
-    MaskedArray<Float> tSysOut;
-    LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width);
-    sc.putTsys(tSysOut.getArray());
-//
-    pTabOut->putSDContainer(sc);
-  }
-  return pTabOut;
-}
-
-SDMemTable* SDMath::resample(const SDMemTable& in, const String& methodStr,
-                             Float width)
+    } else if ( mode == "ADD" ) {
+      spec += val;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts += val;
+        tsysCol.put(i, ts);
+      }
+    }
+  }
+  return out;
+}
+
+MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
+                                        const Vector<uChar>& f)
+{
+  Vector<Bool> mask;
+  mask.resize(f.shape());
+  convertArray(mask, f);
+  return MaskedArray<Float>(s,!mask);
+}
+
+Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
+{
+  const Vector<Bool>& m = ma.getMask();
+  Vector<uChar> flags(m.shape());
+  convertArray(flags, !m);
+  return flags;
+}
+
+CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in,
+                                          const std::string & mode,
+                                          bool preserve )
+{
+  /// @todo make other modes available
+  /// modes should be "nearest", "pair"
+  // make this operation non insitu
+  const Table& tin = in->table();
+  Table ons = tin(tin.col("SRCTYPE") == Int(0));
+  Table offs = tin(tin.col("SRCTYPE") == Int(1));
+  if ( offs.nrow() == 0 )
+    throw(AipsError("No 'off' scans present."));
+  // put all "on" scans into output table
+
+  bool insitu = insitu_;
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(in, true);
+  setInsitu(insitu);
+  Table& tout = out->table();
+
+  TableCopy::copyRows(tout, ons);
+  TableRow row(tout);
+  ROScalarColumn<Double> offtimeCol(offs, "TIME");
+
+  ArrayColumn<Float> outspecCol(tout, "SPECTRA");
+  ROArrayColumn<Float> outtsysCol(tout, "TSYS");
+  ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
+
+  for (uInt i=0; i < tout.nrow(); ++i) {
+    const TableRecord& rec = row.get(i);
+    Double ontime = rec.asDouble("TIME");
+    ROScalarColumn<Double> offtimeCol(offs, "TIME");
+    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
+    Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat
+                       && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
+                       && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
+                       && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
+
+    TableRow offrow(sel);
+    const TableRecord& offrec = offrow.get(0);//should only be one row
+    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
+    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
+    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
+    /// @fixme this assumes tsys is a scalar not vector
+    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
+    Vector<Float> specon, tsyson;
+    outtsysCol.get(i, tsyson);
+    outspecCol.get(i, specon);
+    Vector<uChar> flagon;
+    outflagCol.get(i, flagon);
+    MaskedArray<Float> mon = maskedArray(specon, flagon);
+    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
+    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
+    if (preserve) {
+      quot -= tsysoffscalar;
+    } else {
+      quot -= tsyson[0];
+    }
+    outspecCol.put(i, quot.getArray());
+    outflagCol.put(i, flagsFromMA(quot));
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
+{
+  // make copy or reference
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tout = out->table();
+  Block<String> cols(3);
+  cols[0] = String("SCANNO");
+  cols[1] = String("BEAMNO");
+  cols[2] = String("POLNO");
+  TableIterator iter(tout, cols);
+  while (!iter.pastEnd()) {
+    Table subt = iter.table();
+    // this should leave us with two rows for the two IFs....if not ignore
+    if (subt.nrow() != 2 ) {
+      continue;
+    }
+    ArrayColumn<Float> specCol(tout, "SPECTRA");
+    ArrayColumn<Float> tsysCol(tout, "TSYS");
+    ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
+    Vector<Float> onspec,offspec, ontsys, offtsys;
+    Vector<uChar> onflag, offflag;
+    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
+    specCol.get(0, onspec);   specCol.get(1, offspec);
+    flagCol.get(0, onflag);   flagCol.get(1, offflag);
+    MaskedArray<Float> on  = maskedArray(onspec, onflag);
+    MaskedArray<Float> off = maskedArray(offspec, offflag);
+    MaskedArray<Float> oncopy = on.copy();
+
+    on /= off; on -= 1.0f;
+    on *= ontsys[0];
+    off /= oncopy; off -= 1.0f;
+    off *= offtsys[0];
+    specCol.put(0, on.getArray());
+    const Vector<Bool>& m0 = on.getMask();
+    Vector<uChar> flags0(m0.shape());
+    convertArray(flags0, !m0);
+    flagCol.put(0, flags0);
+
+    specCol.put(1, off.getArray());
+    const Vector<Bool>& m1 = off.getMask();
+    Vector<uChar> flags1(m1.shape());
+    convertArray(flags1, !m1);
+    flagCol.put(1, flags1);
+
+  }
+
+  return out;
+}
+
+std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
+                                        const std::vector< bool > & mask,
+                                        const std::string& which )
+{
+
+  Vector<Bool> m(mask);
+  const Table& tab = in->table();
+  ROArrayColumn<Float> specCol(tab, "SPECTRA");
+  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+  std::vector<float> out;
+  for (uInt i=0; i < tab.nrow(); ++i ) {
+    Vector<Float> spec; specCol.get(i, spec);
+    MaskedArray<Float> ma  = maskedArray(spec, flagCol(i));
+    float outstat;
+    if ( spec.nelements() == m.nelements() ) {
+      outstat = mathutil::statistics(which, ma(m));
+    } else {
+      outstat = mathutil::statistics(which, ma);
+    }
+    out.push_back(outstat);
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
+                                     int width )
+{
+  if ( !in->selection().empty() ) throw(AipsError("Can't bin subset of the data."));
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tout = out->table();
+  out->frequencies().rescale(width, "BIN");
+  ArrayColumn<Float> specCol(tout, "SPECTRA");
+  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
+  for (uInt i=0; i < tout.nrow(); ++i ) {
+    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
+    MaskedArray<Float> maout;
+    LatticeUtilities::bin(maout, main, 0, Int(width));
+    /// @todo implement channel based tsys binning
+    specCol.put(i, maout.getArray());
+    flagCol.put(i, flagsFromMA(maout));
+    // take only the first binned spectrum's length for the deprecated
+    // global header item nChan
+    if (i==0) tout.rwKeywordSet().define(String("nChan"),
+                                       Int(maout.getArray().nelements()));
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
+                                          const std::string& method,
+                                          float width )
 //
 // Should add the possibility of width being specified in km/s. This means
@@ -627 +402 @@
 //
 {
-   Bool doVel = False;
-   if (doVel) {
-      for (uInt j=0; j<in.nCoordinates(); ++j) {
-         SpectralCoordinate sC = in.getSpectralCoordinate(j);
-      }
-   }
-
-// Interpolation method
-
   InterpolateArray1D<Double,Float>::InterpolationMethod interp;
-  convertInterpString(interp, methodStr);
-  Int interpMethod(interp);
-
-// Make output table
-
-  SDMemTable* pTabOut = new SDMemTable(in, True);
+  Int interpMethod(stringToIMethod(method));
+
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tout = out->table();
 
 // Resample SpectralCoordinates (one per freqID)
-
-  const uInt nCoord = in.nCoordinates();
-  Vector<Float> offset(1,0.0);
-  Vector<Float> factors(1,1.0/width);
-  Vector<Int> newShape;
-  for (uInt j=0; j<in.nCoordinates(); ++j) {
-    CoordinateSystem cSys;
-    cSys.addCoordinate(in.getSpectralCoordinate(j));
-    CoordinateSystem cSys2 = cSys.subImage(offset, factors, newShape);
-    SpectralCoordinate sC = cSys2.spectralCoordinate(0);
-//
-    pTabOut->setCoordinate(sC, j);
-  }
-
-// Get header
-
-  SDHeader sh = in.getSDHeader();
-
-// Generate resampling vectors
-
-  const uInt nChanIn = sh.nchan;
-  Vector<Float> xIn(nChanIn);
-  indgen(xIn);
-//
-  Int fac =  Int(nChanIn/width);
-  Vector<Float> xOut(fac+10);          // 10 to be safe - resize later
-  uInt i = 0;
-  Float x = 0.0;
-  Bool more = True;
-  while (more) {
-    xOut(i) = x;
-//
-    i++;
-    x += width;
-    if (x>nChanIn-1) more = False;
-  }
-  const uInt nChanOut = i;
-  xOut.resize(nChanOut,True);
-//
-  IPosition shapeIn(in.rowAsMaskedArray(0).shape());
-  sh.nchan = nChanOut;
-  pTabOut->putSDHeader(sh);
-
-// Loop over rows and resample along channel axis
-
-  Array<Float> valuesOut;
-  Array<Bool> maskOut;
-  Array<Float> tSysOut;
-  Array<Bool> tSysMaskIn(shapeIn,True);
-  Array<Bool> tSysMaskOut;
-  for (uInt i=0; i < in.nRow(); ++i) {
-
-// Get container
-
-     SDContainer sc = in.getSDContainer(i);
-
-// Get data and Tsys
-
-     const Array<Float>& tSysIn = sc.getTsys();
-     const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i));
-     Array<Float> valuesIn = dataIn.getArray();
-     Array<Bool> maskIn = dataIn.getMask();
-
-// Interpolate data
-
-     InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut,
-                                                  xIn, valuesIn, maskIn,
-                                                  interpMethod, True, True);
-     sc.resize(valuesOut.shape());
-     putDataInSDC(sc, valuesOut, maskOut);
-
-// Interpolate TSys
-
-     InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut,
-                                                  xIn, tSysIn, tSysMaskIn,
-                                                  interpMethod, True, True);
-    sc.putTsys(tSysOut);
-
-// Put container in output
-
-    pTabOut->putSDContainer(sc);
-  }
-//
-  return pTabOut;
-}
-
-SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
-                                 uInt what, Bool doTSys)
-//
-// what = 0   Multiply
-//        1   Add
-{
-   SDMemTable* pOut = new SDMemTable(in,False);
-   const Table& tOut = pOut->table();
-   ArrayColumn<Float> specCol(tOut,"SPECTRA");
-   ArrayColumn<Float> tSysCol(tOut,"TSYS");
-   Array<Float> tSysArr;
-
-// Get data slice bounds
-
-   IPosition start, end;
-   setCursorSlice (start, end, doAll, in);
-//
-   for (uInt i=0; i<tOut.nrow(); i++) {
-
-// Modify data
-
-      MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
-      MaskedArray<Float> dataIn2 = dataIn(start,end);    // Reference
-      if (what==0) {
-         dataIn2 *= val;
-      } else if (what==1) {
-         dataIn2 += val;
-      }
-      specCol.put(i, dataIn.getArray());
-
-// Modify Tsys
-
-      if (doTSys) {
-         tSysCol.get(i, tSysArr);
-         Array<Float> tSysArr2 = tSysArr(start,end);     // Reference
-         if (what==0) {
-            tSysArr2 *= val;
-         } else if (what==1) {
-            tSysArr2 += val;
-         }
-         tSysCol.put(i, tSysArr);
-      }
-   }
-//
-   return pOut;
-}
-
-SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask,
-                               const String& weightStr)
-//
-// Average all polarizations together, weighted by variance
-//
-{
-   WeightType wtType = NONE;
-   convertWeightString(wtType, weightStr, True);
-
-// Create output Table and reshape number of polarizations
-
-  Bool clear=True;
-  SDMemTable* pTabOut = new SDMemTable(in, clear);
-  SDHeader header = pTabOut->getSDHeader();
-  header.npol = 1;
-  pTabOut->putSDHeader(header);
-//
-  const Table& tabIn = in.table();
-
-// Shape of input and output data
-
-  const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
-  IPosition shapeOut(shapeIn);
-  shapeOut(asap::PolAxis) = 1;                          // Average all polarizations
-  if (shapeIn(asap::PolAxis)==1) {
-    delete  pTabOut;
-    throw(AipsError("The input has only one polarisation"));
-  }
-//
-  const uInt nRows = in.nRow();
-  const uInt nChan = shapeIn(asap::ChanAxis);
-  AlwaysAssert(asap::nAxes==4,AipsError);
-  const IPosition vecShapeOut(4,1,1,1,nChan);     // A multi-dim form of a Vector shape
-  IPosition start(4), end(4);
-
-// Output arrays
-
-  Array<Float> outData(shapeOut, 0.0);
-  Array<Bool> outMask(shapeOut, True);
-  const IPosition axes(2, asap::PolAxis, asap::ChanAxis);              // pol-channel plane
-
-// Attach Tsys column if needed
-
-  ROArrayColumn<Float> tSysCol;
-  Array<Float> tSys;
-  if (wtType==TSYS) {
-     tSysCol.attach(tabIn,"TSYS");
-  }
-//
-  const Bool useMask = (mask.nelements() == shapeIn(asap::ChanAxis));
-
-// Loop over rows
-
-   for (uInt iRow=0; iRow<nRows; iRow++) {
-
-// Get data for this row
-
-      MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
-      Array<Float>& arr = marr.getRWArray();
-      const Array<Bool>& barr = marr.getMask();
-
-// Get Tsys
-
-      if (wtType==TSYS) {
-         tSysCol.get(iRow,tSys);
-      }
-
-// Make iterators to iterate by pol-channel planes
-// The tSys array is empty unless wtType=TSYS so only
-// access the iterator is that is the case
-
-      ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
-      ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
-      ReadOnlyArrayIterator<Float>* pItTsysPlane = 0;
-      if (wtType==TSYS)
-        pItTsysPlane = new ReadOnlyArrayIterator<Float>(tSys, axes);
-
-// Accumulations
-
-      Float fac = 1.0;
-      Vector<Float> vecSum(nChan,0.0);
-
-// Iterate through data by pol-channel planes
-
-      while (!itDataPlane.pastEnd()) {
-
-// Iterate through plane by polarization  and accumulate Vectors
-
-        Vector<Float> t1(nChan); t1 = 0.0;
-        Vector<Bool> t2(nChan); t2 = True;
-        Float tSys = 0.0;
-        MaskedArray<Float> vecSum(t1,t2);
-        Float norm = 0.0;
-        {
-           ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
-           ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
-//
-           ReadOnlyVectorIterator<Float>* pItTsysVec = 0;
-           if (wtType==TSYS) {
-              pItTsysVec =
-                new ReadOnlyVectorIterator<Float>(pItTsysPlane->array(), 1);
-           }
-//
-           while (!itDataVec.pastEnd()) {
-
-// Create MA of data & mask (optionally including OTF mask) and  get variance for this spectrum
-
-              if (useMask) {
-                 const MaskedArray<Float> spec(itDataVec.vector(),
-                                               mask&&itMaskVec.vector());
-                 if (wtType==VAR) {
-                    fac = 1.0 / variance(spec);
-                 } else if (wtType==TSYS) {
-                    tSys = pItTsysVec->vector()[0];      // Drop pseudo channel dependency
-                    fac = 1.0 / tSys / tSys;
-                 }
-              } else {
-                 const MaskedArray<Float> spec(itDataVec.vector(),
-                                               itMaskVec.vector());
-                 if (wtType==VAR) {
-                    fac = 1.0 / variance(spec);
-                 } else if (wtType==TSYS) {
-                    tSys = pItTsysVec->vector()[0];      // Drop pseudo channel dependency
-                    fac = 1.0 / tSys / tSys;
-                 }
-              }
-
-// Normalize spectrum (without OTF mask) and accumulate
-
-              const MaskedArray<Float> spec(fac*itDataVec.vector(),
-                                            itMaskVec.vector());
-              vecSum += spec;
-              norm += fac;
-
-// Next
-
-              itDataVec.next();
-              itMaskVec.next();
-              if (wtType==TSYS) pItTsysVec->next();
-           }
-
-// Clean up
-
-           if (pItTsysVec) {
-              delete pItTsysVec;
-              pItTsysVec = 0;
-           }
-        }
-
-// Normalize summed spectrum
-
-        vecSum /= norm;
-
-// FInd position in input data array.  We are iterating by pol-channel
-// plane so all that will change is beam and IF and that's what we want.
-
-        IPosition pos = itDataPlane.pos();
-
-// Write out data. This is a bit messy. We have to reform the Vector
-// accumulator into an Array of shape (1,1,1,nChan)
-
-        start = pos;
-        end = pos;
-        end(asap::ChanAxis) = nChan-1;
-        outData(start,end) = vecSum.getArray().reform(vecShapeOut);
-        outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
-
-// Step to next beam/IF combination
-
-        itDataPlane.next();
-        itMaskPlane.next();
-        if (wtType==TSYS) pItTsysPlane->next();
-      }
-
-// Generate output container and write it to output table
-
-      SDContainer sc = in.getSDContainer();
-      sc.resize(shapeOut);
-//
-      putDataInSDC(sc, outData, outMask);
-      pTabOut->putSDContainer(sc);
-//
-      if (wtType==TSYS) {
-         delete pItTsysPlane;
-         pItTsysPlane = 0;
-      }
-   }
-
-// Set polarization cursor to 0
-
-  pTabOut->setPol(0);
-//
-  return pTabOut;
-}
-
-
-SDMemTable* SDMath::smooth(const SDMemTable& in,
-                           const casa::String& kernelType,
-                           casa::Float width, Bool doAll)
-//
-// Should smooth TSys as well
-//
-{
-
-  // Number of channels
-   const uInt nChan = in.nChan();
-
-   // Generate Kernel
-   VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
-   Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
-
-   // Generate Convolver
-   IPosition shape(1,nChan);
-   Convolver<Float> conv(kernel, shape);
-
-   // New Table
-   SDMemTable* pTabOut = new SDMemTable(in,True);
-
-   // Output Vectors
-   Vector<Float> valuesOut(nChan);
-   Vector<Bool> maskOut(nChan);
-
-   // Get data slice bounds
-   IPosition start, end;
-   setCursorSlice (start, end, doAll, in);
-
-   // Loop over rows in Table
-   for (uInt ri=0; ri < in.nRow(); ++ri) {
-
-     // Get slice of data
-      MaskedArray<Float> dataIn = in.rowAsMaskedArray(ri);
-
-      // Deconstruct and get slices which reference these arrays
-      Array<Float> valuesIn = dataIn.getArray();
-      Array<Bool> maskIn = dataIn.getMask();
-
-      Array<Float> valuesIn2 = valuesIn(start,end);       // ref to valuesIn
-      Array<Bool> maskIn2 = maskIn(start,end);
-
-      // Iterate through by spectra
-      VectorIterator<Float> itValues(valuesIn2, asap::ChanAxis);
-      VectorIterator<Bool> itMask(maskIn2, asap::ChanAxis);
-      while (!itValues.pastEnd()) {
-
-        // Smooth
-        if (kernelType==VectorKernel::HANNING) {
-          mathutil::hanning(valuesOut, maskOut, itValues.vector(),
-                            itMask.vector());
-          itMask.vector() = maskOut;
-        } else {
-          mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
-          conv.linearConv(valuesOut, itValues.vector());
-        }
-
-        itValues.vector() = valuesOut;
-        itValues.next();
-        itMask.next();
-      }
-
-      // Create and put back
-      SDContainer sc = in.getSDContainer(ri);
-      putDataInSDC(sc, valuesIn, maskIn);
-
-      pTabOut->putSDContainer(sc);
-   }
-
-  return pTabOut;
-}
-
-
-
-SDMemTable* SDMath::convertFlux(const SDMemTable& in, Float D, Float etaAp,
-                                Float JyPerK, Bool doAll)
-//
-// etaAp = aperture efficiency (-1 means find)
-// D     = geometric diameter (m)  (-1 means find)
-// JyPerK
-//
-{
-  SDHeader sh = in.getSDHeader();
-  SDMemTable* pTabOut = new SDMemTable(in, True);
-
-  // Find out how to convert values into Jy and K (e.g. units might be
-  // mJy or mK) Also automatically find out what we are converting to
-  // according to the flux unit
-  Unit fluxUnit(sh.fluxunit);
+  out->frequencies().rescale(width, "RESAMPLE");
+  TableIterator iter(tout, "IFNO");
+  TableRow row(tout);
+  while ( !iter.pastEnd() ) {
+    Table tab = iter.table();
+    ArrayColumn<Float> specCol(tab, "SPECTRA");
+    //ArrayColumn<Float> tsysCol(tout, "TSYS");
+    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+    Vector<Float> spec;
+    Vector<uChar> flag;
+    specCol.get(0,spec); // the number of channels should be constant per IF
+    uInt nChanIn = spec.nelements();
+    Vector<Float> xIn(nChanIn); indgen(xIn);
+    Int fac =  Int(nChanIn/width);
+    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
+    uInt k = 0;
+    Float x = 0.0;
+    while (x < Float(nChanIn) ) {
+      xOut(k) = x;
+      k++;
+      x += width;
+    }
+    uInt nChanOut = k;
+    xOut.resize(nChanOut, True);
+    // process all rows for this IFNO
+    Vector<Float> specOut;
+    Vector<Bool> maskOut;
+    Vector<uChar> flagOut;
+    for (uInt i=0; i < tab.nrow(); ++i) {
+      specCol.get(i, spec);
+      flagCol.get(i, flag);
+      Vector<Bool> mask(flag.nelements());
+      convertArray(mask, flag);
+
+      IPosition shapeIn(spec.shape());
+      //sh.nchan = nChanOut;
+      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
+                                                   xIn, spec, mask,
+                                                   interpMethod, True, True);
+      /// @todo do the same for channel based Tsys
+      flagOut.resize(maskOut.nelements());
+      convertArray(flagOut, maskOut);
+      specCol.put(i, specOut);
+      flagCol.put(i, flagOut);
+    }
+    ++iter;
+  }
+
+  return out;
+}
+
+STMath::imethod STMath::stringToIMethod(const std::string& in)
+{
+  static STMath::imap lookup;
+
+  // initialize the lookup table if necessary
+  if ( lookup.empty() ) {
+    lookup["NEAR"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
+    lookup["LIN"] = InterpolateArray1D<Double,Float>::linear;
+    lookup["CUB"]  = InterpolateArray1D<Double,Float>::cubic;
+    lookup["SPL"]  = InterpolateArray1D<Double,Float>::spline;
+  }
+
+  STMath::imap::const_iterator iter = lookup.find(in);
+
+  if ( lookup.end() == iter ) {
+    std::string message = in;
+    message += " is not a valid interpolation mode";
+    throw(AipsError(message));
+  }
+  return iter->second;
+}
+
+WeightType STMath::stringToWeight(const std::string& in)
+{
+  static std::map<std::string, WeightType> lookup;
+
+  // 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;
+  }
+
+  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
+
+  if ( lookup.end() == iter ) {
+    std::string message = in;
+    message += " is not a valid weighting mode";
+    throw(AipsError(message));
+  }
+  return iter->second;
+}
+
+CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
+                                               const Vector< Float > & coeffs,
+                                               const std::string & filename,
+                                               const std::string& method)
+{
+  // Get elevation data from Scantable and convert to degrees
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tab = in->table();
+  ROScalarColumn<Float> elev(tab, "ELEVATION");
+  Vector<Float> x = elev.getColumn();
+  x *= Float(180 / C::pi);                        // Degrees
+
+  const uInt nc = coeffs.nelements();
+  if ( filename.length() > 0 && nc > 0 ) {
+    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
+  }
+
+  // Correct
+  if ( nc > 0 || filename.length() == 0 ) {
+    // Find instrument
+    Bool throwit = True;
+    Instrument inst =
+      SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
+                                throwit);
+
+    // Set polynomial
+    Polynomial<Float>* ppoly = 0;
+    Vector<Float> coeff;
+    String msg;
+    if ( nc > 0 ) {
+      ppoly = new Polynomial<Float>(nc);
+      coeff = coeffs;
+      msg = String("user");
+    } else {
+      SDAttr sdAttr;
+      coeff = sdAttr.gainElevationPoly(inst);
+      ppoly = new Polynomial<Float>(3);
+      msg = String("built in");
+    }
+
+    if ( coeff.nelements() > 0 ) {
+      ppoly->setCoefficients(coeff);
+    } else {
+      delete ppoly;
+      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
+    }
+    ostringstream oss;
+    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
+    oss << "   " <<  coeff;
+    pushLog(String(oss));
+    const uInt nrow = tab.nrow();
+    Vector<Float> factor(nrow);
+    for ( uInt i=0; i < nrow; ++i ) {
+      factor[i] = 1.0 / (*ppoly)(x[i]);
+    }
+    delete ppoly;
+    scaleByVector(tab, factor, true);
+
+  } else {
+    // Read and correct
+    pushLog("Making correction from ascii Table");
+    scaleFromAsciiTable(tab, filename, method, x, true);
+  }
+  return out;
+}
+
+void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
+                                 const std::string& method,
+                                 const Vector<Float>& xout, bool dotsys)
+{
+
+// Read gain-elevation ascii file data into a Table.
+
+  String formatString;
+  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
+  scaleFromTable(in, tbl, method, xout, dotsys);
+}
+
+void STMath::scaleFromTable(Table& in,
+                            const Table& table,
+                            const std::string& method,
+                            const Vector<Float>& xout, bool dotsys)
+{
+
+  ROScalarColumn<Float> geElCol(table, "ELEVATION");
+  ROScalarColumn<Float> geFacCol(table, "FACTOR");
+  Vector<Float> xin = geElCol.getColumn();
+  Vector<Float> yin = geFacCol.getColumn();
+  Vector<Bool> maskin(xin.nelements(),True);
+
+  // Interpolate (and extrapolate) with desired method
+
+   //InterpolateArray1D<Double,Float>::InterpolationMethod method;
+   Int intmethod(stringToIMethod(method));
+
+   Vector<Float> yout;
+   Vector<Bool> maskout;
+   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
+                                                xin, yin, maskin, intmethod,
+                                                True, True);
+
+   scaleByVector(in, Float(1.0)/yout, dotsys);
+}
+
+void STMath::scaleByVector( Table& in,
+                            const Vector< Float >& factor,
+                            bool dotsys )
+{
+  uInt nrow = in.nrow();
+  if ( factor.nelements() != nrow ) {
+    throw(AipsError("factors.nelements() != table.nelements()"));
+  }
+  ArrayColumn<Float> specCol(in, "SPECTRA");
+  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
+  ArrayColumn<Float> tsysCol(in, "TSYS");
+  for (uInt i=0; i < nrow; ++i) {
+    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
+    ma *= factor[i];
+    specCol.put(i, ma.getArray());
+    flagCol.put(i, flagsFromMA(ma));
+    if ( dotsys ) {
+      Vector<Float> tsys;
+      tsysCol.get(i, tsys);
+      tsys *= factor[i];
+      specCol.put(i,tsys);
+    }
+  }
+}
+
+CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
+                                             float d, float etaap,
+                                             float jyperk )
+{
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tab = in->table();
+  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
   Unit K(String("K"));
   Unit JY(String("Jy"));
 
-  Bool toKelvin = True;
-  Double cFac = 1.0;
-
-  if (fluxUnit==JY) {
+  bool tokelvin = true;
+  Double cfac = 1.0;
+
+  if ( fluxUnit == JY ) {
     pushLog("Converting to K");
     Quantum<Double> t(1.0,fluxUnit);
     Quantum<Double> t2 = t.get(JY);
-    cFac = (t2 / t).getValue();               // value to Jy
-
-    toKelvin = True;
-    sh.fluxunit = "K";
-  } else if (fluxUnit==K) {
+    cfac = (t2 / t).getValue();               // value to Jy
+
+    tokelvin = true;
+    out->setFluxUnit("K");
+  } else if ( fluxUnit == K ) {
     pushLog("Converting to Jy");
     Quantum<Double> t(1.0,fluxUnit);
     Quantum<Double> t2 = t.get(K);
-    cFac = (t2 / t).getValue();              // value to K
-
-    toKelvin = False;
-    sh.fluxunit = "Jy";
+    cfac = (t2 / t).getValue();              // value to K
+
+    tokelvin = false;
+    out->setFluxUnit("Jy");
   } else {
     throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
   }
-
-  pTabOut->putSDHeader(sh);
-
   // Make sure input values are converted to either Jy or K first...
-  Float factor = cFac;
+  Float factor = cfac;
 
   // Select method
-  if (JyPerK>0.0) {
-    factor *= JyPerK;
-    if (toKelvin) factor = 1.0 / JyPerK;
+  if (jyperk > 0.0) {
+    factor *= jyperk;
+    if ( tokelvin ) factor = 1.0 / jyperk;
     ostringstream oss;
-    oss << "Jy/K = " << JyPerK;
+    oss << "Jy/K = " << jyperk;
     pushLog(String(oss));
-    Vector<Float> factors(in.nRow(), factor);
-    scaleByVector(pTabOut, in, doAll, factors, False);
-  } else if (etaAp>0.0) {
-    Bool throwIt = True;
-    Instrument inst = SDAttr::convertInstrument(sh.antennaname, throwIt);
+    Vector<Float> factors(tab.nrow(), factor);
+    scaleByVector(tab,factors, false);
+  } else if ( etaap > 0.0) {
+    Instrument inst =
+      SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
     SDAttr sda;
-    if (D < 0) D = sda.diameter(inst);
-    Float JyPerK = SDAttr::findJyPerK(etaAp,D);
+    if (d < 0) d = sda.diameter(inst);
+    Float jyPerk = SDAttr::findJyPerK(etaap, d);
     ostringstream oss;
-    oss << "Jy/K = " << JyPerK;
+    oss << "Jy/K = " << jyperk;
     pushLog(String(oss));
-    factor *= JyPerK;
-    if (toKelvin) {
+    factor *= jyperk;
+    if ( tokelvin ) {
       factor = 1.0 / factor;
     }
-
-    Vector<Float> factors(in.nRow(), factor);
-    scaleByVector(pTabOut, in, doAll, factors, False);
+    Vector<Float> factors(tab.nrow(), factor);
+    scaleByVector(tab, factors, False);
   } else {
 
@@ -1128 +701 @@
 
     pushLog("Looking up conversion factors");
-    convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
-  }
-  return pTabOut;
-}
-
-
-SDMemTable* SDMath::gainElevation(const SDMemTable& in,
-                                  const Vector<Float>& coeffs,
-                                  const String& fileName,
-                                  const String& methodStr, Bool doAll)
-{
-
-  // Get header and clone output table
-  SDHeader sh = in.getSDHeader();
-  SDMemTable* pTabOut = new SDMemTable(in, True);
-
-  // Get elevation data from SDMemTable and convert to degrees
-  const Table& tab = in.table();
+    convertBrightnessUnits(out, tokelvin, cfac);
+  }
+
+  return out;
+}
+
+void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
+                                     bool tokelvin, float cfac )
+{
+  Table& table = in->table();
+  Instrument inst =
+    SDAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
+  TableIterator iter(table, "FREQ_ID");
+  STFrequencies stfreqs = in->frequencies();
+  SDAttr sdAtt;
+  while (!iter.pastEnd()) {
+    Table tab = iter.table();
+    ArrayColumn<Float> specCol(tab, "SPECTRA");
+    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
+    MEpoch::ROScalarColumn timeCol(tab, "TIME");
+
+    uInt freqid; freqidCol.get(0, freqid);
+    Vector<Float> tmpspec; specCol.get(0, tmpspec);
+    // SDAttr.JyPerK has a Vector interface... change sometime.
+    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
+    for ( uInt i=0; i<tab.nrow(); ++i) {
+      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
+      Float factor = cfac * jyperk;
+      if ( tokelvin ) factor = Float(1.0) / factor;
+      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
+      ma *= factor;
+      specCol.put(i, ma.getArray());
+      flagCol.put(i, flagsFromMA(ma));
+    }
+  }
+}
+
+CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
+                                         float tau )
+{
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tab = in->table();
   ROScalarColumn<Float> elev(tab, "ELEVATION");
-  Vector<Float> x = elev.getColumn();
-  x *= Float(180 / C::pi);                        // Degrees
-
-  const uInt nC = coeffs.nelements();
-  if (fileName.length()>0 && nC>0) {
-    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
-  }
-
-  // Correct
-  if (nC>0 || fileName.length()==0) {
-    // Find instrument
-     Bool throwIt = True;
-     Instrument inst = SDAttr::convertInstrument (sh.antennaname, throwIt);
-
-     // Set polynomial
-     Polynomial<Float>* pPoly = 0;
-     Vector<Float> coeff;
-     String msg;
-     if (nC>0) {
-       pPoly = new Polynomial<Float>(nC);
-       coeff = coeffs;
-       msg = String("user");
-     } else {
-       SDAttr sdAttr;
-       coeff = sdAttr.gainElevationPoly(inst);
-       pPoly = new Polynomial<Float>(3);
-       msg = String("built in");
-     }
-
-     if (coeff.nelements()>0) {
-       pPoly->setCoefficients(coeff);
-     } else {
-       delete pPoly;
-       throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
-     }
-     ostringstream oss;
-     oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
-     oss << "   " <<  coeff;
-     pushLog(String(oss));
-     const uInt nRow = in.nRow();
-     Vector<Float> factor(nRow);
-     for (uInt i=0; i<nRow; i++) {
-       factor[i] = 1.0 / (*pPoly)(x[i]);
-     }
-     delete pPoly;
-     scaleByVector (pTabOut, in, doAll, factor, True);
-
-  } else {
-
-    // Indicate which columns to read from ascii file
-    String col0("ELEVATION");
-    String col1("FACTOR");
-
-    // Read and correct
-
-    pushLog("Making correction from ascii Table");
-    scaleFromAsciiTable (pTabOut, in, fileName, col0, col1,
-                         methodStr, doAll, x, True);
-  }
-
-  return pTabOut;
-}
-
-
-SDMemTable* SDMath::opacity(const SDMemTable& in, Float tau, Bool doAll)
-{
-
-  // Get header and clone output table
-
-  SDHeader sh = in.getSDHeader();
-  SDMemTable* pTabOut = new SDMemTable(in, True);
-
-// Get elevation data from SDMemTable and convert to degrees
-
-  const Table& tab = in.table();
-  ROScalarColumn<Float> elev(tab, "ELEVATION");
-  Vector<Float> zDist = elev.getColumn();
-  zDist = Float(C::pi_2) - zDist;
-
-// Generate correction factor
-
-  const uInt nRow = in.nRow();
-  Vector<Float> factor(nRow);
-  Vector<Float> factor2(nRow);
-  for (uInt i=0; i<nRow; i++) {
-     factor[i] = exp(tau)/cos(zDist[i]);
-  }
-
-// Correct
-
-  scaleByVector (pTabOut, in, doAll, factor, True);
-
-  return pTabOut;
-}
-
-
-void SDMath::rotateXYPhase(SDMemTable& in, Float value, Bool doAll)
-//
-// phase in degrees
-// assumes linear correlations
-//
-{
-  if (in.nPol() != 4) {
-    throw(AipsError("You must have 4 polarizations to run this function"));
-  }
-
-   SDHeader sh = in.getSDHeader();
-   Instrument inst = SDAttr::convertInstrument(sh.antennaname, False);
-   SDAttr sdAtt;
-   if (sdAtt.feedPolType(inst) != LINEAR) {
-      throw(AipsError("Only linear polarizations are supported"));
-   }
-//
-   const Table& tabIn = in.table();
-   ArrayColumn<Float> specCol(tabIn,"SPECTRA");
-   IPosition start(asap::nAxes,0);
-   IPosition end(asap::nAxes);
-
-// Set cursor slice. Assumes shape the same for all rows
-
-   setCursorSlice (start, end, doAll, in);
-   IPosition start3(start);
-   start3(asap::PolAxis) = 2;                 // Real(XY)
-   IPosition end3(end);
-   end3(asap::PolAxis) = 2;
-//
-   IPosition start4(start);
-   start4(asap::PolAxis) = 3;                 // Imag (XY)
-   IPosition end4(end);
-   end4(asap::PolAxis) = 3;
-//
-   uInt nRow = in.nRow();
-   Array<Float> data;
-   for (uInt i=0; i<nRow;++i) {
-      specCol.get(i,data);
-      IPosition shape = data.shape();
-
-      // Get polarization slice references
-      Array<Float> C3 = data(start3,end3);
-      Array<Float> C4 = data(start4,end4);
-
-      // Rotate
-      SDPolUtil::rotatePhase(C3, C4, value);
-
-      // Put
-      specCol.put(i,data);
-   }
-}
-
-
-void SDMath::rotateLinPolPhase(SDMemTable& in, Float value, Bool doAll)
-//
-// phase in degrees
-// assumes linear correlations
-//
-{
-   if (in.nPol() != 4) {
-      throw(AipsError("You must have 4 polarizations to run this function"));
-   }
-//
-   SDHeader sh = in.getSDHeader();
-   Instrument inst = SDAttr::convertInstrument(sh.antennaname, False);
-   SDAttr sdAtt;
-   if (sdAtt.feedPolType(inst) != LINEAR) {
-      throw(AipsError("Only linear polarizations are supported"));
-   }
-//
-   const Table& tabIn = in.table();
-   ArrayColumn<Float> specCol(tabIn,"SPECTRA");
-   ROArrayColumn<Float> stokesCol(tabIn,"STOKES");
-   IPosition start(asap::nAxes,0);
-   IPosition end(asap::nAxes);
-
-// Set cursor slice. Assumes shape the same for all rows
-
-   setCursorSlice (start, end, doAll, in);
-//
-   IPosition start1(start);
-   start1(asap::PolAxis) = 0;                // C1 (XX)
-   IPosition end1(end);
-   end1(asap::PolAxis) = 0;
-//
-   IPosition start2(start);
-   start2(asap::PolAxis) = 1;                 // C2 (YY)
-   IPosition end2(end);
-   end2(asap::PolAxis) = 1;
-//
-   IPosition start3(start);
-   start3(asap::PolAxis) = 2;                 // C3 ( Real(XY) )
-   IPosition end3(end);
-   end3(asap::PolAxis) = 2;
-//
-   IPosition startI(start);
-   startI(asap::PolAxis) = 0;                 // I
-   IPosition endI(end);
-   endI(asap::PolAxis) = 0;
-//
-   IPosition startQ(start);
-   startQ(asap::PolAxis) = 1;                 // Q
-   IPosition endQ(end);
-   endQ(asap::PolAxis) = 1;
-//
-   IPosition startU(start);
-   startU(asap::PolAxis) = 2;                 // U
-   IPosition endU(end);
-   endU(asap::PolAxis) = 2;
-
-//
-   uInt nRow = in.nRow();
-   Array<Float> data, stokes;
-   for (uInt i=0; i<nRow;++i) {
-      specCol.get(i,data);
-      stokesCol.get(i,stokes);
-      IPosition shape = data.shape();
-
-// Get linear polarization slice references
-
-      Array<Float> C1 = data(start1,end1);
-      Array<Float> C2 = data(start2,end2);
-      Array<Float> C3 = data(start3,end3);
-
-// Get STokes slice references
-
-      Array<Float> I = stokes(startI,endI);
-      Array<Float> Q = stokes(startQ,endQ);
-      Array<Float> U = stokes(startU,endU);
-
-// Rotate
-
-      SDPolUtil::rotateLinPolPhase(C1, C2, C3, I, Q, U, value);
-
-// Put
-
-      specCol.put(i,data);
-   }
-}
-
-// 'private' functions
-
-void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
-                                     Bool toKelvin, Float cFac, Bool doAll)
-{
-
-// Get header
-
-   SDHeader sh = in.getSDHeader();
-   const uInt nChan = sh.nchan;
-
-// Get instrument
-
-   Bool throwIt = True;
-   Instrument inst = SDAttr::convertInstrument(sh.antennaname, throwIt);
-
-// Get Diameter (m)
-
-   SDAttr sdAtt;
-// Get epoch of first row
-
-   MEpoch dateObs = in.getEpoch(0);
-
-// Generate a Vector of correction factors. One per FreqID
-
-   SDFrequencyTable sdft = in.getSDFreqTable();
-   Vector<uInt> freqIDs;
-//
-   Vector<Float> freqs(sdft.length());
-   for (uInt i=0; i<sdft.length(); i++) {
-      freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
-   }
-//
-   Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
-   ostringstream oss;
-   oss << "Jy/K = " << JyPerK;
-   pushLog(String(oss));
-   Vector<Float> factors = cFac * JyPerK;
-   if (toKelvin) factors = Float(1.0) / factors;
-
-// Get data slice bounds
-
-   IPosition start, end;
-   setCursorSlice (start, end, doAll, in);
-   const uInt ifAxis = in.getIF();
-
-// Iteration axes
-
-   IPosition axes(asap::nAxes-1,0);
-   for (uInt i=0,j=0; i<asap::nAxes; i++) {
-      if (i!=asap::IFAxis) {
-         axes(j++) = i;
+  ArrayColumn<Float> specCol(tab, "SPECTRA");
+  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+  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));
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
+                                        const std::string& kernel, float width )
+{
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& table = in->table();
+  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
+  // same IFNO should have same no of channels
+  // this saves overhead
+  TableIterator iter(table, "IFNO");
+  while (!iter.pastEnd()) {
+    Table tab = iter.table();
+    ArrayColumn<Float> specCol(tab, "SPECTRA");
+    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+    Vector<Float> tmpspec; specCol.get(0, tmpspec);
+    uInt nchan = tmpspec.nelements();
+    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
+    Convolver<Float> conv(kvec, IPosition(1,nchan));
+    Vector<Float> spec;
+    Vector<uChar> flag;
+    for ( uInt i=0; i<tab.nrow(); ++i) {
+      specCol.get(i, spec);
+      flagCol.get(i, flag);
+      Vector<Bool> mask(flag.nelements());
+      convertArray(mask, flag);
+      Vector<Float> specout;
+      if ( type == VectorKernel::HANNING ) {
+        Vector<Bool> maskout;
+        mathutil::hanning(specout, maskout, spec , mask);
+        convertArray(flag, maskout);
+        flagCol.put(i, flag);
+      } else {
+        mathutil::replaceMaskByZero(specout, mask);
+        conv.linearConv(specout, spec);
       }
-   }
-
-// Loop over rows and apply correction factor
-
-   Float factor = 1.0;
-   const uInt axis = asap::ChanAxis;
-   for (uInt i=0; i < in.nRow(); ++i) {
-
-// Get data
-
-      MaskedArray<Float> dataIn = in.rowAsMaskedArray(i);
-      Array<Float>& values = dataIn.getRWArray();           // Ref to dataIn
-      Array<Float> values2 = values(start,end);             // Ref to values to dataIn
-
-// Get SDCOntainer
-
-      SDContainer sc = in.getSDContainer(i);
-
-// Get FreqIDs
-
-      freqIDs = sc.getFreqMap();
-
-// Now the conversion factor depends only upon frequency
-// So we need to iterate through by IF only giving
-// us BEAM/POL/CHAN cubes
-
-      ArrayIterator<Float> itIn(values2, axes);
-      uInt ax = 0;
-      while (!itIn.pastEnd()) {
-        itIn.array() *= factors(freqIDs(ax));         // Writes back to dataIn
-        itIn.next();
-      }
-
-// Write out
-
-      putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
-//
-      pTabOut->putSDContainer(sc);
-   }
-}
-
-
-
-SDMemTable* SDMath::frequencyAlign(const SDMemTable& in,
-                                   MFrequency::Types freqSystem,
-                                   const String& refTime,
-                                   const String& methodStr,
-                                   Bool perFreqID)
-{
-// Get Header
-
-   SDHeader sh = in.getSDHeader();
-   const uInt nChan = sh.nchan;
-   const uInt nRows = in.nRow();
-   const uInt nIF = sh.nif;
-
-// 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();
-
-// Generate DataDesc table
-
-   Matrix<uInt> ddIdx;
-   SDDataDesc dDesc;
-   generateDataDescTable(ddIdx, dDesc, nIF, in, tabIn, srcCol,
-                          fqIDCol, perFreqID);
-
-// Get reference Epoch to time of first row or given String
-
-   Unit DAY(String("d"));
-   MEpoch::Ref epochRef(in.getTimeReference());
-   MEpoch refEpoch;
-   if (refTime.length()>0) {
-     refEpoch = epochFromString(refTime, in.getTimeReference());
-   } else {
-     refEpoch = in.getEpoch(0);
-   }
-   ostringstream oss;
-   oss << "Aligned at reference Epoch " << formatEpoch(refEpoch) << " in frame " << MFrequency::showType(freqSystem);
-   pushLog(String(oss));
-
-   // Get Reference Position
-
-   MPosition refPos = in.getAntennaPosition();
-
-   // Create FrequencyAligner Block. One FA for each possible
-   // source/freqID (perFreqID=True) or source/IF (perFreqID=False)
-   // combination
-
-   PtrBlock<FrequencyAligner<Float>* > a(dDesc.length());
-   generateFrequencyAligners(a, dDesc, in, nChan, freqSystem, refPos,
-                              refEpoch, perFreqID);
-
-   // Generate and fill output Frequency Table.  WHen perFreqID=True,
-   // there is one output FreqID for each entry in the SDDataDesc
-   // table.  However, in perFreqID=False mode, there may be some
-   // degeneracy, so we need a little translation map
-
-   SDFrequencyTable freqTabOut = in.getSDFreqTable();
-   freqTabOut.setLength(0);
-   Vector<String> units(1);
-   units = String("Hz");
-   Bool linear=True;
-   //
-   Vector<uInt> ddFQTrans(dDesc.length(),0);
-   for (uInt i=0; i<dDesc.length(); i++) {
-
-     // Get Aligned SC in Hz
-
-     SpectralCoordinate sC = a[i]->alignedSpectralCoordinate(linear);
-     sC.setWorldAxisUnits(units);
-
-     // Add FreqID
-
-     uInt idx = freqTabOut.addFrequency(sC.referencePixel()[0],
-                                        sC.referenceValue()[0],
-                                        sC.increment()[0]);
-     // output FreqID = ddFQTrans(ddIdx)
-     ddFQTrans(i) = idx;
-   }
-
-   // Interpolation method
-
-   InterpolateArray1D<Double,Float>::InterpolationMethod interp;
-   convertInterpString(interp, methodStr);
-
-   // New output Table
-
-   //pushLog("Create output table");
-   SDMemTable* pTabOut = new SDMemTable(in,True);
-   pTabOut->putSDFreqTable(freqTabOut);
-
-   // Loop over rows in Table
-
-   Bool extrapolate=False;
-   const IPosition polChanAxes(2, asap::PolAxis, asap::ChanAxis);
-   Bool useCachedAbcissa = False;
-   Bool first = True;
-   Bool ok;
-   Vector<Float> yOut;
-   Vector<Bool> maskOut;
-   Vector<uInt> freqID(nIF);
-   uInt ifIdx, faIdx;
-   Vector<Double> xIn;
-   //
-   for (uInt iRow=0; iRow<nRows; ++iRow) {
-//      if (iRow%10==0) {
-//        ostringstream oss;
-//        oss << "Processing row " << iRow;
-//        pushLog(String(oss));
-//      }
-
-     // Get EPoch
-
-     Quantum<Double> tQ2(times[iRow],DAY);
-     MVEpoch mv2(tQ2);
-     MEpoch epoch(mv2, epochRef);
-
-     // Get copy of data
-
-     const MaskedArray<Float>& mArrIn(in.rowAsMaskedArray(iRow));
-     Array<Float> values = mArrIn.getArray();
-     Array<Bool> mask = mArrIn.getMask();
-
-     // For each row, the Frequency 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.  Probably won't work for
-     // multiple beams at this point.
-
-     ArrayIterator<Float> itValuesPlane(values, polChanAxes);
-     ArrayIterator<Bool> itMaskPlane(mask, polChanAxes);
-     while (!itValuesPlane.pastEnd()) {
-
-       // Find the IF index and then the FA PtrBlock index
-
-       const IPosition& pos = itValuesPlane.pos();
-       ifIdx = pos(asap::IFAxis);
-       faIdx = ddIdx(iRow,ifIdx);
-
-       // Generate abcissa for perIF.  Could cache this in a Matrix on
-       // a per scan basis.  Pretty expensive doing it for every row.
-
-       if (!perFreqID) {
-         xIn.resize(nChan);
-         uInt fqID = dDesc.secID(ddIdx(iRow,ifIdx));
-         SpectralCoordinate sC = in.getSpectralCoordinate(fqID);
-           Double w;
-           for (uInt i=0; i<nChan; i++) {
-             sC.toWorld(w,Double(i));
-              xIn[i] = w;
-           }
-       }
-
-       VectorIterator<Float> itValuesVec(itValuesPlane.array(), 1);
-       VectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
-
-       // Iterate through the plane by vector and align
-
-        first = True;
-        useCachedAbcissa=False;
-        while (!itValuesVec.pastEnd()) {
-          if (perFreqID) {
-            ok = a[faIdx]->align (yOut, maskOut, itValuesVec.vector(),
-                                  itMaskVec.vector(), epoch, useCachedAbcissa,
-                                  interp, extrapolate);
-          } else {
-            ok = a[faIdx]->align (yOut, maskOut, xIn, itValuesVec.vector(),
-                                  itMaskVec.vector(), epoch, useCachedAbcissa,
-                                  interp, extrapolate);
-          }
-          //
-          itValuesVec.vector() = yOut;
-          itMaskVec.vector() = maskOut;
-          //
-          itValuesVec.next();
-          itMaskVec.next();
-          //
-          if (first) {
-            useCachedAbcissa = True;
-            first = False;
-          }
-        }
-        //
-        itValuesPlane.next();
-        itMaskPlane.next();
-     }
-
-     // Create SDContainer and put back
-
-     SDContainer sc = in.getSDContainer(iRow);
-     putDataInSDC(sc, values, mask);
-
-     // Set output FreqIDs
-
-     for (uInt i=0; i<nIF; i++) {
-       uInt idx = ddIdx(iRow,i);               // Index into SDDataDesc table
-       freqID(i) = ddFQTrans(idx);             // FreqID in output FQ table
-     }
-     sc.putFreqMap(freqID);
-     //
-     pTabOut->putSDContainer(sc);
-   }
-
-   // Now we must set the base and extra frames to the input frame
-   std::vector<string> info = pTabOut->getCoordInfo();
-   info[1] = MFrequency::showType(freqSystem);   // Conversion frame
-   info[3] = info[1];                            // Base frame
-   pTabOut->setCoordInfo(info);
-
-   // Clean up PointerBlock
-   for (uInt i=0; i<a.nelements(); i++) delete a[i];
-
-   return pTabOut;
-}
-
-
-SDMemTable* SDMath::frequencySwitch(const SDMemTable& in)
-{
-  if (in.nIF() != 2) {
-    throw(AipsError("nIF != 2 "));
-  }
-  Bool clear = True;
-  SDMemTable* pTabOut = new SDMemTable(in, clear);
-  const Table& tabIn = in.table();
-
-  // Shape of input and output data
-  const IPosition& shapeIn = in.rowAsMaskedArray(0).shape();
-
-  const uInt nRows = in.nRow();
-  AlwaysAssert(asap::nAxes==4,AipsError);
-
-  ROArrayColumn<Float> tSysCol;
-  Array<Float> tsys;
-  tSysCol.attach(tabIn,"TSYS");
-
-  for (uInt iRow=0; iRow<nRows; iRow++) {
-    // Get data for this row
-    MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
-    tSysCol.get(iRow, tsys);
-
-    // whole Array for IF 0
-    IPosition start(asap::nAxes,0);
-    IPosition end = shapeIn-1;
-    end(asap::IFAxis) = 0;
-
-    MaskedArray<Float> on = marr(start,end);
-    Array<Float> ton = tsys(start,end);
-    // Make a copy as "src" is a refrence which is manipulated.
-    // oncopy is needed for the inverse quotient
-    MaskedArray<Float> oncopy = on.copy();
-
-    // whole Array for IF 1
-    start(asap::IFAxis) = 1;
-    end(asap::IFAxis) = 1;
-
-    MaskedArray<Float> off = marr(start,end);
-    Array<Float> toff = tsys(start,end);
-
-    on /= off; on -= 1.0f;
-    on *= ton;
-    off /= oncopy; off -= 1.0f;
-    off *= toff;
-
-    SDContainer sc = in.getSDContainer(iRow);
-    putDataInSDC(sc, marr.getArray(), marr.getMask());
-    pTabOut->putSDContainer(sc);
-  }
-  return pTabOut;
-}
-
-void SDMath::fillSDC(SDContainer& sc,
-                     const Array<Bool>& mask,
-                     const Array<Float>& data,
-                     const Array<Float>& tSys,
-                     Int scanID, Double timeStamp,
-                     Double interval, const String& sourceName,
-                     const Vector<uInt>& freqID)
-{
-// Data and mask
-
-  putDataInSDC(sc, data, mask);
-
-// TSys
-
-  sc.putTsys(tSys);
-
-// Time things
-
-  sc.timestamp = timeStamp;
-  sc.interval = interval;
-  sc.scanid = scanID;
-//
-  sc.sourcename = sourceName;
-  sc.putFreqMap(freqID);
-}
-
-void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
-                        MaskedArray<Float>& sum, Array<Float>& sumSq,
-                        Array<Float>& nPts, Array<Float>& tSysSum,
-                        Array<Float>& tSysSqSum,
-                        const Array<Float>& tSys, const Array<Float>& nInc,
-                        const Vector<Bool>& mask, Double time, Double interval,
-                        const std::vector<CountedPtr<SDMemTable> >& in,
-                        uInt iTab, uInt iRow, uInt axis,
-                        uInt nAxesSub, Bool useMask,
-                        WeightType wtType)
-{
-
-// Get data
-
-   MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
-   Array<Float>& valuesIn = dataIn.getRWArray();           // writable reference
-   const Array<Bool>& maskIn = dataIn.getMask();          // RO reference
-//
-   if (wtType==NONE) {
-      const MaskedArray<Float> n(nInc,dataIn.getMask());
-      nPts += n;                               // Only accumulates where mask==T
-   } else if (wtType==TINT) {
-
-// We are weighting the data by integration time.
-
-     valuesIn *= Float(interval);
-
-   } else if (wtType==VAR) {
-
-// We are going to average the data, weighted by the noise for each pol, beam and IF.
-// So therefore we need to iterate through by spectrum (axis 3)
-
-      VectorIterator<Float> itData(valuesIn, axis);
-      ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
-      Float fac = 1.0;
-      IPosition pos(nAxesSub,0);
-//
-      while (!itData.pastEnd()) {
-
-// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
-
-         if (useMask) {
-            MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
-            fac = 1.0/variance(tmp);
-         } else {
-            MaskedArray<Float> tmp(itData.vector(),itMask.vector());
-            fac = 1.0/variance(tmp);
-         }
-
-// Scale data
-
-         itData.vector() *= fac;     // Writes back into 'dataIn'
-//
-// Accumulate variance per if/pol/beam averaged over spectrum
-// This method to get pos2 from itData.pos() is only valid
-// because the spectral axis is the last one (so we can just
-// copy the first nAXesSub positions out)
-
-         pos = itData.pos().getFirst(nAxesSub);
-         sumSq(pos) += fac;
-//
-         itData.next();
-         itMask.next();
-      }
-   } else if (wtType==TSYS || wtType==TINTSYS) {
-
-// We are going to average the data, weighted by 1/Tsys**2 for each pol, beam and IF.
-// So therefore we need to iterate through by spectrum (axis 3).  Although
-// Tsys is stored as a vector of length nChan, the values are replicated.
-// We will take a short cut and just use the value from the first channel
-// for now.
-//
-      VectorIterator<Float> itData(valuesIn, axis);
-      ReadOnlyVectorIterator<Float> itTSys(tSys, axis);
-      IPosition pos(nAxesSub,0);
-//
-      Float fac = 1.0;
-      if (wtType==TINTSYS) fac *= interval;
-      while (!itData.pastEnd()) {
-         Float t = itTSys.vector()[0];
-         fac *= 1.0/t/t;
-
-// Scale data
-
-         itData.vector() *= fac;     // Writes back into 'dataIn'
-//
-// Accumulate Tsys  per if/pol/beam averaged over spectrum
-// This method to get pos2 from itData.pos() is only valid
-// because the spectral axis is the last one (so we can just
-// copy the first nAXesSub positions out)
-
-         pos = itData.pos().getFirst(nAxesSub);
-         tSysSqSum(pos) += fac;
-//
-         itData.next();
-         itTSys.next();
-      }
-   }
-
-// Accumulate sum of (possibly scaled) data
-
-   sum += dataIn;
-
-// Accumulate Tsys, time, and interval
-
-   tSysSum += tSys;
-   timeSum += time;
-   intSum += interval;
-   nAccum += 1;
-}
-
-
-void SDMath::normalize(MaskedArray<Float>& sum,
-                       const Array<Float>& sumSq,
-                       const Array<Float>& tSysSqSum,
-                       const Array<Float>& nPts,
-                       Double intSum,
-                       WeightType wtType, Int axis,
-                       Int nAxesSub)
-{
-   IPosition pos2(nAxesSub,0);
-//
-   if (wtType==NONE) {
-
-// We just average by the number of points accumulated.
-// We need to make a MA out of nPts so that no divide by
-// zeros occur
-
-      MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
-      sum /= t;
-   } else if (wtType==TINT) {
-
-// Average by sum of Tint
-
-      sum /= Float(intSum);
-   } else if (wtType==VAR) {
-
-// Normalize each spectrum by sum(1/var) where the variance
-// is worked out for each spectrum
-
-      Array<Float>& data = sum.getRWArray();
-      VectorIterator<Float> itData(data, axis);
-      while (!itData.pastEnd()) {
-         pos2 = itData.pos().getFirst(nAxesSub);
-         itData.vector() /= sumSq(pos2);
-         itData.next();
-      }
-   } else if (wtType==TSYS || wtType==TINTSYS) {
-
-// Normalize each spectrum by sum(1/Tsys**2) (TSYS) or
-// sum(Tint/Tsys**2) (TINTSYS) where the pseudo
-// replication over channel for Tsys has been dropped.
-
-      Array<Float>& data = sum.getRWArray();
-      VectorIterator<Float> itData(data, axis);
-      while (!itData.pastEnd()) {
-         pos2 = itData.pos().getFirst(nAxesSub);
-         itData.vector() /= tSysSqSum(pos2);
-         itData.next();
-      }
-   }
-}
-
-
-
-
-void SDMath::setCursorSlice (IPosition& start, IPosition& end, Bool doAll, const SDMemTable& in) const
-{
-  const uInt nDim = asap::nAxes;
-  DebugAssert(nDim==4,AipsError);
-//
-  start.resize(nDim);
-  end.resize(nDim);
-  if (doAll) {
-     start = 0;
-     end(0) = in.nBeam()-1;
-     end(1) = in.nIF()-1;
-     end(2) = in.nPol()-1;
-     end(3) = in.nChan()-1;
-  } else {
-     start(0) = in.getBeam();
-     end(0) = start(0);
-//
-     start(1) = in.getIF();
-     end(1) = start(1);
-//
-     start(2) = in.getPol();
-     end(2) = start(2);
-//
-     start(3) = 0;
-     end(3) = in.nChan()-1;
-   }
-}
-
-
-void SDMath::convertWeightString(WeightType& wtType, const String& weightStr,
-                                 Bool listType)
-{
-  String tStr(weightStr);
-  tStr.upcase();
-  String msg;
-  if (tStr.contains(String("NONE"))) {
-     wtType = NONE;
-     msg = String("Weighting type selected : None");
-  } else if (tStr.contains(String("VAR"))) {
-     wtType = VAR;
-     msg = String("Weighting type selected : Variance");
-  } else if (tStr.contains(String("TINTSYS"))) {
-       wtType = TINTSYS;
-       msg = String("Weighting type selected : Tint&Tsys");
-  } else if (tStr.contains(String("TINT"))) {
-     wtType = TINT;
-     msg = String("Weighting type selected : Tint");
-  } else if (tStr.contains(String("TSYS"))) {
-     wtType = TSYS;
-     msg = String("Weighting type selected : Tsys");
-  } else {
-     msg = String("Weighting type selected : None");
-     throw(AipsError("Unrecognized weighting type"));
-  }
-//
-  if (listType) pushLog(msg);
-}
-
-
-void SDMath::convertInterpString(casa::InterpolateArray1D<Double,Float>::InterpolationMethod& type,
-                                 const casa::String& interp)
-{
-  String tStr(interp);
-  tStr.upcase();
-  if (tStr.contains(String("NEAR"))) {
-     type = InterpolateArray1D<Double,Float>::nearestNeighbour;
-  } else if (tStr.contains(String("LIN"))) {
-     type = InterpolateArray1D<Double,Float>::linear;
-  } else if (tStr.contains(String("CUB"))) {
-     type = InterpolateArray1D<Double,Float>::cubic;
-  } else if (tStr.contains(String("SPL"))) {
-     type = InterpolateArray1D<Double,Float>::spline;
-  } else {
-    throw(AipsError("Unrecognized interpolation type"));
-  }
-}
-
-void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
-                          const Array<Bool>& mask)
-{
-    sc.putSpectrum(data);
-//
-    Array<uChar> outflags(data.shape());
-    convertArray(outflags,!mask);
-    sc.putFlags(outflags);
-}
-
-Table SDMath::readAsciiFile(const String& fileName) const
-{
-   String formatString;
-   Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
-   return tbl;
-}
-
-
-
-void SDMath::scaleFromAsciiTable(SDMemTable* pTabOut,
-                                 const SDMemTable& in, const String& fileName,
-                                 const String& col0, const String& col1,
-                                 const String& methodStr, Bool doAll,
-                                 const Vector<Float>& xOut, Bool doTSys)
-{
-
-// Read gain-elevation ascii file data into a Table.
-
-  Table geTable = readAsciiFile (fileName);
-//
-  scaleFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut, doTSys);
-}
-
-void SDMath::scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in,
-                            const Table& tTable, const String& col0,
-                            const String& col1,
-                            const String& methodStr, Bool doAll,
-                            const Vector<Float>& xOut, Bool doTsys)
-{
-
-// Get data from Table
-
-  ROScalarColumn<Float> geElCol(tTable, col0);
-  ROScalarColumn<Float> geFacCol(tTable, col1);
-  Vector<Float> xIn = geElCol.getColumn();
-  Vector<Float> yIn = geFacCol.getColumn();
-  Vector<Bool> maskIn(xIn.nelements(),True);
-
-// Interpolate (and extrapolate) with desired method
-
-   InterpolateArray1D<Double,Float>::InterpolationMethod method;
-   convertInterpString(method, methodStr);
-   Int intMethod(method);
-//
-   Vector<Float> yOut;
-   Vector<Bool> maskOut;
-   InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
-                                                xIn, yIn, maskIn, intMethod,
-                                                True, True);
-// Apply
-
-   scaleByVector(pTabOut, in, doAll, Float(1.0)/yOut, doTsys);
-}
-
-
-void SDMath::scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
-                           Bool doAll, const Vector<Float>& factor,
-                           Bool doTSys)
-{
-
-// Set up data slice
-
-  IPosition start, end;
-  setCursorSlice (start, end, doAll, in);
-
-// Get Tsys column
-
-  const Table& tIn = in.table();
-  ArrayColumn<Float> tSysCol(tIn, "TSYS");
-  Array<Float> tSys;
-
-// Loop over rows and apply correction factor
-
-  const uInt axis = asap::ChanAxis;
-  for (uInt i=0; i < in.nRow(); ++i) {
-
-// Get data
-
-     MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
-     MaskedArray<Float> dataIn2 = dataIn(start,end);  // reference to dataIn
-//
-     if (doTSys) {
-        tSysCol.get(i, tSys);
-        Array<Float> tSys2 = tSys(start,end) * factor[i];
-        tSysCol.put(i, tSys);
-     }
-
-// Apply factor
-
-     dataIn2 *= factor[i];
-
-// Write out
-
-     SDContainer sc = in.getSDContainer(i);
-     putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
-//
-     pTabOut->putSDContainer(sc);
-  }
-}
-
-
-
-
-void SDMath::generateDataDescTable (Matrix<uInt>& ddIdx,
-                                    SDDataDesc& dDesc,
-                                    uInt nIF,
-                                    const SDMemTable& in,
-                                    const Table& tabIn,
-                                    const ROScalarColumn<String>& srcCol,
-                                    const ROArrayColumn<uInt>& fqIDCol,
-                                    Bool perFreqID)
-{
-   const uInt nRows = tabIn.nrow();
-   ddIdx.resize(nRows,nIF);
-//
-   String srcName;
-   Vector<uInt> freqIDs;
-   for (uInt iRow=0; iRow<nRows; iRow++) {
-      srcCol.get(iRow, srcName);
-      fqIDCol.get(iRow, freqIDs);
-      const MDirection& dir = in.getDirection(iRow);
-//
-      if (perFreqID) {
-
-// One entry per source/freqID pair
-
-         for (uInt iIF=0; iIF<nIF; iIF++) {
-            ddIdx(iRow,iIF) = dDesc.addEntry(srcName, freqIDs[iIF], dir, 0);
-         }
-      } else {
-
-// One entry per source/IF pair.  Hang onto the FreqID as well
-
-         for (uInt iIF=0; iIF<nIF; iIF++) {
-            ddIdx(iRow,iIF) = dDesc.addEntry(srcName, iIF, dir, freqIDs[iIF]);
-         }
-      }
-   }
-}
-
-
-
-
-
-MEpoch SDMath::epochFromString(const String& str, MEpoch::Types timeRef)
-{
-   Quantum<Double> qt;
-   if (MVTime::read(qt,str)) {
-      MVEpoch mv(qt);
-      MEpoch me(mv, timeRef);
-      return me;
-   } else {
-      throw(AipsError("Invalid format for Epoch string"));
-   }
-}
-
-
-String SDMath::formatEpoch(const MEpoch& epoch)  const
-{
-   MVTime mvt(epoch.getValue());
-   return mvt.string(MVTime::YMD) + String(" (") + epoch.getRefString() + String(")");
-}
-
-
-
-void SDMath::generateFrequencyAligners(PtrBlock<FrequencyAligner<Float>* >& a,
-                                       const SDDataDesc& dDesc,
-                                       const SDMemTable& in, uInt nChan,
-                                       MFrequency::Types system,
-                                       const MPosition& refPos,
-                                       const MEpoch& refEpoch,
-                                       Bool perFreqID)
-{
-   for (uInt i=0; i<dDesc.length(); i++) {
-      uInt ID = dDesc.ID(i);
-      uInt secID = dDesc.secID(i);
-      const MDirection& refDir = dDesc.secDir(i);
-//
-      if (perFreqID) {
-
-// One aligner per source/FreqID pair.
-
-         SpectralCoordinate sC = in.getSpectralCoordinate(ID);
-         a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
-      } else {
-
-// One aligner per source/IF pair.  But we still need the FreqID to
-// get the right SC.  Hence the messing about with the secondary ID
-
-         SpectralCoordinate sC = in.getSpectralCoordinate(secID);
-         a[i] = new FrequencyAligner<Float>(sC, nChan, refEpoch, refDir, refPos, system);
-      }
-   }
-}
-
-Vector<uInt> SDMath::getRowRange(const SDMemTable& in) const
-{
-   Vector<uInt> range(2);
-   range[0] = 0;
-   range[1] = in.nRow()-1;
-   return range;
-}
-
-
-Bool SDMath::rowInRange(uInt i, const Vector<uInt>& range) const
-{
-   return (i>=range[0] && i<=range[1]);
-}
-
+      specCol.put(i, specout);
+    }
+  }
+  return out;
+}

trunk/src/STMath.h

@@ -1 @@
-//#---------------------------------------------------------------------------
-//# SDMath.h: A collection of single dish mathematical operations
-//#---------------------------------------------------------------------------
-//# Copyright (C) 2004
-//# ATNF
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the Free
-//# Software Foundation; either version 2 of the License, or (at your option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
-//# WITHOUT ANY WARRANTY; without even the implied warranty of
-//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-//# Public License for more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning this software should be addressed as follows:
-//#        Internet email: Malte.Marquarding@csiro.au
-//#        Postal address: Malte Marquarding,
-//#                        Australia Telescope National Facility,
-//#                        P.O. Box 76,
-//#                        Epping, NSW, 2121,
-//#                        AUSTRALIA
-//#
-//# $Id:
-//#---------------------------------------------------------------------------
-#ifndef SDMATH_H
-#define SDMATH_H
+//
+// C++ Interface: STMath
+//
+// Description:
+//
+//
+// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+#ifndef ASAPSTMATH_H
+#define ASAPSTMATH_H
 
 #include <string>
-#include <vector>
+#include <map>
+
 #include <casa/aips.h>
 #include <casa/Utilities/CountedPtr.h>
-#include <coordinates/Coordinates/FrequencyAligner.h>
+#include <casa/BasicSL/String.h>
+#include <casa/Arrays/Vector.h>
+#include <scimath/Mathematics/InterpolateArray1D.h>
 
+#include "Scantable.h"
 #include "SDDefs.h"
 #include "SDLog.h"
 
-class casa::Table;
-class casa::MEpoch;
-class casa::MPosition;
-template<class T> class casa::PtrBlock;
-template<class T> class casa::Matrix;
-template<class T> class casa::ROScalarColumn;
-template<class T> class casa::ROArrayColumn;
-
 namespace asap {
 
-class SDMemTable;
-class SDDataDesc;
-  
-class SDMath : private SDLog {
+/**
+Mathmatical operations on Scantable objects
 
+@author Malte Marquarding
+*/
+class STMath : private SDLog {
 public:
-  
-  // Default constructor 
-  SDMath();
-  
-  // Copy Constructor (copy semantics)
-  SDMath(const SDMath& other);
-  
-  // Assignment  (copy semantics)
-  SDMath &operator=(const SDMath& other);
-  
-  // Destructor
-  ~SDMath();
-  
-  // Binary Table operators. op=ADD, SUB, MUL, DIV, QUOTIENT
-  casa::CountedPtr<SDMemTable> binaryOperate(const casa::CountedPtr<SDMemTable>& left,
-                                             const casa::CountedPtr<SDMemTable>& right,
-                                             const casa::String& op, casa::Bool preserve,
-                                             casa::Bool tSys);
-    
-    // Average in time
-    casa::CountedPtr<SDMemTable> average(const std::vector<casa::CountedPtr<SDMemTable> >& in,
-                                         const casa::Vector<casa::Bool>& mask,
-                                         casa::Bool scanAverage, 
-                                         const casa::String& weightStr,
-                                         casa::Bool align=casa::False);
-    
-  // Statistics. If row<0, all rows are done otherwise, just the 
-  // specified row.
-  std::vector<float> statistic(const casa::CountedPtr<SDMemTable>& in, 
-                               const casa::Vector<casa::Bool>& mask, 
-                               const casa::String& which, casa::Int row) const;
-                               
-// Bin up spectra
-   SDMemTable* bin(const SDMemTable& in, casa::Int width);
 
-  // Resample spectra
-   SDMemTable* resample(const SDMemTable& in, const casa::String& method, 
-                        casa::Float factor);
+  typedef casa::InterpolateArray1D<casa::Double,
+                                   casa::Float>::InterpolationMethod imethod;
 
-// Smooth
-  SDMemTable* smooth(const SDMemTable& in, const casa::String& kernel,
-                      casa::Float width, casa::Bool doAll);
+  typedef std::map<std::string, imethod> imap;
 
-// Flux conversion between Jansky and Kelvin
-   SDMemTable* convertFlux(const SDMemTable& in, casa::Float D, 
-                           casa::Float etaAp, 
-                           casa::Float JyPerK, casa::Bool doAll);
+  STMath(bool insitu=true);
 
-// Gain-elevation correction
-   SDMemTable* gainElevation(const SDMemTable& in, 
-                             const casa::Vector<casa::Float>& coeffs,
-                             const casa::String& fileName,
-                             const casa::String& method, 
-                             casa::Bool doAll);
+  ~STMath();
 
-// Frequency Alignment
-   SDMemTable* frequencyAlignment(const SDMemTable& in, const casa::String& refTime,
-                                   const casa::String& method, casa::Bool perFreqID);
+  /**
+   * set the @attr insitu attribute
+   * @param b
+   */
+  bool insitu() const { return insitu_;};
+  void setInsitu(bool b) { insitu_ = b; };
 
-// Opacity correction
-   SDMemTable* opacity(const SDMemTable& in, casa::Float tau, casa::Bool doAll);
+  casa::CountedPtr<Scantable>
+    average( const std::vector<casa::CountedPtr<Scantable> >& in,
+             const casa::Vector<casa::Bool>& mask,
+             const std::string& weight,
+             const std::string& avmode = "",
+             bool alignfreq = false );
 
-// Simple unary mathematical operations.  what=0 (mul) or 1 (add)
-   SDMemTable* unaryOperate(const SDMemTable& in, casa::Float offset, 
-                            casa::Bool doAll, casa::uInt what, casa::Bool tSys);
+  casa::CountedPtr<Scantable>
+    unaryOperate( const casa::CountedPtr<Scantable>& in, float val,
+                  const std::string& mode, bool tsys=false );
 
-  // Average polarizations.
-  SDMemTable* averagePol(const SDMemTable& in, 
-                         const casa::Vector<casa::Bool>& mask,
-                         const casa::String& wtStr);
+  casa::CountedPtr<Scantable> quotient( const casa::CountedPtr<Scantable>& in,
+                                        const std::string& mode = "NEAREST",
+                                        bool preserve = true );
 
-  SDMemTable* frequencySwitch(const SDMemTable& in);
-  
+  casa::CountedPtr<Scantable>
+    freqSwitch( const casa::CountedPtr<Scantable>& in );
 
-// Rotate XY phase. Value in degrees.
-   void rotateXYPhase(SDMemTable& in, casa::Float value, casa::Bool doAll);
+  std::vector<float> statistic(const casa::CountedPtr<Scantable>& in,
+                               const std::vector<bool>& mask,
+                               const std::string& which);
 
-// Rotate Q & U by operating on the raw correlations.Value in degrees.
-   void rotateLinPolPhase(SDMemTable& in, casa::Float value, casa::Bool doAll);
+  casa::CountedPtr<Scantable> bin( const casa::CountedPtr<Scantable>& in,
+                                   int width=5);
+  casa::CountedPtr<Scantable>
+    resample(const casa::CountedPtr<Scantable>& in,
+             const std::string& method, float width);
 
-  
- private:
+  casa::CountedPtr<Scantable>
+    smooth(const casa::CountedPtr<Scantable>& in, const std::string& kernel,
+                      float width);
 
-// Weighting type for time averaging
+  casa::CountedPtr<Scantable>
+    gainElevation(const casa::CountedPtr<Scantable>& in,
+                  const casa::Vector<casa::Float>& coeffs,
+                  const std::string& fileName,
+                  const std::string& method);
+  casa::CountedPtr<Scantable>
+    convertFlux(const casa::CountedPtr<Scantable>& in, float d,
+                float etaap, float jyperk);
 
-  enum WeightType {NONE=0,VAR,TSYS,TINT,TINTSYS,nWeightTypes};
+  casa::CountedPtr<Scantable> opacity(const casa::CountedPtr<Scantable>& in,
+                                      float tau);
 
-// Function to use accumulate data during time averaging
+  /// @todo frequency alignment
 
-  void accumulate(casa::Double& timeSum, casa::Double& intSum, 
-                  casa::Int& nAccum,
-                  casa::MaskedArray<casa::Float>& sum, 
-                  casa::Array<casa::Float>& sumSq,
-                  casa::Array<casa::Float>& nPts, 
-                  casa::Array<casa::Float>& tSysSum,
-                  casa::Array<casa::Float>& tSysSqSum,
-                  const casa::Array<casa::Float>& tSys,  
-                  const casa::Array<casa::Float>& nInc,
-                  const casa::Vector<casa::Bool>& mask, 
-                  casa::Double time, casa::Double interval,
-                  const std::vector<casa::CountedPtr<SDMemTable> >& in,
-                  casa::uInt iTab, casa::uInt iRow, 
-                  casa::uInt axis, casa::uInt nAxesSub,
-                  casa::Bool useMask, WeightType wtType);
+private:
+  static imethod stringToIMethod(const std::string& in);
+  static WeightType stringToWeight(const std::string& in);
 
-// Work out conversion factor for converting Jy<->K per IF per row and apply
-    void convertBrightnessUnits(SDMemTable* pTabOut, const SDMemTable& in, 
-                                casa::Bool toKelvin, casa::Float sFac, casa::Bool doAll);
+  void scaleByVector(casa::Table& in,
+                     const casa::Vector<casa::Float>& factor,
+                     bool dotsys);
 
-// Convert weight string to enum value
+  void scaleFromAsciiTable(casa::Table& in, const std::string& filename,
+                           const std::string& method,
+                           const casa::Vector<casa::Float>& xout,
+                           bool dotsys);
 
-   void convertWeightString(WeightType& wt, const casa::String& weightStr,
-                             casa::Bool listType);;
+  void scaleFromTable(casa::Table& in, const casa::Table& table,
+                      const std::string& method,
+                      const casa::Vector<casa::Float>& xout, bool dotsys);
 
-// Convert interpolation type string
-//   void convertInterpString(casa::Int& type, const casa::String& interp);
-   void convertInterpString(casa::InterpolateArray1D<casa::Double,casa::Float>::InterpolationMethod& method,  
-                             const casa::String& interp);
+  void convertBrightnessUnits(casa::CountedPtr<Scantable>& in,
+                              bool tokelvin, float cfac);
 
-// Scale data with values from an ascii Table
-   void scaleFromAsciiTable(SDMemTable* pTabOut, const SDMemTable& in, 
-                              const casa::String& fileName,
-                              const casa::String& col0, const casa::String& col1,
-                              const casa::String& methodStr, casa::Bool doAll,
-                              const casa::Vector<casa::Float>& xOut, casa::Bool doTSys);
+  casa::CountedPtr< Scantable >
+    getScantable(const casa::CountedPtr< Scantable >& in, bool droprows);
 
-// Scale data with values from a Table
-   void scaleFromTable(SDMemTable* pTabOut, const SDMemTable& in, const casa::Table& tTable,
-                         const casa::String& col0, const casa::String& col1,
-                         const casa::String& methodStr, casa::Bool doAll,
-                         const casa::Vector<casa::Float>& xOut, casa::Bool doTSys);
+  casa::MaskedArray<casa::Float>
+    maskedArray( const casa::Vector<casa::Float>& s,
+                 const casa::Vector<casa::uChar>& f );
+  casa::Vector<casa::uChar>
+    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
 
-// Scale data and optionally TSys by values in a Vector
-   void scaleByVector(SDMemTable* pTabOut, const SDMemTable& in,
-                       casa::Bool doAll, const casa::Vector<casa::Float>& factor,
-                       casa::Bool doTSys);
+  bool insitu_;
+};
 
-// Convert time String to Epoch
-   casa::MEpoch epochFromString(const casa::String& str, casa::MEpoch::Types timeRef);
-
-// Function to fill Scan Container when averaging in time
-
-  void fillSDC(SDContainer& sc, const casa::Array<casa::Bool>& mask,
-                const casa::Array<casa::Float>& data,
-                const casa::Array<casa::Float>& tSys,
-                casa::Int scanID, casa::Double timeStamp,
-                casa::Double interval, const casa::String& sourceName,
-                const casa::Vector<casa::uInt>& freqID);
-
-// Format EPoch
-   casa::String formatEpoch(const casa::MEpoch& epoch) const; 
-
-// Align in Frequency
-   SDMemTable* frequencyAlign(const SDMemTable& in,
-                              casa::MFrequency::Types system,
-                              const casa::String& timeRef,
-                              const casa::String& method, 
-                               casa::Bool perIF);
-
-// Generate frequency aligners
-   void generateFrequencyAligners(casa::PtrBlock<casa::FrequencyAligner<casa::Float>* >& a,
-                                   const SDDataDesc& dDesc,
-                                   const SDMemTable& in, casa::uInt nChan,
-                                   casa::MFrequency::Types system,
-                                   const casa::MPosition& refPos,
-                                   const casa::MEpoch& refEpoch, 
-                                   casa::Bool perFreqID);
-
-// Generate data description table(combines source and freqID);
-   void generateDataDescTable(casa::Matrix<casa::uInt>& ddIdx,
-                               SDDataDesc& dDesc,
-                               casa::uInt nIF,
-                               const SDMemTable& in,
-                               const casa::Table& tabIn,
-                               const casa::ROScalarColumn<casa::String>& srcCol,
-                               const casa::ROArrayColumn<casa::uInt>& fqIDCol,
-                               casa::Bool perFreqID);
-
-// Get row range from SDMemTable state
-  casa::Vector<casa::uInt> getRowRange(const SDMemTable& in) const;
-
-// Is row in the row range ?
-  casa::Bool rowInRange(casa::uInt i, const casa::Vector<casa::uInt>& range) const;
-
-// Set slice to cursor or all axes
-    void setCursorSlice(casa::IPosition& start, casa::IPosition& end, 
-                        casa::Bool doAll, const SDMemTable& in) const;
-
-// Function to normalize data when averaging in time
-
-  void normalize(casa::MaskedArray<casa::Float>& data,
-                  const casa::Array<casa::Float>& sumSq,
-                  const casa::Array<casa::Float>& tSysSumSq,
-                  const casa::Array<casa::Float>& nPts,
-                  casa::Double intSum,
-                 WeightType wtType, casa::Int axis, casa::Int nAxes);
-
-// Put the data and mask into the SDContainer
-   void putDataInSDC(SDContainer& sc, const casa::Array<casa::Float>& data,
-                     const casa::Array<casa::Bool>& mask);
-
-// Read ascii file into a Table
-
-  casa::Table readAsciiFile(const casa::String& fileName) const;
-
-}; // class
-
-} // namespace
-
+}
 #endif

trunk/src/Scantable.cpp

@@ -1 @@
-//#---------------------------------------------------------------------------
-//# SDMemTable.cc: A MemoryTable container for single dish integrations
-//#---------------------------------------------------------------------------
-//# Copyright (C) 2004
-//# ATNF
-//#
-//# This program is free software; you can redistribute it and/or modify it
-//# under the terms of the GNU General Public License as published by the Free
-//# Software Foundation; either version 2 of the License, or (at your option)
-//# any later version.
-//#
-//# This program is distributed in the hope that it will be useful, but
-//# WITHOUT ANY WARRANTY; without even the implied warranty of
-//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
-//# Public License for more details.
-//#
-//# You should have received a copy of the GNU General Public License along
-//# with this program; if not, write to the Free Software Foundation, Inc.,
-//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
-//#
-//# Correspondence concerning this software should be addressed as follows:
-//#        Internet email: Malte.Marquarding@csiro.au
-//#        Postal address: Malte Marquarding,
-//#                        Australia Telescope National Facility,
-//#                        P.O. Box 76,
-//#                        Epping, NSW, 2121,
-//#                        AUSTRALIA
-//#
-//# $Id:
-//#---------------------------------------------------------------------------
-
+//
+// C++ Implementation: Scantable
+//
+// Description:
+//
+//
+// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
 #include <map>
 
@@ -35 +15 @@
 #include <casa/iostream.h>
 #include <casa/iomanip.h>
+#include <casa/OS/Path.h>
+#include <casa/OS/File.h>
 #include <casa/Arrays/Array.h>
 #include <casa/Arrays/ArrayMath.h>
@@ -45 +27 @@
 #include <casa/BasicSL/Constants.h>
 #include <casa/Quanta/MVAngle.h>
+#include <casa/Containers/RecordField.h>
 
 #include <tables/Tables/TableParse.h>
@@ -52 +35 @@
 #include <tables/Tables/ScaColDesc.h>
 #include <tables/Tables/ArrColDesc.h>
+#include <tables/Tables/TableRow.h>
+#include <tables/Tables/TableVector.h>
+#include <tables/Tables/TableIter.h>
 
 #include <tables/Tables/ExprNode.h>
 #include <tables/Tables/TableRecord.h>
 #include <measures/Measures/MFrequency.h>
+#include <measures/Measures/MEpoch.h>
 #include <measures/Measures/MeasTable.h>
+#include <measures/Measures/MeasRef.h>
+#include <measures/TableMeasures/TableMeasRefDesc.h>
+#include <measures/TableMeasures/TableMeasValueDesc.h>
+#include <measures/TableMeasures/TableMeasDesc.h>
+#include <measures/TableMeasures/ScalarMeasColumn.h>
 #include <coordinates/Coordinates/CoordinateUtil.h>
 #include <casa/Quanta/MVTime.h>
 #include <casa/Quanta/MVAngle.h>
 
-#include "SDDefs.h"
-#include "SDContainer.h"
-#include "MathUtils.h"
-#include "SDPol.h"
+#include "Scantable.h"
 #include "SDAttr.h"
 
-#include "SDMemTable.h"
-
 using namespace casa;
-using namespace asap;
-
-SDMemTable::SDMemTable() :
-  IFSel_(0),
-  beamSel_(0),
-  polSel_(0)
-{
-  setup();
+
+namespace asap {
+
+Scantable::Scantable(Table::TableType ttype) :
+  type_(ttype),
+  freqTable_(ttype),
+  focusTable_(ttype),
+  weatherTable_(ttype),
+  tcalTable_(ttype),
+  moleculeTable_(ttype)
+{
+  setupMainTable();
+  table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
+  table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
+  table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
+  table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
+  table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
+  setupHistoryTable();
+  setupFitTable();
+
+  originalTable_ = table_;
   attach();
 }
 
-SDMemTable::SDMemTable(const std::string& name) :
-  IFSel_(0),
-  beamSel_(0),
-  polSel_(0)
+Scantable::Scantable(const std::string& name, Table::TableType ttype) :
+  type_(ttype),
+  freqTable_(ttype)
 {
   Table tab(name);
@@ -92 +91 @@
     throw(AipsError("Unsupported version of ASAP file."));
   }
-  table_ = tab.copyToMemoryTable("dummy");
+  if ( type_ == Table::Memory )
+    table_ = tab.copyToMemoryTable("dummy");
+  else
+    table_ = tab;
+
+  originalTable_ = table_;
   attach();
 }
 
-SDMemTable::SDMemTable(const SDMemTable& other, Bool clear)
-{
-  table_ = other.table_.copyToMemoryTable(String("dummy"));
-  // clear all rows()
+
+Scantable::Scantable( const Scantable& other, bool clear )
+{
+  // with or without data
   if (clear) {
-    table_.removeRow(this->table_.rowNumbers());
-    IFSel_= 0;
-    beamSel_= 0;
-    polSel_= 0;
-  } else {
-    IFSel_= other.IFSel_;
-    beamSel_= other.beamSel_;
-    polSel_= other.polSel_;
-  }
-
+    table_ = TableCopy::makeEmptyMemoryTable(String("dummy"), other.table_,
+                                             True);
+  } else {
+    table_ = other.table_.copyToMemoryTable(String("dummy"));
+  }
+
+  originalTable_ = table_;
   attach();
 }
 
-SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
-  IFSel_(0),
-  beamSel_(0),
-  polSel_(0)
-{
-  Table t = tableCommand(exprs,tab);
-  if (t.nrow() == 0)
-      throw(AipsError("Query unsuccessful."));
-  table_ = t.copyToMemoryTable("dummy");
-  attach();
-  renumber();
-}
-
-SDMemTable::~SDMemTable()
-{
-  //cerr << "goodbye from SDMemTable @ " << this << endl;
-}
-
-SDMemTable SDMemTable::getScan(Int scanID) const
-{
-  String cond("SELECT * from $1 WHERE SCANID == ");
-  cond += String::toString(scanID);
-  return SDMemTable(table_, cond);
-}
-
-SDMemTable &SDMemTable::operator=(const SDMemTable& other)
-{
-  if (this != &other) {
-     IFSel_= other.IFSel_;
-     beamSel_= other.beamSel_;
-     polSel_= other.polSel_;
-     table_ = other.table_.copyToMemoryTable(String("dummy"));
-     attach();
-  }
-  return *this;
-}
-
-SDMemTable SDMemTable::getSource(const std::string& source) const
-{
-  String cond("SELECT * from $1 WHERE SRCNAME == ");
-  cond += source;
-  return SDMemTable(table_, cond);
-}
-
-void SDMemTable::setup()
+Scantable::~Scantable()
+{
+  cout << "~Scantable() " << this << endl;
+}
+
+void Scantable::setupMainTable()
 {
   TableDesc td("", "1", TableDesc::Scratch);
-  td.comment() = "A SDMemTable";
+  td.comment() = "An ASAP Scantable";
   td.rwKeywordSet().define("VERSION", Int(version_));
 
+  // n Cycles
+  td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
+  // new index every nBeam x nIF x nPol
+  td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
+
+  td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
+  td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
+  td.rwKeywordSet().define("POLTYPE", String("linear"));
+  td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
+
+  td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
+  td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
+  // linear, circular, stokes [I Q U V], stokes1 [I Plinear Pangle V]
+  td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
+
   td.addColumn(ScalarColumnDesc<Double>("TIME"));
+  TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
+  TableMeasValueDesc measVal(td, "TIME");
+  TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
+  mepochCol.write(td);
+
+  td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
+
   td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
+  // Type of source (on=0, off=1, other=-1)
+  td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1)));
+  td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
+
+  //The actual Data Vectors
   td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
   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"));
-  td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
-  td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID"));
-  td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
-  td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
-  td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
-  td.addColumn(ArrayColumnDesc<Float>("TCAL"));
-  td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
-  td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
+
+  td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
+                                       IPosition(1,2),
+                                       ColumnDesc::Direct));
+  TableMeasRefDesc mdirRef(MDirection::J2000); // default
+  TableMeasValueDesc tmvdMDir(td, "DIRECTION");
+  // the TableMeasDesc gives the column a type
+  TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
+  // writing create the measure column
+  mdirCol.write(td);
+  td.addColumn(ScalarColumnDesc<Double>("AZIMUTH"));
+  td.addColumn(ScalarColumnDesc<Double>("ELEVATION"));
   td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
-  td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
-  td.addColumn(ArrayColumnDesc<Int>("FITID"));
+
+  td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
+  td.addColumn(ScalarColumnDesc<uInt>("FIT_ID"));
+
+  td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
+  td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
+
+  td.rwKeywordSet().define("OBSMODE", String(""));
 
   // 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::registerClass();
-  SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));
-  aNewTab.bindColumn("STOKES", stokesEngine);
-
-  // Create Table
   table_ = Table(aNewTab, Table::Memory, 0);
-  // add subtable
-  TableDesc tdf("", "1", TableDesc::Scratch);
-  tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
-  tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
-  tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
-  tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
-  tdf.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
-  SetupNewTable fittab("fits", tdf, Table::New);
-  Table fitTable(fittab, Table::Memory);
-  table_.rwKeywordSet().defineTable("FITS", fitTable);
-
+
+  originalTable_ = table_;
+
+}
+
+void asap::Scantable::setupHistoryTable( )
+{
   TableDesc tdh("", "1", TableDesc::Scratch);
   tdh.addColumn(ScalarColumnDesc<String>("ITEM"));
-  SetupNewTable histtab("hist", tdh, Table::New);
+  SetupNewTable histtab("history", tdh, Table::New);
   Table histTable(histtab, Table::Memory);
   table_.rwKeywordSet().defineTable("HISTORY", histTable);
 }
 
-void SDMemTable::attach()
-{
-  timeCol_.attach(table_, "TIME");
-  srcnCol_.attach(table_, "SRCNAME");
-  specCol_.attach(table_, "SPECTRA");
-  flagsCol_.attach(table_, "FLAGTRA");
-  tsCol_.attach(table_, "TSYS");
-  stokesCol_.attach(table_, "STOKES");
-  scanCol_.attach(table_, "SCANID");
-  integrCol_.attach(table_, "INTERVAL");
-  freqidCol_.attach(table_, "FREQID");
-  restfreqidCol_.attach(table_, "RESTFREQID");
-  dirCol_.attach(table_, "DIRECTION");
-  fldnCol_.attach(table_, "FIELDNAME");
-  tcaltCol_.attach(table_, "TCALTIME");
-  tcalCol_.attach(table_, "TCAL");
-  azCol_.attach(table_, "AZIMUTH");
-  elCol_.attach(table_, "ELEVATION");
-  paraCol_.attach(table_, "PARANGLE");
-  rbeamCol_.attach(table_, "REFBEAM");
-  fitCol_.attach(table_,"FITID");
-}
-
-
-std::string SDMemTable::getSourceName(Int whichRow) const
-{
-  String name;
-  srcnCol_.get(whichRow, name);
-  return name;
-}
-
-float SDMemTable::getElevation(Int whichRow) const
-{
-  float elevation;
-  elCol_.get(whichRow, elevation);
-  return elevation;
-}
-
-float SDMemTable::getAzimuth(Int whichRow) const
-{
-  float azimuth;
-  azCol_.get(whichRow, azimuth);
-  return azimuth;
-}
-
-float SDMemTable::getParAngle(Int whichRow) const
-{
-  float parangle;
-  paraCol_.get(whichRow, parangle);
-  return parangle;
-}
-
-std::string SDMemTable::getTime(Int whichRow, Bool showDate) const
-{
-  Double tm;
-  if (whichRow > -1) {
-    timeCol_.get(whichRow, tm);
-  } else {
-    table_.keywordSet().get("UTC",tm);
-  }
-  MVTime mvt(tm);
-  if (showDate)
-    mvt.setFormat(MVTime::YMD);
-  else
-    mvt.setFormat(MVTime::TIME);
-  ostringstream oss;
-  oss << mvt;
-  return String(oss);
-}
-
-double SDMemTable::getInterval(Int whichRow) const
-{
-  Double intval;
-  integrCol_.get(whichRow, intval);
-  return intval;
-}
-
-bool SDMemTable::setIF(Int whichIF)
-{
-  if ( whichIF >= 0 && whichIF < nIF()) {
-    IFSel_ = whichIF;
-    return true;
-  }
-  return false;
-}
-
-bool SDMemTable::setBeam(Int whichBeam)
-{
-  if ( whichBeam >= 0 && whichBeam < nBeam()) {
-    beamSel_ = whichBeam;
-    return true;
-  }
-  return false;
-}
-
-bool SDMemTable::setPol(Int whichPol)
-{
-  if ( whichPol >= 0 && whichPol < nPol()) {
-    polSel_ = whichPol;
-    return true;
-  }
-  return false;
-}
-
-void SDMemTable::resetCursor()
-{
-   polSel_ = 0;
-   IFSel_ = 0;
-   beamSel_ = 0;
-}
-
-std::vector<bool> SDMemTable::getMask(Int whichRow) const
-{
-
-  std::vector<bool> mask;
-
-  Array<uChar> arr;
-  flagsCol_.get(whichRow, arr);
-
-  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
-  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
-  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
-  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
-
-  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
-    bool out =!static_cast<bool>(*i);
-    mask.push_back(out);
-  }
-  return mask;
-}
-
-
-
-std::vector<float> SDMemTable::getSpectrum(Int whichRow) const
-{
-  Array<Float> arr;
-  specCol_.get(whichRow, arr);
-  return getFloatSpectrum(arr);
-}
-
-
-int SDMemTable::nStokes() const
-{
-   return stokesCol_.shape(0).nelements();        // All rows same shape
-}
-
-
-std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow,
-                                                 Bool doPol) const
-  //
-  // Gets one STokes parameter depending on cursor polSel location
-  //  doPol=False  : I,Q,U,V
-  //  doPol=True   : I,P,PA,V   ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
-  //
-{
-  AlwaysAssert(asap::nAxes==4,AipsError);
-  if (nPol()!=1 && nPol()!=2 && nPol()!=4) {
-    throw (AipsError("You must have 1,2 or 4 polarizations to get the Stokes parameters"));
-  }
-
-  // For full conversion we are only supporting linears at the moment
-
-  if (nPol() > 2) {
-    String antName;
-    table_.keywordSet().get("AntennaName", antName);
-    Instrument inst = SDAttr::convertInstrument (antName, True);
-    SDAttr sdAtt;
-    if (sdAtt.feedPolType(inst) != LINEAR) {
-      throw(AipsError("Only linear polarizations are supported"));
-    }
-  }
-
-  Array<Float> arr;
-  stokesCol_.get(whichRow, arr);
-
-  if (doPol && (polSel_==1 || polSel_==2)) {       //   Q,U --> P, P.A.
-
-    // Set current cursor location
-
-    const IPosition& shape = arr.shape();
-    IPosition start, end;
-    getCursorSlice(start, end, shape);
-
-    // Get Q and U slices
-
-    Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q");
-    Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U");
-
-    // Compute output
-
-    Array<Float> out;
-    if (polSel_==1) {                                        // P
-      out = SDPolUtil::polarizedIntensity(Q,U);
-    } else if (polSel_==2) {                                 // P.A.
-      out = SDPolUtil::positionAngle(Q,U);
-    }
-
-    // Copy to output
-
-    IPosition vecShape(1,shape(asap::ChanAxis));
-    Vector<Float> outV = out.reform(vecShape);
-    std::vector<float> stlout;
-    outV.tovector(stlout);
-    return stlout;
-
-  } else {
-    // Selects at the cursor location
-    return getFloatSpectrum(arr);
-  }
-}
-
-std::string SDMemTable::getPolarizationLabel(Bool linear, Bool stokes,
-                                              Bool linPol, Int polIdx) const
-{
-   uInt idx = polSel_;
-   if (polIdx >=0) idx = polIdx;
-   return SDPolUtil::polarizationLabel(idx, linear, stokes, linPol);
-}
-
-
-
-std::vector<float> SDMemTable::stokesToPolSpectrum(Int whichRow,
-                                                   Bool toLinear,
-                                                   Int polIdx) const
-//
-// polIdx
-//   0:3 -> RR,LL,Real(RL),Imag(RL)
-//          XX,YY,Real(XY),Image(XY)
-//
-// Gets only
-//  RR = I + V
-//  LL = I - V
-// at the moment
-//
-{
-  AlwaysAssert(asap::nAxes==4,AipsError);
-  if (nStokes()!=4) {
-     throw (AipsError("You must have 4 Stokes to convert to linear or circular"));
-  }
-//
-  Array<Float> arr, out;
-  stokesCol_.get(whichRow, arr);
-
-// Set current cursor location
-
-  const IPosition& shape = arr.shape();
-  IPosition start, end;
-  getCursorSlice(start, end, shape);
-
-// Get the slice
-
-  if (toLinear) {
-     throw(AipsError("Conversion to linears not yet supported"));
-  } else {
-    uInt selection = polSel_;
-    if (polIdx > -1) selection = polIdx;
-    Bool doRR = (selection==0);
-    if (selection>1) {
-      throw(AipsError("Only conversion to RR & LL is currently supported"));
-    }
-
-    // Get I and V slices
-    Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I");
-    Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V");
-
-    // Compute output
-    out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR);
-  }
-
-  // Copy to output
-   IPosition vecShape(1,shape(asap::ChanAxis));
-   Vector<Float> outV = out.reform(vecShape);
-   std::vector<float> stlout;
-   outV.tovector(stlout);
-//
-   return stlout;
-}
-
-
-
-
-Array<Float> SDMemTable::getStokesSpectrum(Int whichRow, Int iBeam, Int iIF) const
-{
-
-// Get data
-
-  Array<Float> arr;
-  stokesCol_.get(whichRow, arr);
-
-// Set current cursor location and overwrite polarization axis
-
-  const IPosition& shape = arr.shape();
-  IPosition start(shape.nelements(),0);
-  IPosition end(shape-1);
-  if (iBeam!=-1) {
-     start(asap::BeamAxis) = iBeam;
-     end(asap::BeamAxis) = iBeam;
-  }
-  if (iIF!=-1) {
-     start(asap::IFAxis) = iIF;
-     end(asap::IFAxis) = iIF;
-  }
-
-// Get slice
-
-  return arr(start,end);
-}
-
-
-std::vector<string> SDMemTable::getCoordInfo() const
-{
-  String un;
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-  String sunit;
-  t.keywordSet().get("UNIT",sunit);
-  String dpl;
-  t.keywordSet().get("DOPPLER",dpl);
-  if (dpl == "") dpl = "RADIO";
-  String rfrm;
-  t.keywordSet().get("REFFRAME",rfrm);
-  std::vector<string> inf;
-  inf.push_back(sunit);
-  inf.push_back(rfrm);
-  inf.push_back(dpl);
-  t.keywordSet().get("BASEREFFRAME",rfrm);
-  inf.push_back(rfrm);
-  return inf;
-}
-
-void SDMemTable::setCoordInfo(std::vector<string> theinfo)
-{
-  std::vector<string>::iterator it;
-  String un,rfrm, brfrm,dpl;
-  un = theinfo[0];              // Abcissa unit
-  rfrm = theinfo[1];            // Active (or conversion) frame
-  dpl = theinfo[2];             // Doppler
-  brfrm = theinfo[3];           // Base frame
-  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
-
-  Vector<Double> rstf;
-  t.keywordSet().get("RESTFREQS",rstf);
-
-  Bool canDo = True;
-  Unit u1("km/s");Unit u2("Hz");
-  if (Unit(un) == u1) {
-    Vector<Double> rstf;
-    t.keywordSet().get("RESTFREQS",rstf);
-    if (rstf.nelements() == 0) {
-        throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
-    }
-  } else if (Unit(un) != u2 && un != "") {
-        throw(AipsError("Unit not conformant with Spectral Coordinates"));
-  }
-  t.rwKeywordSet().define("UNIT", un);
-
-  MFrequency::Types mdr;
-  if (!MFrequency::getType(mdr, rfrm)) {
-
-    Int a,b;const uInt* c;
-    const String* valid = MFrequency::allMyTypes(a, b, c);
-    String pfix = "Please specify a legal frame type. Types are\n";
-    throw(AipsError(pfix+(*valid)));
-  } else {
-    t.rwKeywordSet().define("REFFRAME",rfrm);
-  }
-
-  MDoppler::Types dtype;
-  dpl.upcase();
-  if (!MDoppler::getType(dtype, dpl)) {
-    throw(AipsError("Doppler type unknown"));
-  } else {
-    t.rwKeywordSet().define("DOPPLER",dpl);
-  }
-
-  if (!MFrequency::getType(mdr, brfrm)) {
-     Int a,b;const uInt* c;
-     const String* valid = MFrequency::allMyTypes(a, b, c);
-     String pfix = "Please specify a legal frame type. Types are\n";
-     throw(AipsError(pfix+(*valid)));
-   } else {
-    t.rwKeywordSet().define("BASEREFFRAME",brfrm);
-   }
-}
-
-
-std::vector<double> SDMemTable::getAbcissa(Int whichRow) const
-{
-  std::vector<double> abc(nChan());
-
-  // Get header units keyword
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-  String sunit;
-  t.keywordSet().get("UNIT",sunit);
-  if (sunit == "") sunit = "pixel";
-  Unit u(sunit);
-
-  // Easy if just wanting pixels
-  if (sunit==String("pixel")) {
-    // assume channels/pixels
-    std::vector<double>::iterator it;
-    uInt i=0;
-    for (it = abc.begin(); it != abc.end(); ++it) {
-      (*it) = Double(i++);
-    }
-    return abc;
-  }
-
-  // Continue with km/s or Hz.  Get FreqIDs
-  Vector<uInt> freqIDs;
-  freqidCol_.get(whichRow, freqIDs);
-  uInt freqID = freqIDs(IFSel_);
-  restfreqidCol_.get(whichRow, freqIDs);
-  uInt restFreqID = freqIDs(IFSel_);
-
-  // Get SpectralCoordinate, set reference frame conversion,
-  // velocity conversion, and rest freq state
-
-  SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
-  Vector<Double> pixel(nChan());
-  indgen(pixel);
-
-  if (u == Unit("km/s")) {
-     Vector<Double> world;
-     spc.pixelToVelocity(world,pixel);
-     std::vector<double>::iterator it;
-     uInt i = 0;
-     for (it = abc.begin(); it != abc.end(); ++it) {
-       (*it) = world[i];
-       i++;
-     }
-  } else if (u == Unit("Hz")) {
-
-    // Set world axis units
-    Vector<String> wau(1); wau = u.getName();
-    spc.setWorldAxisUnits(wau);
-
-    std::vector<double>::iterator it;
-    Double tmp;
-    uInt i = 0;
-    for (it = abc.begin(); it != abc.end(); ++it) {
-      spc.toWorld(tmp,pixel[i]);
-      (*it) = tmp;
-      i++;
-    }
-  }
-  return abc;
-}
-
-std::string SDMemTable::getAbcissaString(Int whichRow) const
-{
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-
-  String sunit;
-  t.keywordSet().get("UNIT",sunit);
-  if (sunit == "") sunit = "pixel";
-  Unit u(sunit);
-
-  Vector<uInt> freqIDs;
-  freqidCol_.get(whichRow, freqIDs);
-  uInt freqID = freqIDs(IFSel_);
-  restfreqidCol_.get(whichRow, freqIDs);
-  uInt restFreqID = freqIDs(IFSel_);
-
-  // Get SpectralCoordinate, with frame, velocity, rest freq state set
-  SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
-
-  String s = "Channel";
-  if (u == Unit("km/s")) {
-    s = CoordinateUtil::axisLabel(spc,0,True,True,True);
-  } else if (u == Unit("Hz")) {
-    Vector<String> wau(1);wau = u.getName();
-    spc.setWorldAxisUnits(wau);
-
-    s = CoordinateUtil::axisLabel(spc,0,True,True,False);
-  }
-  return s;
-}
-
-void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow)
-{
-  Array<Float> arr;
-  specCol_.get(whichRow, arr);
-  if (spectrum.size() != arr.shape()(asap::ChanAxis)) {
-    throw(AipsError("Attempting to set spectrum with incorrect length."));
-  }
-
-  // Setup accessors
-  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));                   // Beam selection
-  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));                     // IF selection
-  ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
-  aa2.reset(aa2.begin(uInt(polSel_)));                    // Pol selection
-
-  // Iterate
-  std::vector<float>::iterator it = spectrum.begin();
-  for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
-    (*i) = Float(*it);
-    it++;
-  }
-  specCol_.put(whichRow, arr);
-}
-
-void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const
-{
-  Array<Float> arr;
-  specCol_.get(whichRow, arr);
-
-  // Iterate and extract
-  spectrum.resize(arr.shape()(3));
-  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
-  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
-  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) {
-    (*va) = (*i);
-    va++;
-  }
-}
-
-
-/*
-void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
-  Array<uChar> arr;
-  flagsCol_.get(whichRow, arr);
-  mask.resize(arr.shape()(3));
-
-  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
-  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
-  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
-  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
-
-  Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
-
-  ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask);
-  std::vector<bool> tmp;
-  tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
-                   // iterator should work on chanMask_??
-  std::vector<bool>::iterator miter;
-  miter = tmp.begin();
-
-  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
-    bool out =!static_cast<bool>(*i);
-    if (useUserMask) {
-      out = out && (*miter);
-      miter++;
-    }
-    (*va) = out;
-    va++;
-  }
-}
-*/
-
-MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
-                                                Bool toStokes) const
-{
-  // Get flags
-  Array<uChar> farr;
-  flagsCol_.get(whichRow, farr);
-
-  // Get data and convert mask to Bool
-  Array<Float> arr;
-  Array<Bool> mask;
-  if (toStokes) {
-     stokesCol_.get(whichRow, arr);
-
-     Array<Bool> tMask(farr.shape());
-     convertArray(tMask, farr);
-     mask = SDPolUtil::stokesData (tMask, True);
-  } else {
-     specCol_.get(whichRow, arr);
-     mask.resize(farr.shape());
-     convertArray(mask, farr);
-  }
-
-  return MaskedArray<Float>(arr,!mask);
-}
-
-Float SDMemTable::getTsys(Int whichRow) const
-{
-  Array<Float> arr;
-  tsCol_.get(whichRow, arr);
-  Float out;
-
-  IPosition ip(arr.shape());
-  ip(asap::BeamAxis) = beamSel_;
-  ip(asap::IFAxis) = IFSel_;
-  ip(asap::PolAxis) = polSel_;
-  ip(asap::ChanAxis)=0;               // First channel only
-
-  out = arr(ip);
-  return out;
-}
-
-MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const
-{
-  MDirection::Types mdr = getDirectionReference();
-  Array<Double> posit;
-  dirCol_.get(whichRow,posit);
-  Vector<Double> wpos(2);
-  Int rb;
-  rbeamCol_.get(whichRow,rb);
-  wpos[0] = posit(IPosition(2,beamSel_,0));
-  wpos[1] = posit(IPosition(2,beamSel_,1));
-  if (refBeam && rb > -1) {  // use refBeam instead if it exists
-    wpos[0] = posit(IPosition(2,rb,0));
-    wpos[1] = posit(IPosition(2,rb,1));
-  }
-
-  Quantum<Double> lon(wpos[0],Unit(String("rad")));
-  Quantum<Double> lat(wpos[1],Unit(String("rad")));
-  return MDirection(lon, lat, mdr);
-}
-
-MEpoch SDMemTable::getEpoch(Int whichRow) const
-{
-  MEpoch::Types met = getTimeReference();
-
-  Double obstime;
-  timeCol_.get(whichRow,obstime);
-  MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
-  return MEpoch(tm2, met);
-}
-
-MPosition SDMemTable::getAntennaPosition () const
-{
-  Vector<Double> antpos;
-  table_.keywordSet().get("AntennaPosition", antpos);
-  MVPosition mvpos(antpos(0),antpos(1),antpos(2));
-  return MPosition(mvpos);
-}
-
-
-SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const
-{
-
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-  if (freqID> t.nrow() ) {
-    throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range"));
-  }
-
-  Double rp,rv,inc;
-  String rf;
-  ROScalarColumn<Double> rpc(t, "REFPIX");
-  ROScalarColumn<Double> rvc(t, "REFVAL");
-  ROScalarColumn<Double> incc(t, "INCREMENT");
-  t.keywordSet().get("BASEREFFRAME",rf);
-
-  // Create SpectralCoordinate (units Hz)
-  MFrequency::Types mft;
-  if (!MFrequency::getType(mft, rf)) {
-    ostringstream oss;
-    pushLog("WARNING: Frequency type unknown assuming TOPO");
-    mft = MFrequency::TOPO;
-  }
-  rpc.get(freqID, rp);
-  rvc.get(freqID, rv);
-  incc.get(freqID, inc);
-
-  SpectralCoordinate spec(mft,rv,inc,rp);
-  return spec;
-}
-
-
-SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID,
-                                                     uInt restFreqID,
-                                                     uInt whichRow) const
-{
-
-  // Create basic SC
-  SpectralCoordinate spec = getSpectralCoordinate (freqID);
-
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-
-  // Set rest frequency
-  Vector<Double> restFreqIDs;
-  t.keywordSet().get("RESTFREQS",restFreqIDs);
-  if (restFreqID < restFreqIDs.nelements()) {    // Should always be true
-    spec.setRestFrequency(restFreqIDs(restFreqID),True);
-  }
-
-  // Set up frame conversion layer
-  String frm;
-  t.keywordSet().get("REFFRAME",frm);
-  if (frm == "") frm = "TOPO";
-  MFrequency::Types mtype;
-  if (!MFrequency::getType(mtype, frm)) {
-    // Should never happen
-    pushLog("WARNING: Frequency type unknown assuming TOPO");
-    mtype = MFrequency::TOPO;
-  }
-
-  // Set reference frame conversion  (requires row)
-  MDirection dir = getDirection(whichRow);
-  MEpoch epoch = getEpoch(whichRow);
-  MPosition pos = getAntennaPosition();
-
-  if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) {
-    throw(AipsError("Couldn't convert frequency frame."));
-  }
-
-  // Now velocity conversion if appropriate
-  String unitStr;
-  t.keywordSet().get("UNIT",unitStr);
-
-  String dpl;
-  t.keywordSet().get("DOPPLER",dpl);
-  MDoppler::Types dtype;
-  MDoppler::getType(dtype, dpl);
-
-  // Only set velocity unit if non-blank and non-Hz
-  if (!unitStr.empty()) {
-     Unit unitU(unitStr);
-     if (unitU==Unit("Hz")) {
-     } else {
-        spec.setVelocity(unitStr, dtype);
-     }
-  }
-
-  return spec;
-}
-
-
-Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
-                               uInt freqID) {
-  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
-  if (freqID > t.nrow() ) {
-    throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
-  }
-  ScalarColumn<Double> rpc(t, "REFPIX");
-  ScalarColumn<Double> rvc(t, "REFVAL");
-  ScalarColumn<Double> incc(t, "INCREMENT");
-
-  rpc.put(freqID, speccord.referencePixel()[0]);
-  rvc.put(freqID, speccord.referenceValue()[0]);
-  incc.put(freqID, speccord.increment()[0]);
-
-  return True;
-}
-
-Int SDMemTable::nCoordinates() const
-{
-  return table_.keywordSet().asTable("FREQUENCIES").nrow();
-}
-
-
-std::vector<double> SDMemTable::getRestFreqs() const
-{
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-  Vector<Double> tvec;
-  t.keywordSet().get("RESTFREQS",tvec);
-  std::vector<double> stlout;
-  tvec.tovector(stlout);
-  return stlout;
-}
-
-bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft)
+void Scantable::setupFitTable()
 {
   TableDesc td("", "1", TableDesc::Scratch);
-  td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
-  td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
-  td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
-  SetupNewTable aNewTab("freqs", td, Table::New);
-  Table aTable (aNewTab, Table::Memory, sdft.length());
-  ScalarColumn<Double> sc0(aTable, "REFPIX");
-  ScalarColumn<Double> sc1(aTable, "REFVAL");
-  ScalarColumn<Double> sc2(aTable, "INCREMENT");
-  for (uInt i=0; i < sdft.length(); ++i) {
-    sc0.put(i,sdft.referencePixel(i));
-    sc1.put(i,sdft.referenceValue(i));
-    sc2.put(i,sdft.increment(i));
-  }
-  String rf = sdft.refFrame();
-  if (rf.contains("TOPO")) rf = "TOPO";
-  String brf = sdft.baseRefFrame();
-  if (brf.contains("TOPO")) brf = "TOPO";
-
-  aTable.rwKeywordSet().define("BASEREFFRAME", brf);
-  aTable.rwKeywordSet().define("REFFRAME", rf);
-  aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
-  aTable.rwKeywordSet().define("UNIT", sdft.unit());
-  aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
-  Vector<Double> rfvec;
-  String rfunit;
-  sdft.restFrequencies(rfvec,rfunit);
-  Quantum<Vector<Double> > q(rfvec, rfunit);
-  rfvec.resize();
-  rfvec = q.getValue("Hz");
-  aTable.rwKeywordSet().define("RESTFREQS", rfvec);
-  table_.rwKeywordSet().defineTable("FREQUENCIES", aTable);
-  return true;
-}
-
-bool SDMemTable::putSDFitTable(const SDFitTable& sdft)
-{
-  TableDesc td("", "1", TableDesc::Scratch);
+  td.addColumn(ScalarColumnDesc<uInt>("FIT_ID"));
   td.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
   td.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
@@ -1026 +208 @@
   SetupNewTable aNewTab("fits", td, Table::New);
   Table aTable(aNewTab, Table::Memory);
-  ArrayColumn<String> sc0(aTable, "FUNCTIONS");
-  ArrayColumn<Int> sc1(aTable, "COMPONENTS");
-  ArrayColumn<Double> sc2(aTable, "PARAMETERS");
-  ArrayColumn<Bool> sc3(aTable, "PARMASK");
-  ArrayColumn<String> sc4(aTable, "FRAMEINFO");
-  for (uInt i; i<sdft.length(); ++i) {
-    const Vector<Double>& parms = sdft.getParameters(i);
-    const Vector<Bool>& parmask = sdft.getParameterMask(i);
-    const Vector<String>& funcs = sdft.getFunctions(i);
-    const Vector<Int>& comps = sdft.getComponents(i);
-    const Vector<String>& finfo = sdft.getFrameInfo(i);
-    sc0.put(i,funcs);
-    sc1.put(i,comps);
-    sc3.put(i,parmask);
-    sc2.put(i,parms);
-    sc4.put(i,finfo);
-  }
   table_.rwKeywordSet().defineTable("FITS", aTable);
-  return true;
-}
-
-SDFitTable SDMemTable::getSDFitTable() const
-{
-  const Table& t = table_.keywordSet().asTable("FITS");
-  Vector<Double> parms;
-  Vector<Bool> parmask;
-  Vector<String> funcs;
-  Vector<String> finfo;
-  Vector<Int> comps;
-  ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
-  ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
-  ROArrayColumn<Int> compsCol(t, "COMPONENTS");
-  ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
-  ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
-  uInt n = t.nrow();
-  SDFitTable sdft;
-  for (uInt i=0; i<n; ++i) {
-    parmaskCol.get(i, parmask);
-    parmsCol.get(i, parms);
-    funcsCol.get(i, funcs);
-    compsCol.get(i, comps);
-    finfoCol.get(i, finfo);
-    sdft.addFit(parms, parmask, funcs, comps, finfo);
-  }
-  return sdft;
-}
-
-SDFitTable SDMemTable::getSDFitTable(uInt whichRow) const {
-  const Table& t = table_.keywordSet().asTable("FITS");
-  if (t.nrow() == 0 || whichRow >= t.nrow()) return SDFitTable();
-  Array<Int> fitid;
-  fitCol_.get(whichRow, fitid);
-  if (fitid.nelements() == 0) return SDFitTable();
-
-  IPosition shp = fitid.shape();
-  IPosition start(4, beamSel_, IFSel_, polSel_,0);
-  IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
-
-  // reform the output array slice to be of dim=1
-  Vector<Int> tmp = (fitid(start, end)).reform(IPosition(1,shp[3]));
-
-  Vector<Double> parms;
-  Vector<Bool> parmask;
-  Vector<String> funcs;
-  Vector<String> finfo;
-  Vector<Int> comps;
-  ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
-  ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
-  ROArrayColumn<Int> compsCol(t, "COMPONENTS");
-  ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
-  ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
-  SDFitTable sdft;
-  Int k=-1;
-  for (uInt i=0; i< tmp.nelements(); ++i) {
-    k = tmp[i];
-    if ( k > -1 && k < t.nrow() ) {
-      parms.resize();
-      parmsCol.get(k, parms);
-      parmask.resize();
-      parmaskCol.get(k, parmask);
-      funcs.resize();
-      funcsCol.get(k, funcs);
-      comps.resize();
-      compsCol.get(k, comps);
-      finfo.resize();
-      finfoCol.get(k, finfo);
-      sdft.addFit(parms, parmask, funcs, comps, finfo);
-    }
-  }
-  return sdft;
-}
-
-void SDMemTable::addFit(uInt whichRow,
-                        const Vector<Double>& p, const Vector<Bool>& m,
-                        const Vector<String>& f, const Vector<Int>& c)
-{
-  if (whichRow >= nRow()) {
-    throw(AipsError("Specified row out of range"));
-  }
-  Table t = table_.keywordSet().asTable("FITS");
-  uInt nrow = t.nrow();
-  t.addRow();
-  ArrayColumn<Double> parmsCol(t, "PARAMETERS");
-  ArrayColumn<Bool> parmaskCol(t, "PARMASK");
-  ArrayColumn<Int> compsCol(t, "COMPONENTS");
-  ArrayColumn<String> funcsCol(t, "FUNCTIONS");
-  ArrayColumn<String> finfoCol(t, "FRAMEINFO");
-  parmsCol.put(nrow, p);
-  parmaskCol.put(nrow, m);
-  compsCol.put(nrow, c);
-  funcsCol.put(nrow, f);
-  Vector<String> fi = mathutil::toVectorString(getCoordInfo());
-  finfoCol.put(nrow, fi);
-
-  Array<Int> fitarr;
-  fitCol_.get(whichRow, fitarr);
-
-  Array<Int> newarr;               // The new Array containing the fitid
-  Int pos =-1;                     // The fitid position in the array
-  if ( fitarr.nelements() == 0 ) { // no fits at all in this row
-    Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1));
-    arr = -1;
-    newarr.reference(arr);
-    pos = 0;
-  } else {
-    IPosition shp = fitarr.shape();
-    IPosition start(4, beamSel_, IFSel_, polSel_,0);
-    IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
-    // reform the output array slice to be of dim=1
-    Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3]));
-    const Vector<Int>& fits = tmp;
-    VectorSTLIterator<Int> it(fits);
-    Int i = 0;
-    while (it != fits.end()) {
-      if (*it == -1) {
-        pos = i;
-        break;
-      }
-      ++i;
-      ++it;
-    };
-  }
-  if (pos == -1) {
-      mathutil::extendLastArrayAxis(newarr,fitarr, -1);
-      pos = fitarr.shape()[3];     // the new element position
-  } else {
-    if (fitarr.nelements() > 0)
-      newarr = fitarr;
-  }
-  newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow);
-  fitCol_.put(whichRow, newarr);
-
-}
-
-SDFrequencyTable SDMemTable::getSDFreqTable() const
-{
-  const Table& t = table_.keywordSet().asTable("FREQUENCIES");
-  SDFrequencyTable sdft;
-
-  // Add refpix/refval/incr.  What are the units ? Hz I suppose
-  // but it's nowhere described...
-  Vector<Double> refPix, refVal, incr;
-  ScalarColumn<Double> refPixCol(t, "REFPIX");
-  ScalarColumn<Double> refValCol(t, "REFVAL");
-  ScalarColumn<Double> incrCol(t, "INCREMENT");
-  refPix = refPixCol.getColumn();
-  refVal = refValCol.getColumn();
-  incr = incrCol.getColumn();
-
-  uInt n = refPix.nelements();
-  for (uInt i=0; i<n; i++) {
-     sdft.addFrequency(refPix[i], refVal[i], incr[i]);
-  }
-
-  // Frequency reference frame.  I don't know if this
-  // is the correct frame.  It might be 'REFFRAME'
-  // rather than 'BASEREFFRAME' ?
-  String frame;
-  t.keywordSet().get("REFFRAME",frame);
-  sdft.setRefFrame(frame);
-  t.keywordSet().get("BASEREFFRAME",frame);
-  sdft.setBaseRefFrame(frame);
-
-  // Equinox
-  Float equinox;
-  t.keywordSet().get("EQUINOX", equinox);
-  sdft.setEquinox(equinox);
-
-  String unit;
-  t.keywordSet().get("UNIT", unit);
-  sdft.setUnit(unit);
-
-  // Rest Frequency
-  Vector<Double> restFreqs;
-  t.keywordSet().get("RESTFREQS", restFreqs);
-  for (uInt i=0; i<restFreqs.nelements(); i++) {
-     sdft.addRestFrequency(restFreqs[i]);
-  }
-  sdft.setRestFrequencyUnit(String("Hz"));
-
-  return sdft;
-}
-
-bool SDMemTable::putSDContainer(const SDContainer& sdc)
-{
-  uInt rno = table_.nrow();
-  table_.addRow();
-
-  timeCol_.put(rno, sdc.timestamp);
-  srcnCol_.put(rno, sdc.sourcename);
-  fldnCol_.put(rno, sdc.fieldname);
-  specCol_.put(rno, sdc.getSpectrum());
-  flagsCol_.put(rno, sdc.getFlags());
-  tsCol_.put(rno, sdc.getTsys());
-  scanCol_.put(rno, sdc.scanid);
-  integrCol_.put(rno, sdc.interval);
-  freqidCol_.put(rno, sdc.getFreqMap());
-  restfreqidCol_.put(rno, sdc.getRestFreqMap());
-  dirCol_.put(rno, sdc.getDirection());
-  rbeamCol_.put(rno, sdc.refbeam);
-  tcalCol_.put(rno, sdc.tcal);
-  tcaltCol_.put(rno, sdc.tcaltime);
-  azCol_.put(rno, sdc.azimuth);
-  elCol_.put(rno, sdc.elevation);
-  paraCol_.put(rno, sdc.parangle);
-  fitCol_.put(rno, sdc.getFitMap());
-  return true;
-}
-
-SDContainer SDMemTable::getSDContainer(uInt whichRow) const
-{
-  SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
-  timeCol_.get(whichRow, sdc.timestamp);
-  srcnCol_.get(whichRow, sdc.sourcename);
-  integrCol_.get(whichRow, sdc.interval);
-  scanCol_.get(whichRow, sdc.scanid);
-  fldnCol_.get(whichRow, sdc.fieldname);
-  rbeamCol_.get(whichRow, sdc.refbeam);
-  azCol_.get(whichRow, sdc.azimuth);
-  elCol_.get(whichRow, sdc.elevation);
-  paraCol_.get(whichRow, sdc.parangle);
-  Vector<Float> tc;
-  tcalCol_.get(whichRow, tc);
-  sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
-  tcaltCol_.get(whichRow, sdc.tcaltime);
-
-  Array<Float> spectrum;
-  Array<Float> tsys;
-  Array<uChar> flagtrum;
-  Vector<uInt> fmap;
-  Array<Double> direction;
-  Array<Int> fits;
-
-  specCol_.get(whichRow, spectrum);
-  sdc.putSpectrum(spectrum);
-  flagsCol_.get(whichRow, flagtrum);
-  sdc.putFlags(flagtrum);
-  tsCol_.get(whichRow, tsys);
-  sdc.putTsys(tsys);
-  freqidCol_.get(whichRow, fmap);
-  sdc.putFreqMap(fmap);
-  restfreqidCol_.get(whichRow, fmap);
-  sdc.putRestFreqMap(fmap);
-  dirCol_.get(whichRow, direction);
-  sdc.putDirection(direction);
-  fitCol_.get(whichRow, fits);
-  sdc.putFitMap(fits);
-  return sdc;
-}
-
-bool SDMemTable::putSDHeader(const SDHeader& sdh)
+}
+
+void Scantable::attach()
+{
+  historyTable_ = table_.keywordSet().asTable("HISTORY");
+  fitTable_ = table_.keywordSet().asTable("FITS");
+
+  timeCol_.attach(table_, "TIME");
+  srcnCol_.attach(table_, "SRCNAME");
+  specCol_.attach(table_, "SPECTRA");
+  flagsCol_.attach(table_, "FLAGTRA");
+  tsCol_.attach(table_, "TSYS");
+  cycleCol_.attach(table_,"CYCLENO");
+  scanCol_.attach(table_, "SCANNO");
+  beamCol_.attach(table_, "BEAMNO");
+  integrCol_.attach(table_, "INTERVAL");
+  azCol_.attach(table_, "AZIMUTH");
+  elCol_.attach(table_, "ELEVATION");
+  dirCol_.attach(table_, "DIRECTION");
+  paraCol_.attach(table_, "PARANGLE");
+  fldnCol_.attach(table_, "FIELDNAME");
+  rbeamCol_.attach(table_, "REFBEAMNO");
+
+  mfitidCol_.attach(table_,"FIT_ID");
+  fitidCol_.attach(fitTable_,"FIT_ID");
+
+  mfreqidCol_.attach(table_, "FREQ_ID");
+
+  mtcalidCol_.attach(table_, "TCAL_ID");
+
+  mfocusidCol_.attach(table_, "FOCUS_ID");
+
+  mmolidCol_.attach(table_, "MOLECULE_ID");
+}
+
+void Scantable::putSDHeader(const SDHeader& sdh)
 {
   table_.rwKeywordSet().define("nIF", sdh.nif);
@@ -1314 +262 @@
   table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
   table_.rwKeywordSet().define("Epoch", sdh.epoch);
-  return true;
-}
-
-SDHeader SDMemTable::getSDHeader() const
+}
+
+SDHeader Scantable::getSDHeader() const
 {
   SDHeader sdh;
@@ -1338 +285 @@
   return sdh;
 }
-void SDMemTable::makePersistent(const std::string& filename)
-{
-  table_.deepCopy(filename,Table::New);
-
-}
-
-Int SDMemTable::nScan() const {
-  Int n = 0;
-  Int previous = -1;Int current=0;
+
+int Scantable::rowToScanIndex( int therow )
+{
+  int therealrow = -1;
+
+  return therealrow;
+}
+
+int Scantable::nScan() const {
+  int n = 0;
+  int previous = -1; int current = 0;
   for (uInt i=0; i< scanCol_.nrow();i++) {
     scanCol_.getScalar(i,current);
@@ -1357 +306 @@
 }
 
-String SDMemTable::formatSec(Double x) const
+std::string Scantable::formatSec(Double x) const
 {
   Double xcop = x;
@@ -1374 +323 @@
 };
 
-String SDMemTable::formatDirection(const MDirection& md) const
+std::string Scantable::formatDirection(const MDirection& md) const
 {
   Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
@@ -1387 +336 @@
 
 
-std::string SDMemTable::getFluxUnit() const
+std::string Scantable::getFluxUnit() const
 {
   String tmp;
@@ -1394 +343 @@
 }
 
-void SDMemTable::setFluxUnit(const std::string& unit)
+void Scantable::setFluxUnit(const std::string& unit)
 {
   String tmp(unit);
@@ -1405 +354 @@
 }
 
-
-void SDMemTable::setInstrument(const std::string& name)
-{
-  Bool throwIt = True;
+void Scantable::setInstrument(const std::string& name)
+{
+  bool throwIt = true;
   Instrument ins = SDAttr::convertInstrument(name, throwIt);
   String nameU(name);
@@ -1415 +363 @@
 }
 
-std::string SDMemTable::summary(bool verbose) const  {
-
-  // Format header info
-  ostringstream oss;
-  oss << endl;
-  oss << "--------------------------------------------------------------------------------" << endl;
-  oss << " Scan Table Summary" << endl;
-  oss << "--------------------------------------------------------------------------------" << endl;
-  oss.flags(std::ios_base::left);
-  oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
-      << setw(15) << "IFs:" << setw(4) << nIF() << endl
-      << setw(15) << "Polarisations:" << setw(4) << nPol() << endl
-      << setw(15) << "Channels:"  << setw(4) << nChan() << endl;
-  oss << endl;
-  String tmp;
-  table_.keywordSet().get("Observer", tmp);
-  oss << setw(15) << "Observer:" << tmp << endl;
-  oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl;
-  table_.keywordSet().get("Project", tmp);
-  oss << setw(15) << "Project:" << tmp << endl;
-  table_.keywordSet().get("Obstype", tmp);
-  oss << setw(15) << "Obs. Type:" << tmp << endl;
-  table_.keywordSet().get("AntennaName", tmp);
-  oss << setw(15) << "Antenna Name:" << tmp << endl;
-  table_.keywordSet().get("FluxUnit", tmp);
-  oss << setw(15) << "Flux Unit:" << tmp << endl;
-  Table t = table_.keywordSet().asTable("FREQUENCIES");
-  Vector<Double> vec;
-  t.keywordSet().get("RESTFREQS",vec);
-  oss << setw(15) << "Rest Freqs:";
-  if (vec.nelements() > 0) {
-      oss << setprecision(10) << vec << " [Hz]" << endl;
-  } else {
-      oss << "None set" << endl;
-  }
-  oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
-  oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
-      << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
-  oss << endl;
-
-  String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")";
-  oss << setw(5) << "Scan"
-      << setw(15) << "Source"
-      << setw(24) << dirtype
-      << setw(10) << "Time"
-      << setw(18) << "Integration"
-      << setw(7) << "FreqIDs" << endl;
-  oss << "--------------------------------------------------------------------------------" << endl;
-
-  // Generate list of scan start and end integrations
-  Vector<Int> scanIDs = scanCol_.getColumn();
-  Vector<uInt> startInt, endInt;
-  mathutil::scanBoundaries(startInt, endInt, scanIDs);
-
-  const uInt nScans = startInt.nelements();
-  String name;
-  Vector<uInt> freqIDs, listFQ;
-  Vector<uInt> restFreqIDs, listRestFQ;
-  uInt nInt;
-
-  for (uInt i=0; i<nScans; i++) {
-    // Get things from first integration of scan
-    String time = getTime(startInt(i),False);
-    String tInt = formatSec(Double(getInterval(startInt(i))));
-    String posit = formatDirection(getDirection(startInt(i),True));
-    srcnCol_.getScalar(startInt(i),name);
-
-    // Find all the FreqIDs in this scan
-    listFQ.resize(0);
-    listRestFQ.resize(0);
-    for (uInt j=startInt(i); j<endInt(i)+1; j++) {
-      freqidCol_.get(j, freqIDs);
-      for (uInt k=0; k<freqIDs.nelements(); k++) {
-        mathutil::addEntry(listFQ, freqIDs(k));
-      }
-//
-      restfreqidCol_.get(j, restFreqIDs);
-      for (uInt k=0; k<restFreqIDs.nelements(); k++) {
-        mathutil::addEntry(listRestFQ, restFreqIDs(k));
-      }
-    }
-
-    nInt = endInt(i) - startInt(i) + 1;
-    oss << setw(3) << std::right << i << std::left << setw(2) << "  "
-        << setw(15) << name
-        << setw(24) << posit
-        << setw(10) << time
-        << setw(3) << std::right << nInt  << setw(3) << " x " << std::left
-        << setw(6) <<  tInt
-        << " " << listFQ << " " << listRestFQ << endl;
-  }
-  oss << endl;
-  oss << "Table contains " << table_.nrow() << " integration(s) in "
-      << nScans << " scan(s)." << endl;
-
-  // Frequency Table
-  if (verbose) {
-    std::vector<string> info = getCoordInfo();
-    SDFrequencyTable sdft = getSDFreqTable();
-    oss << endl << endl;
-    oss << "FreqID  Frame   RefFreq(Hz)     RefPix   Increment(Hz)" << endl;
-    oss << "--------------------------------------------------------------------------------" << endl;
-    for (uInt i=0; i<sdft.length(); i++) {
-      oss << setw(8) << i << setw(8)
-          << info[3] << setw(16) << setprecision(8)
-          << sdft.referenceValue(i) << setw(10)
-          << sdft.referencePixel(i) << setw(12)
-          << sdft.increment(i) << endl;
-    }
-    oss << "--------------------------------------------------------------------------------" << endl;
-  }
-  return String(oss);
-}
-/*
-std::string SDMemTable::scanSummary(const std::vector<int>& whichScans) {
-  ostringstream oss;
-  Vector<Int> scanIDs = scanCol_.getColumn();
-  Vector<uInt> startInt, endInt;
-  mathutil::scanBoundaries(startInt, endInt, scanIDs);
-  const uInt nScans = startInt.nelements();
-  std::vector<int>::const_iterator it(whichScans);
-  return String(oss);
-}
-*/
-Int SDMemTable::nBeam() const
-{
-  Int n;
-  table_.keywordSet().get("nBeam",n);
-  return n;
-}
-
-Int SDMemTable::nIF() const {
-  Int n;
-  table_.keywordSet().get("nIF",n);
-  return n;
-}
-
-Int SDMemTable::nPol() const {
-  Int n;
-  table_.keywordSet().get("nPol",n);
-  return n;
-}
-
-Int SDMemTable::nChan() const {
-  Int n;
-  table_.keywordSet().get("nChan",n);
-  return n;
-}
-
-Table SDMemTable::getHistoryTable() const
+MPosition Scantable::getAntennaPosition () const
+{
+  Vector<Double> antpos;
+  table_.keywordSet().get("AntennaPosition", antpos);
+  MVPosition mvpos(antpos(0),antpos(1),antpos(2));
+  return MPosition(mvpos);
+}
+
+void Scantable::makePersistent(const std::string& filename)
+{
+  String inname(filename);
+  Path path(inname);
+  inname = path.expandedName();
+  table_.deepCopy(inname, Table::New);
+}
+
+int asap::Scantable::nbeam( int scanno ) const
+{
+  if ( scanno < 0 ) {
+    Int n;
+    table_.keywordSet().get("nBeam",n);
+    return int(n);
+  } else {
+    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
+    Table tab = table_(table_.col("SCANNO") == scanno
+                       && table_.col("POLNO") == 0
+                       && table_.col("IFNO") == 0
+                       && table_.col("CYCLENO") == 0 );
+    ROTableVector<uInt> v(tab, "BEAMNO");
+    return int(v.nelements());
+  }
+  return 0;
+}
+
+int asap::Scantable::nif( int scanno ) const
+{
+  if ( scanno < 0 ) {
+    Int n;
+    table_.keywordSet().get("nIF",n);
+    return int(n);
+  } else {
+    // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
+    Table tab = table_(table_.col("SCANNO") == scanno
+                       && table_.col("POLNO") == 0
+                       && table_.col("BEAMNO") == 0
+                       && table_.col("CYCLENO") == 0 );
+    ROTableVector<uInt> v(tab, "IFNO");
+    return int(v.nelements());
+  }
+  return 0;
+}
+
+int asap::Scantable::npol( int scanno ) const
+{
+  if ( scanno < 0 ) {
+    Int n;
+    table_.keywordSet().get("nPol",n);
+    return n;
+  } else {
+    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
+    Table tab = table_(table_.col("SCANNO") == scanno
+                       && table_.col("IFNO") == 0
+                       && table_.col("BEAMNO") == 0
+                       && table_.col("CYCLENO") == 0 );
+    ROTableVector<uInt> v(tab, "POLNO");
+    return int(v.nelements());
+  }
+  return 0;
+}
+
+int asap::Scantable::nrow( int scanno ) const
+{
+  if ( scanno < 0 ) {
+    return int(table_.nrow());
+  } else {
+    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
+    Table tab = table_(table_.col("SCANNO") == scanno
+                       && table_.col("BEAMNO") == 0
+                       && table_.col("IFNO") == 0
+                       && table_.col("POLNO") == 0 );
+    return int(tab.nrow());
+  }
+  return 0;
+}
+
+
+int asap::Scantable::nchan( int scanno, int ifno ) const
+{
+  if ( scanno < 0 && ifno < 0 ) {
+    Int n;
+    table_.keywordSet().get("nChan",n);
+    return int(n);
+  } else {
+    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
+    Table tab = table_(table_.col("SCANNO") == scanno
+                       && table_.col("IFNO") == ifno
+                       && table_.col("BEAMNO") == 0
+                       && table_.col("POLNO") == 0
+                       && table_.col("CYCLENO") == 0 );
+    ROArrayColumn<Float> v(tab, "SPECTRA");
+    cout << v.shape(0) << endl;
+    return 0;
+  }
+  return 0;
+}
+
+Table Scantable::getHistoryTable() const
 {
   return table_.keywordSet().asTable("HISTORY");
 }
 
-void SDMemTable::appendToHistoryTable(const Table& otherHist)
+void Scantable::appendToHistoryTable(const Table& otherHist)
 {
   Table t = table_.rwKeywordSet().asTable("HISTORY");
-  const String sep = "--------------------------------------------------------------------------------";
-  addHistory(sep);
+
+  addHistory(asap::SEPERATOR);
   TableCopy::copyRows(t, otherHist, t.nrow(), 0, otherHist.nrow());
-  addHistory(sep);
-}
-
-void SDMemTable::addHistory(const std::string& hist)
+  addHistory(asap::SEPERATOR);
+}
+
+void Scantable::addHistory(const std::string& hist)
 {
   Table t = table_.rwKeywordSet().asTable("HISTORY");
@@ -1587 +492 @@
 }
 
-std::vector<std::string> SDMemTable::getHistory() const
+std::vector<std::string> Scantable::getHistory() const
 {
   Vector<String> history;
@@ -1602 +507 @@
 }
 
-
-/*
-  void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
-
-  std::vector<int>::iterator it;
-  ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_);
-  for (it = whichChans.begin(); it != whichChans.end(); it++) {
-    j.reset(j.begin(uInt(*it)));
-    for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) {
-      for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) {
-        for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii);
-             iii != iii.end(); ++iii) {
-          (*iii) =
-        }
-      }
-    }
-  }
-
-}
-*/
-void SDMemTable::flag(int whichRow)
-  {
-  Array<uChar> arr;
-  flagsCol_.get(whichRow, arr);
-
-  ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
-  ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
-  ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
-  aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
-
-  for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
-    (*i) = uChar(True);
-  }
-
-  flagsCol_.put(whichRow, arr);
-}
-
-MDirection::Types SDMemTable::getDirectionReference() const
-{
-  Float eq;
-  table_.keywordSet().get("Equinox",eq);
-  std::map<float,string> mp;
-  mp[2000.0] = "J2000";
-  mp[1950.0] = "B1950";
-  MDirection::Types mdr;
-  if (!MDirection::getType(mdr, mp[eq])) {
-    mdr = MDirection::J2000;
-    pushLog("WARNING: Unknown equinox using J2000");
-  }
-
-  return mdr;
-}
-
-MEpoch::Types SDMemTable::getTimeReference() const
-{
-  MEpoch::Types met;
-  String ep;
-  table_.keywordSet().get("Epoch",ep);
-  if (!MEpoch::getType(met, ep)) {
-     pushLog("WARNING: Epoch type unknown - using UTC");
-    met = MEpoch::UTC;
-  }
-
-  return met;
-}
-
-
-Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn,
-                              const String& sUnit,
-                              const vector<string>& lines,
-                              const String& source,
-                              Int whichIF)
-{
-   const Int nIFs = nIF();
-   if (whichIF>=nIFs) {
-      throw(AipsError("Illegal IF index"));
-   }
-
-   // Find vector of restfrequencies
-   // Double takes precedence over String
-   Unit unit;
-   Vector<Double> restFreqs;
-   if (restFreqsIn.nelements()>0) {
-      restFreqs.resize(restFreqsIn.nelements());
-      restFreqs = restFreqsIn;
-      unit = Unit(sUnit);
-   } else if (lines.size()>0) {
-      const uInt nLines = lines.size();
-      unit = Unit("Hz");
-      restFreqs.resize(nLines);
-      MFrequency lineFreq;
-      for (uInt i=0; i<nLines; i++) {
-         String tS(lines[i]);
-         tS.upcase();
-         if (MeasTable::Line(lineFreq, tS)) {
-            restFreqs[i] = lineFreq.getValue().getValue();          // Hz
-         } else {
-            String s = String(lines[i]) +
-              String(" is an unrecognized spectral line");
-            throw(AipsError(s));
-         }
-      }
-   } else {
-      throw(AipsError("You have not specified any rest frequencies or lines"));
-   }
-
-   // If multiple restfreqs, must be length nIF. In this
-   // case we will just replace the rest frequencies
-   // We can't disinguish scalar and vector of length 1
-   const uInt nRestFreqs = restFreqs.nelements();
-   Int idx = -1;
-   SDFrequencyTable sdft = getSDFreqTable();
-
-   ostringstream oss;
-   if (nRestFreqs>1) {
-     // Replace restFreqs, one per IF
-     if (nRestFreqs != nIFs) {
-       throw (AipsError("Number of rest frequencies must be equal to the number of IFs"));
-     }
-     ostringstream oss;
-     oss << "Replaced rest frequencies, one per IF, with given list : " << restFreqs;
-     sdft.deleteRestFrequencies();
-     for (uInt i=0; i<nRestFreqs; i++) {
-       Quantum<Double> rf(restFreqs[i], unit);
-       sdft.addRestFrequency(rf.getValue("Hz"));
-     }
-   } else {
-
-     // Add new rest freq
-      Quantum<Double> rf(restFreqs[0], unit);
-      idx = sdft.addRestFrequency(rf.getValue("Hz"));
-      if (whichIF>=0) {
-         oss << "Selected given rest frequency (" << restFreqs[0] << ") for IF " << whichIF << endl;
-      } else {
-         oss << "Selected given rest frequency (" << restFreqs[0] << ") for all IFs" << endl;
-      }
-   }
-   pushLog(String(oss));
-   // Replace
-   Bool empty = source.empty();
-   Bool ok = False;
-   if (putSDFreqTable(sdft)) {
-      const uInt nRow = table_.nrow();
-      String srcName;
-      Vector<uInt> restFreqIDs;
-      for (uInt i=0; i<nRow; i++) {
-         srcnCol_.get(i, srcName);
-         restfreqidCol_.get(i,restFreqIDs);
-         if (idx==-1) {
-           // Replace vector of restFreqs; one per IF.
-           // No selection possible
-            for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i;
-         } else {
-           // Set RestFreqID for selected data
-            if (empty || source==srcName) {
-               if (whichIF<0) {
-                  restFreqIDs = idx;
-               } else {
-                  restFreqIDs[whichIF] = idx;
-               }
-            }
-         }
-         restfreqidCol_.put(i,restFreqIDs);
-      }
-      ok = True;
-   } else {
-     ok = False;
-   }
-
-   return ok;
-}
-
-std::string SDMemTable::spectralLines() const
-{
-   Vector<String> lines = MeasTable::Lines();
-   MFrequency lineFreq;
-   Double freq;
-   ostringstream oss;
-
-   oss.flags(std::ios_base::left);
-   oss << "Line      Frequency (Hz)" << endl;
-   oss << "-----------------------" << endl;
-   for (uInt i=0; i<lines.nelements(); i++) {
-     MeasTable::Line(lineFreq, lines[i]);
-     freq = lineFreq.getValue().getValue();          // Hz
-     oss << setw(11) << lines[i] << setprecision(10) << freq << endl;
-   }
-   return String(oss);
-}
-
-void SDMemTable::renumber()
-{
-  uInt nRow = scanCol_.nrow();
-  Int newscanid = 0;
-  Int cIdx;// the current scanid
-  // get the first scanid
-  scanCol_.getScalar(0,cIdx);
-  Int pIdx = cIdx;// the scanid of the previous row
-  for (uInt i=0; i<nRow;++i) {
-    scanCol_.getScalar(i,cIdx);
-    if (pIdx == cIdx) {
-      // renumber
-      scanCol_.put(i,newscanid);
-    } else {
-      ++newscanid;
-      pIdx = cIdx;   // store scanid
-      --i;           // don't increment next loop
-    }
-  }
-}
-
-
-void SDMemTable::getCursorSlice(IPosition& start, IPosition& end,
-                                const IPosition& shape) const
-{
-  const uInt nDim = shape.nelements();
-  start.resize(nDim);
-  end.resize(nDim);
-
-  start(asap::BeamAxis) = beamSel_;
-  end(asap::BeamAxis) = beamSel_;
-  start(asap::IFAxis) = IFSel_;
-  end(asap::IFAxis) = IFSel_;
-
-  start(asap::PolAxis) = polSel_;
-  end(asap::PolAxis) = polSel_;
-
-  start(asap::ChanAxis) = 0;
-  end(asap::ChanAxis) = shape(asap::ChanAxis) - 1;
-}
-
-
-std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const
-  // Get spectrum at cursor location
-{
-
-  // Setup accessors
-  ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
-  aa0.reset(aa0.begin(uInt(beamSel_)));                    // Beam selection
-
-  ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
-  aa1.reset(aa1.begin(uInt(IFSel_)));                      // IF selection
-
-  ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
-  aa2.reset(aa2.begin(uInt(polSel_)));                     // Pol selection
-
-  std::vector<float> spectrum;
-  for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
-    spectrum.push_back(*i);
-  }
-  return spectrum;
-}
-
-void SDMemTable::calculateAZEL()
+std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
+{
+  MVTime mvt(me.getValue());
+  if (showdate)
+    mvt.setFormat(MVTime::YMD);
+  else
+    mvt.setFormat(MVTime::TIME);
+  ostringstream oss;
+  oss << mvt;
+  return String(oss);
+}
+
+void Scantable::calculateAZEL()
 {
   MPosition mp = getAntennaPosition();
+  MEpoch::ROScalarColumn timeCol(table_, "TIME");
   ostringstream oss;
   oss << "Computed azimuth/elevation using " << endl
       << mp << endl;
-  for (uInt i=0; i<nRow();++i) {
-    MEpoch me = getEpoch(i);
-    MDirection md = getDirection(i,False);
-    oss  << " Time: " << getTime(i,False) << " Direction: " << formatDirection(md)
+  for (uInt i=0; i<nrow(); ++i) {
+    MEpoch me = timeCol(i);
+    MDirection md = dirCol_(i);
+    dirCol_.get(i,md);
+    oss  << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
          << endl << "     => ";
     MeasFrame frame(mp, me);
@@ -1880 +544 @@
   pushLog(String(oss));
 }
+
+double Scantable::getInterval(int whichrow) const
+{
+  if (whichrow < 0) return 0.0;
+  Double intval;
+  integrCol_.get(Int(whichrow), intval);
+  return intval;
+}
+
+std::vector<bool> Scantable::getMask(int whichrow) const
+{
+  Vector<uChar> flags;
+  flagsCol_.get(uInt(whichrow), flags);
+  Vector<Bool> bflag(flags.shape());
+  convertArray(bflag, flags);
+  bflag = !bflag;
+  std::vector<bool> mask;
+  bflag.tovector(mask);
+  return mask;
+}
+
+std::vector<float> Scantable::getSpectrum(int whichrow) const
+{
+  Vector<Float> arr;
+  specCol_.get(whichrow, arr);
+  std::vector<float> out;
+  arr.tovector(out);
+  return out;
+}
+
+String Scantable::generateName()
+{
+  return (File::newUniqueName("./","temp")).baseName();
+}
+
+const casa::Table& Scantable::table( ) const
+{
+  return table_;
+}
+
+casa::Table& Scantable::table( )
+{
+  return table_;
+}
+
+std::string Scantable::getPolarizationLabel(bool linear, bool stokes,
+                                            bool linpol, int polidx) const
+{
+  uInt idx = 0;
+  if (polidx >=0) idx = polidx;
+  return "";
+  //return SDPolUtil::polarizationLabel(idx, linear, stokes, linpol);
+}
+
+void Scantable::unsetSelection()
+{
+  table_ = originalTable_;
+  selector_.reset();
+}
+
+void Scantable::setSelection( const STSelector& selection )
+{
+  Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
+  if ( tab.nrow() == 0 ) {
+    throw(AipsError("Selection contains no data. Not applying it."));
+  }
+  table_ = tab;
+  selector_ = selection;
+}
+
+std::string asap::Scantable::summary( bool verbose )
+{
+  // Format header info
+  ostringstream oss;
+  oss << endl;
+  oss << asap::SEPERATOR << endl;
+  oss << " Scan Table Summary" << endl;
+  oss << asap::SEPERATOR << endl;
+  oss.flags(std::ios_base::left);
+  oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
+      << setw(15) << "IFs:" << setw(4) << nif() << endl
+      << setw(15) << "Polarisations:" << setw(4) << npol() << endl
+      << setw(15) << "Channels:"  << setw(4) << nchan() << endl;
+  oss << endl;
+  String tmp;
+  //table_.keywordSet().get("Observer", tmp);
+  oss << setw(15) << "Observer:" << table_.keywordSet().asString("Observer") << endl;
+  oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
+  table_.keywordSet().get("Project", tmp);
+  oss << setw(15) << "Project:" << tmp << endl;
+  table_.keywordSet().get("Obstype", tmp);
+  oss << setw(15) << "Obs. Type:" << tmp << endl;
+  table_.keywordSet().get("AntennaName", tmp);
+  oss << setw(15) << "Antenna Name:" << tmp << endl;
+  table_.keywordSet().get("FluxUnit", tmp);
+  oss << setw(15) << "Flux Unit:" << tmp << endl;
+  Vector<Float> vec;
+  oss << setw(15) << "Rest Freqs:";
+  if (vec.nelements() > 0) {
+      oss << setprecision(10) << vec << " [Hz]" << endl;
+  } else {
+      oss << "none" << endl;
+  }
+  oss << setw(15) << "Abcissa:" << "channel" << endl;
+  oss << selector_.print() << endl;
+  oss << endl;
+  // main table
+  String dirtype = "Position ("
+                  + MDirection::showType(dirCol_.getMeasRef().getType())
+                  + ")";
+  oss << setw(5) << "Scan"
+      << setw(15) << "Source"
+//      << setw(24) << dirtype
+      << setw(10) << "Time"
+      << setw(18) << "Integration" << endl
+      << setw(5) << "" << setw(10) << "Beam" << dirtype << endl
+      << setw(15) << "" << setw(5) << "IF"
+      << setw(8) << "Frame" << setw(16)
+      << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
+  oss << asap::SEPERATOR << endl;
+  TableIterator iter(table_, "SCANNO");
+  while (!iter.pastEnd()) {
+    Table subt = iter.table();
+    ROTableRow row(subt);
+    MEpoch::ROScalarColumn timeCol(subt,"TIME");
+    const TableRecord& rec = row.get(0);
+    oss << setw(4) << std::right << rec.asuInt("SCANNO")
+        << std::left << setw(1) << ""
+        << setw(15) << rec.asString("SRCNAME")
+        << setw(10) << formatTime(timeCol(0), false);
+    // count the cycles in the scan
+    TableIterator cyciter(subt, "CYCLENO");
+    int nint = 0;
+    while (!cyciter.pastEnd()) {
+      ++nint;
+      ++cyciter;
+    }
+    oss << setw(3) << std::right << nint  << setw(3) << " x " << std::left
+        << setw(6) <<  formatSec(rec.asFloat("INTERVAL")) << endl;
+
+    TableIterator biter(subt, "BEAMNO");
+    while (!biter.pastEnd()) {
+      Table bsubt = biter.table();
+      ROTableRow brow(bsubt);
+      MDirection::ROScalarColumn bdirCol(bsubt,"DIRECTION");
+      const TableRecord& brec = brow.get(0);
+      oss << setw(6) << "" <<  setw(10) << brec.asuInt("BEAMNO");
+      oss  << setw(24) << formatDirection(bdirCol(0)) << endl;
+      TableIterator iiter(bsubt, "IFNO");
+      while (!iiter.pastEnd()) {
+        Table isubt = iiter.table();
+        ROTableRow irow(isubt);
+        const TableRecord& irec = irow.get(0);
+        oss << std::right <<setw(8) << "" << std::left << irec.asuInt("IFNO");
+        oss << frequencies().print(irec.asuInt("FREQ_ID"));
+
+        ++iiter;
+      }
+      ++biter;
+    }
+    ++iter;
+  }
+  /// @todo implement verbose mode
+  return String(oss);
+}
+
+std::string Scantable::getTime(int whichrow, bool showdate) const
+{
+  MEpoch::ROScalarColumn timeCol(table_, "TIME");
+  MEpoch me;
+  if (whichrow > -1) {
+    me = timeCol(uInt(whichrow));
+  } else {
+    Double tm;
+    table_.keywordSet().get("UTC",tm);
+    me = MEpoch(MVEpoch(tm));
+  }
+  return formatTime(me, showdate);
+}
+
+}//namespace asap