Changeset 1819 for trunk/src

Timestamp:

08/02/10 17:28:20 (14 years ago)

Author:

Kana Sugimoto

Message:

New Development: No

JIRA Issue: No (merge alma branch to trunk)

Ready for Test: Yes

Interface Changes: No

Test Programs: regressions may work

Module(s): all single dish modules

Description:

Merged all changes in alma (r1386:1818) and newfiller (r1774:1818) branch.

Location:

trunk/src

Files:

• . (modified) (1 prop)
• AsapLogSink.cpp (copied) (copied from branches/alma/src/AsapLogSink.cpp)
• AsapLogSink.h (copied) (copied from branches/alma/src/AsapLogSink.h)
• FillerBase.cpp (copied) (copied from branches/alma/src/FillerBase.cpp)
• FillerBase.h (copied) (copied from branches/alma/src/FillerBase.h)
• FillerWrapper.h (copied) (copied from branches/alma/src/FillerWrapper.h)
• IndexedCompare.h (copied) (copied from branches/alma/src/IndexedCompare.h)
• Lorentzian1D.h (copied) (copied from branches/alma/src/Lorentzian1D.h)
• Lorentzian1D.tcc (copied) (copied from branches/alma/src/Lorentzian1D.tcc)
• Lorentzian1D2.tcc (copied) (copied from branches/alma/src/Lorentzian1D2.tcc)
• Lorentzian1DParam.h (copied) (copied from branches/alma/src/Lorentzian1DParam.h)
• Lorentzian1DParam.tcc (copied) (copied from branches/alma/src/Lorentzian1DParam.tcc)
• Makefile (copied) (copied from branches/alma/src/Makefile)
• MathUtils.cpp (modified) (2 diffs)
• MathUtils.h (modified) (2 diffs)
• NROFiller.cpp (copied) (copied from branches/alma/src/NROFiller.cpp)
• NROFiller.h (copied) (copied from branches/alma/src/NROFiller.h)
• PKSFiller.cpp (copied) (copied from branches/alma/src/PKSFiller.cpp)
• PKSFiller.h (copied) (copied from branches/alma/src/PKSFiller.h)
• RowAccumulator.cpp (modified) (1 diff)
• RowAccumulator.h (modified) (1 diff)
• SConscript (modified) (1 prop)
• STAsciiWriter.cpp (modified) (3 diffs)
• STAtmosphere.cpp (copied) (copied from branches/alma/src/STAtmosphere.cpp)
• STAtmosphere.h (copied) (copied from branches/alma/src/STAtmosphere.h)
• STCoordinate.h (copied) (copied from branches/alma/src/STCoordinate.h)
• STFITSImageWriter.cpp (copied) (copied from branches/alma/src/STFITSImageWriter.cpp)
• STFITSImageWriter.h (copied) (copied from branches/alma/src/STFITSImageWriter.h)
• STFiller.cpp (modified) (19 diffs)
• STFiller.h (modified) (5 diffs)
• STFitter.cpp (modified) (5 diffs)
• STFrequencies.cpp (modified) (3 diffs)
• STFrequencies.h (modified) (3 diffs)
• STHeader.cpp (modified) (3 diffs)
• STMath.cpp (modified) (23 diffs)
• STMath.h (modified) (9 diffs)
• STMathWrapper.h (modified) (4 diffs)
• STMolecules.cpp (modified) (7 diffs)
• STMolecules.h (modified) (3 diffs)
• STSelector.cpp (modified) (6 diffs)
• STSelector.h (modified) (3 diffs)
• STWriter.cpp (modified) (9 diffs)
• STWriter.h (modified) (2 diffs)
• Scantable.cpp (modified) (24 diffs)
• Scantable.h (modified) (9 diffs)
• ScantableWrapper.h (modified) (7 diffs)
• Templates.cpp (copied) (copied from branches/alma/src/Templates.cpp)
• python_Filler.cpp (copied) (copied from branches/alma/src/python_Filler.cpp)
• python_LogSink.cpp (copied) (copied from branches/alma/src/python_LogSink.cpp)
• python_STAtmosphere.cpp (copied) (copied from branches/alma/src/python_STAtmosphere.cpp)
• python_STCoordinate.cpp (copied) (copied from branches/alma/src/python_STCoordinate.cpp)
• python_STMath.cpp (modified) (3 diffs)
• python_STSelector.cpp (modified) (2 diffs)
• python_Scantable.cpp (modified) (4 diffs)
• python_SrcType.cpp (copied) (copied from branches/alma/src/python_SrcType.cpp)
• python_asap.cpp (modified) (2 diffs)
• python_asap.h (modified) (2 diffs)

trunk/src/MathUtils.cpp

@@ -38 +38 @@
 #include <casa/BasicSL/String.h>
 #include <scimath/Mathematics/MedianSlider.h>
+#include <casa/Exceptions/Error.h>
 
 #include <scimath/Fitting/LinearFit.h>
@@ -53 +54 @@
    String str(which);
    str.upcase();
-   if (str.contains(String("MIN"))) {
+   if (str.matches(String("MIN"))) {
       return min(data);
-   } else if (str.contains(String("MAX"))) {
+   } else if (str.matches(String("MAX"))) {
       return max(data);
-   } else if (str.contains(String("SUMSQ"))) {
+   } else if (str.matches(String("SUMSQ"))) {
       return sumsquares(data);
-   } else if (str.contains(String("SUM"))) {
+   } else if (str.matches(String("SUM"))) {
       return sum(data);
-   } else if (str.contains(String("MEAN"))) {
+   } else if (str.matches(String("MEAN"))) {
       return mean(data);
-   } else if (str.contains(String("VAR"))) {
+   } else if (str.matches(String("VAR"))) {
       return variance(data);
-   } else if (str.contains(String("STDDEV"))) {
+      } else if (str.matches(String("STDDEV"))) {
       return stddev(data);
-   } else if (str.contains(String("AVDEV"))) {
+   } else if (str.matches(String("AVDEV"))) {
       return avdev(data);
-   } else if (str.contains(String("RMS"))) {
+   } else if (str.matches(String("RMS"))) {
       uInt n = data.nelementsValid();
       return sqrt(sumsquares(data)/n);
-   } else if (str.contains(String("MED"))) {
+   } else if (str.matches(String("MEDIAN"))) {
       return median(data);
-   }
+   } else {
+      String msg = str + " is not a valid type of statistics";
+      throw(AipsError(msg));
+   } 
    return 0.0;
 }
 
+IPosition mathutil::minMaxPos(const String& which,
+                           const MaskedArray<Float>& data)
+{
+   Float minVal, maxVal;
+   IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0); 
+   minMax(minVal, maxVal, minPos, maxPos, data);
+   String str(which);
+   str.upcase();
+   if (str.contains(String("MIN"))) {
+     return minPos;
+   } else if (str.contains(String("MAX"))) {
+     return maxPos;
+   } else {
+      String msg = str + " is not a valid type of statistics";
+      throw(AipsError(msg));
+   } 
+   //return 0.0;
+}
 
 void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)

trunk/src/MathUtils.h

@@ -37 +37 @@
 #include <casa/Arrays/Vector.h>
 #include <casa/BasicSL/String.h>
+#include <casa/Arrays/IPosition.h>
 
 namespace mathutil {
@@ -79 +80 @@
                    float hwidth, int order);
 
-
 // Generate specified statistic
 float statistics(const casa::String& which,
+                 const casa::MaskedArray<casa::Float>& data);
+
+// Return a position of min or max value
+ casa::IPosition minMaxPos(const casa::String& which,
                  const casa::MaskedArray<casa::Float>& data);
 

trunk/src/RowAccumulator.cpp

@@ -152 +152 @@
   userMask_ = m;
 }
+
+// Added by TT  check the state of RowAccumulator
+casa::Bool RowAccumulator::state() const
+{
+  return initialized_;
+}
+

trunk/src/RowAccumulator.h

@@ -85 +85 @@
     */
   void reset();
+  /**
+    * check the initialization state 
+    */ 
+  casa::Bool state() const;
 
 private:

trunk/src/STAsciiWriter.cpp

@@ +1 @@
+
 //#---------------------------------------------------------------------------
 //# STAsciiWriter.cc: ASAP class to write out single dish spectra as FITS images
@@ -88 +89 @@
 
    String rootName(fileName);
-   
+
   Block<String> cols(4);
   cols[0] = String("SCANNO");
@@ -132 +133 @@
     String wcs = stable.frequencies().print(rec.asuInt("FREQ_ID"), True);
     addLine(of, "WCS", wcs);
-    addLine(of, "Rest Freq.", 
-            stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
+    std::vector<double> restfreqs= stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID"));
+    int nf = restfreqs.size();
+    //addLine(of, "Rest Freq.", 
+    //        stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
+    addLine(of, "Rest Freq.", restfreqs[0]);
+    for ( unsigned int i=1; i<nf; ++i) {
+      addLine(of, " ", restfreqs[i]);
+    }
+    ostringstream osflagrow;
+    for ( unsigned int i=0; i<t.nrow(); ++i) {
+      osflagrow << "Pol" << i << ":" << ((row.get(i).asuInt("FLAGROW") > 0) ? "True" : "False") << " ";
+    }
+    addLine(of, "Row_Flagged", String(osflagrow));
     of << setfill('#') << setw(70) << "" << setfill(' ') << endl;
 

trunk/src/STFiller.cpp

@@ -25 +25 @@
 #include <tables/Tables/TableRow.h>
 
+#include <measures/Measures/MDirection.h>
+#include <measures/Measures/MeasConvert.h>
+
 #include <atnf/PKSIO/PKSrecord.h>
 #include <atnf/PKSIO/PKSreader.h>
@@ -30 +33 @@
  #include <casa/System/ProgressMeter.h>
 #endif
+#include <casa/System/ProgressMeter.h>
+#include <atnf/PKSIO/NROReader.h>
+#include <casa/Logging/LogIO.h>
+
+#include <time.h>
+
 
 #include "STDefs.h"
@@ -45 +54 @@
   header_(0),
   table_(0),
-  refRx_(".*(e|w|_R)$")
+  refRx_(".*(e|w|_R)$"),
+  nreader_(0)
 {
 }
@@ -53 +63 @@
   header_(0),
   table_(stbl),
-  refRx_(".*(e|w|_R)$")
+  refRx_(".*(e|w|_R)$"),
+  nreader_(0)
 {
 }
@@ -61 +72 @@
   header_(0),
   table_(0),
-  refRx_(".*(e|w|_R)$")
-{
-  open(filename, whichIF, whichBeam);
+  refRx_(".*(e|w|_R)$"),
+  nreader_(0)
+{
+  open(filename, "", whichIF, whichBeam);
 }
 
@@ -71 +83 @@
 }
 
-void STFiller::open( const std::string& filename, int whichIF, int whichBeam )
+void STFiller::open( const std::string& filename, const std::string& antenna, int whichIF, int whichBeam, casa::Bool getPt )
 {
   if (table_.null())  {
@@ -93 +105 @@
   Vector<Bool> beams, ifs;
   Vector<uInt> nchans,npols;
-  if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
+
+  //
+  // if isNRO_ is true, try NROReader
+  //
+  // 2008/11/11 Takeshi Nakazato
+  isNRO_ = fileCheck() ;
+  if ( isNRO_ ) {
+    if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
+      throw(AipsError("Creation of NROReader failed")) ;
+    }
+    else {
+      openNRO( whichIF, whichBeam ) ;
+      return ;
+    }
+  }
+  //
+
+  if ( (reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
                               nchans, npols, haveXPol_,haveBase, haveSpectra
                               )) == 0 )  {
@@ -118 +147 @@
   header_->npol = max(npols);
   header_->nbeam = nBeam_;
-  
+
   Int status = reader_->getHeader(header_->observer, header_->project,
                                   header_->antennaname, header_->antennaposition,
@@ -197 +226 @@
   Vector<Int> start(nIF_, 1);
   Vector<Int> end(nIF_, 0);
-  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0]);
+  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
   table_->setHeader(*header_);
   //For MS, add the location of POINTING of the input MS so one get
   //pointing data from there, if necessary.
-  //Also find nrow in MS 
+  //Also find nrow in MS
   nInDataRow = 0;
   if (format == "MS2") {
-    Path datapath(inName); 
+    Path datapath(inName);
     String ptTabPath = datapath.absoluteName();
     Table inMS(ptTabPath);
@@ -216 +245 @@
     }
   }
+  String freqFrame = header_->freqref;
+  //translate frequency reference frame back to
+  //MS style (as PKSMS2reader converts the original frame
+  //in FITS standard style)
+  if (freqFrame == "TOPOCENT") {
+    freqFrame = "TOPO";
+  } else if (freqFrame == "GEOCENER") {
+    freqFrame = "GEO";
+  } else if (freqFrame == "BARYCENT") {
+    freqFrame = "BARY";
+  } else if (freqFrame == "GALACTOC") {
+    freqFrame = "GALACTO";
+  } else if (freqFrame == "LOCALGRP") {
+    freqFrame = "LGROUP";
+  } else if (freqFrame == "CMBDIPOL") {
+    freqFrame = "CMB";
+  } else if (freqFrame == "SOURCE") {
+    freqFrame = "REST";
+  }
+  // set both "FRAME" and "BASEFRAME"
+  table_->frequencies().setFrame(freqFrame, false);
+  table_->frequencies().setFrame(freqFrame,true);
   //table_->focus().setParallactify(true);
 }
@@ -222 +273 @@
 {
   delete reader_;reader_=0;
+  delete nreader_;nreader_=0;
   delete header_;header_=0;
   table_ = 0;
@@ -229 +281 @@
 {
   int status = 0;
+
+  //
+  // for NRO data
+  //
+  // 2008/11/12 Takeshi Nakazato
+  if ( isNRO_ ) {
+    status = readNRO() ;
+    return status ;
+  }
+  //
+
+/**
+  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
+  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
+    humidity, parAngle, pressure, temperature, windAz, windSpeed;
+  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
+  String          fieldName, srcName, tcalTime, obsType;
+  Vector<Float>   calFctr, sigma, tcal, tsys;
+  Matrix<Float>   baseLin, baseSub;
+  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
+  Matrix<Float>   spectra;
+  Matrix<uChar>   flagtra;
+  Complex         xCalFctr;
+  Vector<Complex> xPol;
+**/
 
   Double min = 0.0;
@@ -236 +313 @@
 #endif
   PKSrecord pksrec;
+  pksrec.srcType=-1;
   int n = 0;
+  bool isGBTFITS = false ;
+  if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
+    FILE *fp = fopen( filename_.c_str(), "r" ) ;
+    fseek( fp, 640, SEEK_SET ) ;
+    char buf[81] ;
+    fread( buf, 80, 1, fp ) ;
+    buf[80] = '\0' ;
+    if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
+      isGBTFITS = true ;
+    }
+    fclose( fp ) ;
+  }
   while ( status == 0 ) {
     status = reader_->read(pksrec);
@@ -288 +378 @@
     //*srcnCol = pksrec.srcName;//.before(rx2);
     *srctCol = match;
+    if ( pksrec.srcType != -1 ) {
+      *srctCol = pksrec.srcType ;
+    }
     RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
     *beamCol = pksrec.beamNo-beamOffset_-1;
@@ -296 +389 @@
     RecordFieldPtr<uInt> ifCol(rec, "IFNO");
     *ifCol = pksrec.IFno-ifOffset_- 1;
-    uInt id;
-    /// @todo this has to change when nchan isn't global anymore
-    id = table_->frequencies().addEntry(Double(header_->nchan/2),
-                                        pksrec.refFreq, pksrec.freqInc);
+    uInt id = table_->frequencies().addEntry(Double(pksrec.spectra.nrow()/2),
+                                             pksrec.refFreq, pksrec.freqInc);
     RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
     *mfreqidCol = id;
+    //*ifCol = id;
 
     id = table_->molecules().addEntry(pksrec.restFreq);
@@ -317 +409 @@
 
     RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
-    id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi, 
+    id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
                                   pksrec.focusTan, pksrec.focusRot);
     *mfocusidCol = id;
@@ -335 +427 @@
     // into 2-4 rows in the scantable
     Vector<Float> tsysvec(1);
-    // Why is pksrec.spectra.ncolumn() == 3 for haveXPol_ == True
+    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
     uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
     for ( uInt i=0; i< npol; ++i ) {
       tsysvec = pksrec.tsys(i);
       *tsysCol = tsysvec;
-      *polnoCol = i;
+      if (isGBTFITS)
+        *polnoCol = pksrec.polNo ;
+      else
+        *polnoCol = i;
 
       *specCol = pksrec.spectra.column(i);
@@ -347 +442 @@
       row.put(table_->table().nrow()-1, rec);
     }
+
+    RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
+    *flagrowCol = pksrec.flagrow;
+
     if ( haveXPol_[0] ) {
       // no tsys given for xpol, so emulate it
@@ -381 +480 @@
 }
 
+/**
+ * For NRO data
+ *
+ * 2008/11/11 Takeshi Nakazato
+ **/
+void STFiller::openNRO( int whichIF, int whichBeam )
+{
+  // open file
+  // DEBUG
+  time_t t0 ;
+  time( &t0 ) ;
+  tm *ttm = localtime( &t0 ) ;
+  LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
+//   cout << "STFiller::openNRO()  Start time = " << t0
+//        << " ("
+//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+//        << " "
+//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+//        << ")" << endl ;
+  os << "Start time = " << t0
+     << " ("
+     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+     << " "
+     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+     << ")" << LogIO::POST ;
+
+  // fill STHeader
+  header_ = new STHeader() ;
+
+  if ( nreader_->getHeaderInfo( header_->nchan,
+                                header_->npol,
+                                nIF_,
+                                nBeam_,
+                                header_->observer,
+                                header_->project,
+                                header_->obstype,
+                                header_->antennaname,
+                                header_->antennaposition,
+                                header_->equinox,
+                                header_->freqref,
+                                header_->reffreq,
+                                header_->bandwidth,
+                                header_->utc,
+                                header_->fluxunit,
+                                header_->epoch,
+                                header_->poltype ) ) {
+//     cout << "STFiller::openNRO()  Failed to get header information." << endl ;
+//     return ;
+    throw( AipsError("Failed to get header information.") ) ;
+  }
+
+  // set FRAME and BASEFRAME keyword of FREQUENCIES table
+  if ( header_->freqref != "TOPO" ) {
+    table_->frequencies().setFrame( header_->freqref, false ) ;
+    table_->frequencies().setFrame( header_->freqref, true ) ;
+  }
+
+  ifOffset_ = 0;
+  vector<Bool> ifs = nreader_->getIFs() ;
+  if ( whichIF >= 0 ) {
+    if ( whichIF >= 0 && whichIF < nIF_ ) {
+      for ( int i = 0 ; i < nIF_ ; i++ )
+        ifs[i] = False ;
+      ifs[whichIF] = True ;
+      header_->nif = 1;
+      nIF_ = 1;
+      ifOffset_ = whichIF;
+    } else {
+      delete reader_;
+      reader_ = 0;
+      delete header_;
+      header_ = 0;
+      throw(AipsError("Illegal IF selection"));
+    }
+  }
+
+  beamOffset_ = 0;
+  vector<Bool> beams = nreader_->getBeams() ;
+  if (whichBeam>=0) {
+    if (whichBeam>=0 && whichBeam<nBeam_) {
+      for ( int i = 0 ; i < nBeam_ ; i++ )
+        beams[i] = False ;
+      beams[whichBeam] = True;
+      header_->nbeam = 1;
+      nBeam_ = 1;
+      beamOffset_ = whichBeam;
+    } else {
+      delete reader_;
+      reader_ = 0;
+      delete header_;
+      header_ = 0;
+      throw(AipsError("Illegal Beam selection"));
+    }
+  }
+
+  header_->nbeam = nBeam_ ;
+  header_->nif = nIF_ ;
+
+  // set header
+  table_->setHeader( *header_ ) ;
+
+  // DEBUG
+  time_t t1 ;
+  time( &t1 ) ;
+  ttm = localtime( &t1 ) ;
+//   cout << "STFiller::openNRO()  End time = " << t1
+//        << " ("
+//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+//        << " "
+//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+//        << ")" << endl ;
+//   cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
+  os << "End time = " << t1
+     << " ("
+     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+     << " "
+     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+     << ")" << endl ;
+  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
+  os.post() ;
+  //
+
+  return ;
+}
+
+int STFiller::readNRO()
+{
+  // DEBUG
+  time_t t0 ;
+  time( &t0 ) ;
+  tm *ttm = localtime( &t0 ) ;
+  LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
+//   cout << "STFiller::readNRO()  Start time = " << t0
+//        << " ("
+//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+//        << " "
+//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+//        << ")" << endl ;
+  os << "Start time = " << t0
+     << " ("
+     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+     << " "
+     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+     << ")" << LogIO::POST ;
+  //
+
+  // fill row
+  uInt id ;
+  uInt imax = nreader_->getRowNum() ;
+  vector< vector<double > > freqs ;
+  uInt i = 0 ;
+  int count = 0 ;
+  uInt scanno ;
+  uInt cycleno ;
+  uInt beamno ;
+  uInt polno ;
+  vector<double> fqs ;
+  Vector<Double> restfreq ;
+  uInt refbeamno ;
+  Double scantime ;
+  Double interval ;
+  String srcname ;
+  String fieldname ;
+  Array<Float> spectra ;
+  Array<uChar> flagtra ;
+  Array<Float> tsys ;
+  Array<Double> direction ;
+  Float azimuth ;
+  Float elevation ;
+  Float parangle ;
+  Float opacity ;
+  uInt tcalid ;
+  Int fitid ;
+  uInt focusid ;
+  Float temperature ;
+  Float pressure ;
+  Float humidity ;
+  Float windvel ;
+  Float winddir ;
+  Double srcvel ;
+  Array<Double> propermotion ;
+  Vector<Double> srcdir ;
+  Array<Double> scanrate ;
+  for ( i = 0 ; i < imax ; i++ ) {
+    string scanType = nreader_->getScanType( i ) ;
+    Int srcType = -1 ;
+    if ( scanType.compare( 0, 2, "ON") == 0 ) {
+      // os << "ON srcType: " << i << LogIO::POST ;
+      srcType = 0 ;
+    }
+    else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
+      //os << "OFF srcType: " << i << LogIO::POST ;
+      srcType = 1 ;
+    }
+    else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
+      //os << "ZERO srcType: " << i << LogIO::POST ;
+      srcType = 2 ;
+    }
+    else {
+      //os << "Undefined srcType: " << i << LogIO::POST ;
+      srcType = 3 ;
+    }
+
+    // if srcType is 2 (ZERO scan), ignore scan
+    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
+      TableRow row( table_->table() ) ;
+      TableRecord& rec = row.record();
+
+      if ( nreader_->getScanInfo( i,
+                                  scanno,
+                                  cycleno,
+                                  beamno,
+                                  polno,
+                                  fqs,
+                                  restfreq,
+                                  refbeamno,
+                                  scantime,
+                                  interval,
+                                  srcname,
+                                  fieldname,
+                                  spectra,
+                                  flagtra,
+                                  tsys,
+                                  direction,
+                                  azimuth,
+                                  elevation,
+                                  parangle,
+                                  opacity,
+                                  tcalid,
+                                  fitid,
+                                  focusid,
+                                  temperature,
+                                  pressure,
+                                  humidity,
+                                  windvel,
+                                  winddir,
+                                  srcvel,
+                                  propermotion,
+                                  srcdir,
+                                  scanrate ) ) {
+//         cerr << "STFiller::readNRO()  Failed to get scan information." << endl ;
+//         return 1 ;
+        throw( AipsError("Failed to get scan information.") ) ;
+      }
+
+      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
+      *scannoCol = scanno ;
+      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
+      *cyclenoCol = cycleno ;
+      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
+      *beamCol = beamno ;
+      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
+      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
+      *polnoCol = polno ;
+      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
+      if ( freqs.size() == 0 ) {
+        id = table_->frequencies().addEntry( Double( fqs[0] ),
+                                             Double( fqs[1] ),
+                                             Double( fqs[2] ) ) ;
+        *mfreqidCol = id ;
+        *ifCol = id ;
+        freqs.push_back( fqs ) ;
+      }
+      else {
+        int iadd = -1 ;
+        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
+          //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
+          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
+          //os << "fdiff = " << fdiff << LogIO::POST ;
+          if ( fdiff < 1.0e-8 ) {
+            iadd = iif ;
+            break ;
+          }
+        }
+        if ( iadd == -1 ) {
+          id = table_->frequencies().addEntry( Double( fqs[0] ),
+                                               Double( fqs[1] ),
+                                               Double( fqs[2] ) ) ;
+          *mfreqidCol = id ;
+          *ifCol = id ;
+          freqs.push_back( fqs ) ;
+        }
+        else {
+          *mfreqidCol = iadd ;
+          *ifCol = iadd ;
+        }
+      }
+      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
+      id = table_->molecules().addEntry( restfreq ) ;
+      *molidCol = id ;
+      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
+      *rbCol = refbeamno ;
+      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
+      *mjdCol = scantime ;
+      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
+      *intervalCol = interval ;
+      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
+      *srcnCol = srcname ;
+      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
+      *srctCol = srcType ;
+      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
+      *fieldnCol = fieldname ;
+      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
+      *specCol = spectra ;
+      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
+      *flagCol = flagtra ;
+      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
+      *tsysCol = tsys ;
+      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
+      *dirCol = direction ;
+      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
+      *azCol = azimuth ;
+      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
+      *elCol = elevation ;
+      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
+      *parCol = parangle ;
+      RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
+      *tauCol = opacity ;
+      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
+      *mcalidCol = tcalid ;
+      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
+      *fitCol = fitid ;
+      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
+      *mfocusidCol = focusid ;
+      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
+      id = table_->weather().addEntry( temperature,
+                                       pressure,
+                                       humidity,
+                                       windvel,
+                                       winddir ) ;
+      *mweatheridCol = id ;
+      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
+      *svelCol = srcvel ;
+      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
+      *spmCol = propermotion ;
+      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
+      *sdirCol = srcdir ;
+      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
+      *srateCol = scanrate ;
+
+      table_->table().addRow() ;
+      row.put(table_->table().nrow()-1, rec) ;
+    }
+    else {
+      count++ ;
+    }
+    // DEBUG
+    //int rownum = nreader_->getRowNum() ;
+    //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
+    //
+  }
+
+  // DEBUG
+  time_t t1 ;
+  time( &t1 ) ;
+  ttm = localtime( &t1 ) ;
+//   cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
+//   cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
+//        << count << " \"ZERO\" scans)" << endl ;
+//   cout << "STFiller::readNRO()  End time = " << t1
+//        << " ("
+//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+//        << " "
+//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+//        << ")" << endl ;
+//   cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
+  os << "Processed " << i << " rows" << endl ;
+  os << "Added " << i - count << " rows (ignored "
+     << count << " \"ZERO\" scans)" << endl ;
+  os.post() ;
+  os << "End time = " << t1
+     << " ("
+     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
+     << " "
+     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
+     << ")" << endl ;
+  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
+  os.post() ;
+  //
+
+  return 0 ;
+}
+
+Bool STFiller::fileCheck()
+{
+  bool bval = false ;
+
+  // if filename_ is directory, return false
+  File inFile( filename_ ) ;
+  if ( inFile.isDirectory() )
+    return bval ;
+
+  // if beginning of header data is "RW", return true
+  // otherwise, return false ;
+  FILE *fp = fopen( filename_.c_str(), "r" ) ;
+  char buf[9] ;
+  char buf2[80] ;
+  fread( buf, 4, 1, fp ) ;
+  buf[4] = '\0' ;
+  fseek( fp, 640, SEEK_SET ) ;
+  fread( buf2, 80, 1, fp ) ;
+  if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
+    bval = true ;
+  }
+  fclose( fp ) ;
+  return bval ;
+}
+
 }//namespace asap

trunk/src/STFiller.h

@@ -27 +27 @@
 
 class PKSreader;
+class NROReader;
 
 namespace asap {
@@ -61 +62 @@
     */
   explicit STFiller( const std::string& filename, int whichIF=-1,
-                  int whichBeam=-1 );
+                     int whichBeam=-1 );
 
   /**
@@ -75 +76 @@
    * @exception AipsError Creation of PKSreader failed
    */
-  void open( const std::string& filename, int whichIF=-1, int whichBeam=-1 );
+  void open( const std::string& filename, const std::string& antenna, int whichIF=-1, int whichBeam=-1, casa::Bool getPt=casa::False );
 
   /**
@@ -93 +94 @@
   casa::CountedPtr<Scantable> getTable() const { return table_;}
 
+  /**
+   * For NRO data
+   *
+   * 2008/11/11 Takeshi Nakazato
+   *
+   * openNRO  : NRO version of open(), which performs to open file and 
+   *            read header data.
+   *  
+   * readNRO  : NRO version of read(), which performs to read scan 
+   *            records.
+   *
+   * fileCheck: Identify a type (NRO data or not) of filename_.
+   **/
+  void openNRO( int whichIF=-1, int whichBeam=-1 ) ;
+  int readNRO() ;
+  casa::Bool fileCheck() ;
+
   void setReferenceExpr(const std::string& rx) { refRx_ = rx; }
 
@@ -105 +123 @@
   casa::Vector<casa::Bool> haveXPol_;
   casa::String refRx_;
+  NROReader *nreader_ ;
+  casa::Bool isNRO_ ;
 };
 

trunk/src/STFitter.cpp

@@ -32 +32 @@
 #include <casa/Arrays/ArrayMath.h>
 #include <casa/Arrays/ArrayLogical.h>
+#include <casa/Logging/LogIO.h>
 #include <scimath/Fitting.h>
 #include <scimath/Fitting/LinearFit.h>
@@ -37 +38 @@
 #include <scimath/Functionals/CompoundFunction.h>
 #include <scimath/Functionals/Gaussian1D.h>
+#include "Lorentzian1D.h"
 #include <scimath/Functionals/Polynomial.h>
 #include <scimath/Mathematics/AutoDiff.h>
@@ -146 +148 @@
     funcs_.resize(1);
     funcs_[0] = new Polynomial<Float>(ncomp);
+  } else if (expr == "lorentz") {
+    if (ncomp < 1) throw (AipsError("Need at least one lorentzian to fit."));
+    funcs_.resize(ncomp);
+    for (Int k=0; k<ncomp; ++k) {
+      funcs_[k] = new Lorentzian1D<Float>();
+    }
   } else {
-    cerr << " compiled functions not yet implemented" << endl;
+    //cerr << " compiled functions not yet implemented" << endl;
+    LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ;
+    os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST;
     //funcs_.resize(1);
     //funcs_[0] = new CompiledFunction<Float>();
@@ -227 +237 @@
             (funcs_[0]->parameters())[i] =  tmppar[i];
         }
+    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
+        uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                (funcs_[j]->parameters())[i] = tmppar[count];
+                parameters_[count] = tmppar[count];
+                ++count;
+            }
+        }
     }
     // reset
@@ -260 +279 @@
             funcs_[0]->mask(i) =  !fixed[i];
         }
+    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
+      uInt count = 0;
+        for (uInt j=0; j < funcs_.nelements(); ++j) {
+            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
+                funcs_[j]->mask(i) = !fixed[count];
+                fixedpar_[count] = fixed[count];
+                ++count;
+            }
+        }
     }
     return true;

trunk/src/STFrequencies.cpp

@@ -128 +128 @@
 }
 
+void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
+{
+  Table t = table_(table_.col("ID") == Int(id) );
+  if (t.nrow() == 0 ) {
+    throw(AipsError("STFrequencies::getEntry - freqID out of range"));
+  }
+  for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
+    uInt fid ;
+    idCol_.get( i, fid ) ;
+    if ( fid == id ) {
+      refpixCol_.put( i, refpix ) ;
+      refvalCol_.put( i, refval ) ;
+      incrCol_.put( i, inc ) ;
+    }
+  }
+}
+
 SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
 {
@@ -145 +162 @@
 }
 
+/**
 SpectralCoordinate
   STFrequencies::getSpectralCoordinate( const MDirection& md,
@@ -150 +168 @@
                                         const MEpoch& me,
                                         Double restfreq, uInt id ) const
+**/
+SpectralCoordinate
+  STFrequencies::getSpectralCoordinate( const MDirection& md,
+                                              const MPosition& mp,
+                                              const MEpoch& me,
+                                              Vector<Double> restfreq, uInt id ) const
 {
   SpectralCoordinate spc = getSpectralCoordinate(id);
-  spc.setRestFrequency(restfreq, True);
+  //spc.setRestFrequency(restfreq, True);
+  // for now just use the first rest frequency 
+  if (restfreq.nelements()==0 ) {
+    restfreq.resize(1);
+    restfreq[0] = 0;
+  }
+  spc.setRestFrequency(restfreq[0], True);
   if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
     throw(AipsError("Couldn't convert frequency frame."));

trunk/src/STFrequencies.h

@@ -59 +59 @@
                  casa::Double& inc, casa::uInt id );
 
+  /***
+   * Set the frequency values for a specific id via references
+   * @param refpix the reference pixel
+   * @param refval the reference value
+   * @param inc    the increment
+   * @param id     the identifier
+   *
+   * 17/09/2008 Takeshi Nakazato
+   ***/
+  void setEntry( casa::Double refpix, casa::Double refval, 
+                 casa::Double inc, casa::uInt id ) ;
+
 
   bool conformant(const STFrequencies& other) const;
@@ -69 +81 @@
   casa::SpectralCoordinate getSpectralCoordinate( casa::uInt freqID ) const;
 
+  /**
   casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
                                                   const casa::MPosition& mp,
@@ -75 +88 @@
                                                   casa::uInt freqID
                                                   ) const;
+  **/
+  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
+                                                  const casa::MPosition& mp,
+                                                  const casa::MEpoch& me,
+                                                  casa::Vector<casa::Double> restfreq,
+                                                  casa::uInt freqID
+                                                  ) const;
 
   /**

trunk/src/STHeader.cpp

@@ -37 +37 @@
 #include <casa/Arrays/IPosition.h>
 #include <casa/Quanta/MVTime.h>
+#include <casa/Logging/LogIO.h>
 
-
+#include <sstream>
 
 #include "STDefs.h"
@@ -79 +80 @@
   MVTime mvt(this->utc);
   mvt.setFormat(MVTime::YMD);
-  cout << "Observer: " << this->observer << endl
-       << "Project: " << this->project << endl
-       << "Obstype: " << this->obstype << endl
-       << "Antenna: " << this->antennaname << endl
-       << "Ant. Position: " << this->antennaposition << endl
-       << "Equinox: " << this->equinox << endl
-       << "Freq. ref.: " << this->freqref << endl
-       << "Ref. frequency: " << this->reffreq << endl
-       << "Bandwidth: "  << this->bandwidth << endl
-       << "Time (utc): "
-       << mvt
-       << endl;
+//   cout << "Observer: " << this->observer << endl
+//        << "Project: " << this->project << endl
+//        << "Obstype: " << this->obstype << endl
+//        << "Antenna: " << this->antennaname << endl
+//        << "Ant. Position: " << this->antennaposition << endl
+//        << "Equinox: " << this->equinox << endl
+//        << "Freq. ref.: " << this->freqref << endl
+//        << "Ref. frequency: " << this->reffreq << endl
+//        << "Bandwidth: "  << this->bandwidth << endl
+//        << "Time (utc): "
+//        << mvt
+//        << endl;
+  LogIO os( LogOrigin( "STHeader", "print()", WHERE ) ) ;
+  os << "Observer: " << this->observer << endl
+     << "Project: " << this->project << endl
+     << "Obstype: " << this->obstype << endl
+     << "Antenna: " << this->antennaname << endl
+     << "Ant. Position: " << this->antennaposition << endl
+     << "Equinox: " << this->equinox << endl
+     << "Freq. ref.: " << this->freqref << endl
+     << "Ref. frequency: " << this->reffreq << endl
+     << "Bandwidth: "  << this->bandwidth << endl
+     << "Time (utc): "
+     << mvt
+     << LogIO::POST ;
   //setprecision(10) << this->utc << endl;
 }
@@ -129 +143 @@
 {
    if (n_>0) {
-      cerr << "Source    ID" << endl;
-      for (uInt i=0; i<n_; i++) {
-         cout << setw(11) << source_(i) << ID_(i) << endl;
-      }
+//       cerr << "Source    ID" << endl;
+//       for (uInt i=0; i<n_; i++) {
+//          cout << setw(11) << source_(i) << ID_(i) << endl;
+     LogIO os( LogOrigin( "SDDataDesc", "summary()", WHERE ) ) ;
+     ostringstream oss ;
+     oss << "Source    ID" << endl;
+     for (uInt i=0; i<n_; i++) {
+       oss << setw(11) << source_(i) << ID_(i) << endl;
+     }
+     os << oss.str() << LogIO::POST ;
    }
 }

trunk/src/STMath.cpp

@@ -44 +44 @@
 #include <scimath/Functionals/Polynomial.h>
 
+#include <atnf/PKSIO/SrcType.h>
+
+#include <casa/Logging/LogIO.h>
+#include <sstream>
+
 #include "MathUtils.h"
 #include "RowAccumulator.h"
@@ -53 +58 @@
 
 using namespace asap;
+
+// tolerance for direction comparison (rad)
+#define TOL_OTF    1.0e-15
+#define TOL_POINT  2.9088821e-4  // 1 arcmin
 
 STMath::STMath(bool insitu) :
@@ -70 +79 @@
                  const std::string& avmode)
 {
+  LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
   if ( avmode == "SCAN" && in.size() != 1 )
     throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
                     "Use merge first."));
   WeightType wtype = stringToWeight(weight);
+
+  // check if OTF observation
+  String obstype = in[0]->getHeader().obstype ;
+  Double tol = 0.0 ;
+  if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
+    tol = TOL_OTF ;
+  }
+  else {
+    tol = TOL_POINT ;
+  }
 
   // output
@@ -113 +133 @@
   }
   if ( avmode == "SCAN"  && in.size() == 1) {
-    cols.resize(4);
-    cols[3] = String("SCANNO");
+    //cols.resize(4);
+    //cols[3] = String("SCANNO");
+    cols.resize(5);
+    cols[3] = String("SRCNAME");
+    cols[4] = String("SCANNO");
   }
   uInt outrowCount = 0;
   TableIterator iter(baset, cols);
+//   int count = 0 ;
   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);
-    // re-index to 0
-    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
-      scanColOut.put(outrowCount, uInt(0));
-    }
-    ++outrowCount;
+//     // copy the first row of this selection into the new table
+//     tout.addRow();
+//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
+//     // re-index to 0
+//     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
+//       scanColOut.put(outrowCount, uInt(0));
+//     }
+//     ++outrowCount;
+    MDirection::ScalarColumn dircol ;
+    dircol.attach( subt, "DIRECTION" ) ;
+    Int length = subt.nrow() ;
+    vector< Vector<Double> > dirs ;
+    vector<int> indexes ;
+    for ( Int i = 0 ; i < length ; i++ ) {
+      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
+      //os << << count++ << ": " ;
+      //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
+      bool adddir = true ;
+      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
+        //if ( allTrue( t == dirs[j] ) ) {
+        Double dx = t[0] - dirs[j][0] ;
+        Double dy = t[1] - dirs[j][1] ;
+        Double dd = sqrt( dx * dx + dy * dy ) ;
+        //if ( allNearAbs( t, dirs[j], tol ) ) {
+        if ( dd <= tol ) {
+          adddir = false ;
+          break ;
+        }
+      }
+      if ( adddir ) {
+        dirs.push_back( t ) ;
+        indexes.push_back( i ) ;
+      }
+    }
+    uInt rowNum = dirs.size() ;
+    tout.addRow( rowNum ) ;
+    for ( uInt i = 0 ; i < rowNum ; i++ ) {
+      TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
+      // re-index to 0
+      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
+        scanColOut.put(outrowCount+i, uInt(0));
+      }        
+    }
+    outrowCount += rowNum ;
     ++iter;
   }
-
   RowAccumulator acc(wtype);
   Vector<Bool> cmask(mask);
@@ -155 +214 @@
         subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
       } else if (avmode == "SCAN") {
-        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
+        //subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
+        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO"))
+                         && basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
       } else {
         subt = basesubt;
       }
+
+      vector<uInt> removeRows ;
+      uInt nrsubt = subt.nrow() ;
+      for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
+        //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
+        Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
+        Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
+        double dx = x0[0] - x1[0] ;
+        double dy = x0[0] - x1[0] ;
+        Double dd = sqrt( dx * dx + dy * dy ) ;
+        //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
+        if ( dd > tol ) {
+          removeRows.push_back( irow ) ;
+        }
+      }
+      if ( removeRows.size() != 0 ) {
+        subt.removeRow( removeRows ) ;
+      }
+      
+      if ( nrsubt == removeRows.size() )
+        throw(AipsError("Averaging data is empty.")) ;
+
       specCol.attach(subt,"SPECTRA");
       flagCol.attach(subt,"FLAGTRA");
@@ -192 +275 @@
     }
     //write out
-    Vector<uChar> flg(msk.shape());
-    convertArray(flg, !msk);
-    flagColOut.put(i, flg);
-    specColOut.put(i, acc.getSpectrum());
-    tsysColOut.put(i, acc.getTsys());
-    intColOut.put(i, acc.getInterval());
-    mjdColOut.put(i, acc.getTime());
-    // we should only have one cycle now -> reset it to be 0
-    // frequency switched data has different CYCLENO for different IFNO
-    // which requires resetting this value
-    cycColOut.put(i, uInt(0));
+    if (acc.state()) {
+      Vector<uChar> flg(msk.shape());
+      convertArray(flg, !msk);
+      flagColOut.put(i, flg);
+      specColOut.put(i, acc.getSpectrum());
+      tsysColOut.put(i, acc.getTsys());
+      intColOut.put(i, acc.getInterval());
+      mjdColOut.put(i, acc.getTime());
+      // we should only have one cycle now -> reset it to be 0
+      // frequency switched data has different CYCLENO for different IFNO
+      // which requires resetting this value
+      cycColOut.put(i, uInt(0));
+    } else {
+      ostringstream oss;
+      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
+      pushLog(String(oss));
+    }
     acc.reset();
   }
   if (rowstodelete.nelements() > 0) {
+    //cout << rowstodelete << endl;
+    os << rowstodelete << LogIO::POST ;
     tout.removeRow(rowstodelete);
     if (tout.nrow() == 0) {
@@ -219 +310 @@
                           const std::string& avmode )
 {
+  // check if OTF observation
+  String obstype = in->getHeader().obstype ;
+  Double tol = 0.0 ;
+  if ( obstype.find( "OTF" ) != String::npos ) {
+    tol = TOL_OTF ;
+  }
+  else {
+    tol = TOL_POINT ;
+  }
+
   // clone as this is non insitu
   bool insitu = insitu_;
@@ -250 +351 @@
     flagCol.attach(subt,"FLAGTRA");
     tsysCol.attach(subt,"TSYS");
-    tout.addRow();
-    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
-    if ( avmode != "SCAN") {
-      scanColOut.put(outrowCount, uInt(0));
-    }
-    Vector<Float> tmp;
-    specCol.get(0, tmp);
-    uInt nchan = tmp.nelements();
-    // have to do channel by channel here as MaskedArrMath
-    // doesn't have partialMedians
-    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
-    Vector<Float> outspec(nchan);
-    Vector<uChar> outflag(nchan,0);
-    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
-    for (uInt i=0; i<nchan; ++i) {
-      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
-      MaskedArray<Float> ma = maskedArray(specs,flags);
-      outspec[i] = median(ma);
-      if ( allEQ(ma.getMask(), False) )
-        outflag[i] = userflag;// flag data
-    }
-    outtsys[0] = median(tsysCol.getColumn());
-    specColOut.put(outrowCount, outspec);
-    flagColOut.put(outrowCount, outflag);
-    tsysColOut.put(outrowCount, outtsys);
-    Double intsum = sum(intCol.getColumn());
-    intColOut.put(outrowCount, intsum);
-    ++outrowCount;
+//     tout.addRow();
+//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
+//     if ( avmode != "SCAN") {
+//       scanColOut.put(outrowCount, uInt(0));
+//     }
+//     Vector<Float> tmp;
+//     specCol.get(0, tmp);
+//     uInt nchan = tmp.nelements();
+//     // have to do channel by channel here as MaskedArrMath
+//     // doesn't have partialMedians
+//     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
+//     Vector<Float> outspec(nchan);
+//     Vector<uChar> outflag(nchan,0);
+//     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
+//     for (uInt i=0; i<nchan; ++i) {
+//       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
+//       MaskedArray<Float> ma = maskedArray(specs,flags);
+//       outspec[i] = median(ma);
+//       if ( allEQ(ma.getMask(), False) )
+//         outflag[i] = userflag;// flag data
+//     }
+//     outtsys[0] = median(tsysCol.getColumn());
+//     specColOut.put(outrowCount, outspec);
+//     flagColOut.put(outrowCount, outflag);
+//     tsysColOut.put(outrowCount, outtsys);
+//     Double intsum = sum(intCol.getColumn());
+//     intColOut.put(outrowCount, intsum);
+//     ++outrowCount;
+//     ++iter;
+    MDirection::ScalarColumn dircol ;
+    dircol.attach( subt, "DIRECTION" ) ;
+    Int length = subt.nrow() ;
+    vector< Vector<Double> > dirs ;
+    vector<int> indexes ;
+    for ( Int i = 0 ; i < length ; i++ ) {
+      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
+      bool adddir = true ;
+      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
+        //if ( allTrue( t == dirs[j] ) ) {
+        Double dx = t[0] - dirs[j][0] ;
+        Double dy = t[1] - dirs[j][1] ;
+        Double dd = sqrt( dx * dx + dy * dy ) ;
+        //if ( allNearAbs( t, dirs[j], tol ) ) {
+        if ( dd <= tol ) {
+          adddir = false ;
+          break ;
+        }
+      }
+      if ( adddir ) {
+        dirs.push_back( t ) ;
+        indexes.push_back( i ) ;
+      }
+    }
+    uInt rowNum = dirs.size() ;
+    tout.addRow( rowNum );
+    for ( uInt i = 0 ; i < rowNum ; i++ ) {
+      TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
+      if ( avmode != "SCAN") {
+        //scanColOut.put(outrowCount+i, uInt(0));
+      }
+    }
+    MDirection::ScalarColumn dircolOut ;
+    dircolOut.attach( tout, "DIRECTION" ) ;
+    for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
+      Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
+      Vector<Float> tmp;
+      specCol.get(0, tmp);
+      uInt nchan = tmp.nelements();
+      // have to do channel by channel here as MaskedArrMath
+      // doesn't have partialMedians
+      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
+      // mask spectra for different DIRECTION
+      for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
+        Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
+        //if ( t[0] != direction[0] || t[1] != direction[1] ) {
+        Double dx = t[0] - direction[0] ;
+        Double dy = t[1] - direction[1] ;
+        Double dd = sqrt( dx * dx + dy * dy ) ;
+        //if ( !allNearAbs( t, direction, tol ) ) {
+        if ( dd > tol ) {
+          flags[jrow] = userflag ;
+        }
+      }
+      Vector<Float> outspec(nchan);
+      Vector<uChar> outflag(nchan,0);
+      Vector<Float> outtsys(1);/// @fixme when tsys is channel based
+      for (uInt i=0; i<nchan; ++i) {
+        Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
+        MaskedArray<Float> ma = maskedArray(specs,flags);
+        outspec[i] = median(ma);
+        if ( allEQ(ma.getMask(), False) )
+          outflag[i] = userflag;// flag data
+      }
+      outtsys[0] = median(tsysCol.getColumn());
+      specColOut.put(outrowCount+irow, outspec);
+      flagColOut.put(outrowCount+irow, outflag);
+      tsysColOut.put(outrowCount+irow, outtsys);
+      Vector<Double> integ = intCol.getColumn() ;
+      MaskedArray<Double> mi = maskedArray( integ, flags ) ;
+      Double intsum = sum(mi);
+      intColOut.put(outrowCount+irow, intsum);
+    }
+    outrowCount += rowNum ;
     ++iter;
   }
@@ -323 +499 @@
       if ( tsys ) {
         ts += val;
+        tsysCol.put(i, ts);
+      }
+    }
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
+                                              const std::vector<float> val,
+                                              const std::string& mode,
+                                              const std::string& opmode, 
+                                              bool tsys )
+{
+  CountedPtr< Scantable > out ;
+  if ( opmode == "channel" ) {
+    out = arrayOperateChannel( in, val, mode, tsys ) ;
+  }
+  else if ( opmode == "row" ) {
+    out = arrayOperateRow( in, val, mode, tsys ) ;
+  }
+  else {
+    throw( AipsError( "Unknown array operation mode." ) ) ;
+  }
+  return out ;
+}
+
+CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
+                                                     const std::vector<float> val,
+                                                     const std::string& mode,
+                                                     bool tsys )
+{
+  if ( val.size() == 1 ){
+    return unaryOperate( in, val[0], mode, tsys ) ;
+  }
+
+  // conformity of SPECTRA and TSYS
+  if ( tsys ) {
+    TableIterator titer(in->table(), "IFNO");
+    while ( !titer.pastEnd() ) {
+      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
+      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
+      Array<Float> spec = specCol.getColumn() ;
+      Array<Float> ts = tsysCol.getColumn() ;
+      if ( !spec.conform( ts ) ) {
+        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
+      }
+      titer.next() ;
+    }
+  }
+
+  // check if all spectra in the scantable have the same number of channel
+  vector<uInt> nchans;
+  vector<uInt> ifnos = in->getIFNos() ;
+  for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
+    nchans.push_back( in->nchan( ifnos[i] ) ) ;
+  }
+  Vector<uInt> mchans( nchans ) ;
+  if ( anyNE( mchans, mchans[0] ) ) {
+    throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
+  }
+
+  // check if vector size is equal to nchan
+  Vector<Float> fact( val ) ;
+  if ( fact.nelements() != mchans[0] ) {
+    throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
+  }
+
+  // check divided by zero
+  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
+    throw( AipsError("Divided by zero is not recommended." ) ) ;
+  }
+
+  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" || mode == "DIV") {
+      if (mode == "DIV") fact = (float)1.0 / fact;
+      spec *= fact;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts *= fact;
+        tsysCol.put(i, ts);
+      }
+    } else if ( mode == "ADD"  || mode == "SUB") {
+      if (mode == "SUB") fact *= (float)-1.0 ;
+      spec += fact;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts += fact;
+        tsysCol.put(i, ts);
+      }
+    }
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
+                                                 const std::vector<float> val,
+                                                 const std::string& mode,
+                                                 bool tsys )
+{
+  if ( val.size() == 1 ) {
+    return unaryOperate( in, val[0], mode, tsys ) ;
+  }
+
+  // conformity of SPECTRA and TSYS
+  if ( tsys ) {
+    TableIterator titer(in->table(), "IFNO");
+    while ( !titer.pastEnd() ) {
+      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
+      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
+      Array<Float> spec = specCol.getColumn() ;
+      Array<Float> ts = tsysCol.getColumn() ;
+      if ( !spec.conform( ts ) ) {
+        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
+      }
+      titer.next() ;
+    }
+  }
+
+  // check if vector size is equal to nrow
+  Vector<Float> fact( val ) ;
+  if ( fact.nelements() != in->nrow() ) {
+    throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
+  }
+
+  // check divided by zero
+  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
+    throw( AipsError("Divided by zero is not recommended." ) ) ;
+  }
+
+  CountedPtr< Scantable > out = getScantable(in, false);
+  Table& tab = out->table();
+  ArrayColumn<Float> specCol(tab,"SPECTRA");
+  ArrayColumn<Float> tsysCol(tab,"TSYS");
+  if (mode == "DIV") fact = (float)1.0 / fact;
+  if (mode == "SUB") fact *= (float)-1.0 ;
+  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" || mode == "DIV") {
+      spec *= fact[i];
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts *= fact[i];
+        tsysCol.put(i, ts);
+      }
+    } else if ( mode == "ADD"  || mode == "SUB") {
+      spec += fact[i];
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts += fact[i];
+        tsysCol.put(i, ts);
+      }
+    }
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
+                                                const std::vector< std::vector<float> > val,
+                                                const std::string& mode,
+                                                bool tsys )
+{
+  // conformity of SPECTRA and TSYS
+  if ( tsys ) {
+    TableIterator titer(in->table(), "IFNO");
+    while ( !titer.pastEnd() ) {
+      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
+      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
+      Array<Float> spec = specCol.getColumn() ;
+      Array<Float> ts = tsysCol.getColumn() ;
+      if ( !spec.conform( ts ) ) {
+        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
+      }
+      titer.next() ;
+    }
+  }
+
+  // some checks 
+  vector<uInt> nchans;
+  for ( uInt i = 0 ; i < in->nrow() ; i++ ) {
+    nchans.push_back( (in->getSpectrum( i )).size() ) ;
+  }
+  //Vector<uInt> mchans( nchans ) ;
+  vector< Vector<Float> > facts ;
+  for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
+    Vector<Float> tmp( val[i] ) ;
+    // check divided by zero
+    if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
+      throw( AipsError("Divided by zero is not recommended." ) ) ;
+    }
+    // conformity check
+    if ( tmp.nelements() != nchans[i] ) {
+      stringstream ss ;
+      ss << "Row " << i << ": Vector size must be same as number of channel." ;
+      throw( AipsError( ss.str() ) ) ;
+    }
+    facts.push_back( tmp ) ;
+  }
+
+
+  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> fact = facts[i] ;
+    Vector<Float> spec;
+    Vector<Float> ts;
+    specCol.get(i, spec);
+    tsysCol.get(i, ts);
+    if (mode == "MUL" || mode == "DIV") {
+      if (mode == "DIV") fact = (float)1.0 / fact;
+      spec *= fact;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts *= fact;
+        tsysCol.put(i, ts);
+      }
+    } else if ( mode == "ADD"  || mode == "SUB") {
+      if (mode == "SUB") fact *= (float)-1.0 ;
+      spec += fact;
+      specCol.put(i, spec);
+      if ( tsys ) {
+        ts += fact;
         tsysCol.put(i, ts);
       }
@@ -386 +796 @@
 }
 
+MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
+                                         const Vector<uChar>& f)
+{
+  Vector<Bool> mask;
+  mask.resize(f.shape());
+  convertArray(mask, f);
+  return MaskedArray<Double>(s,!mask);
+}
+
 Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
 {
@@ -402 +821 @@
   // 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));
+  Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
+  Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
   if ( offs.nrow() == 0 )
     throw(AipsError("No 'off' scans present."));
@@ -641 +1060 @@
          //Debug
          //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
+         //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ; 
+         //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
          meanoff = sum(spoff)/noff;
          meandiff = sum(spdiff)/ndiff;
@@ -760 +1181 @@
       //Debug
       //cerr<<"Tsys used="<<tsysrefscalar<<endl;
+      //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
+      //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
       // fill the result, replay signal tsys by reference tsys
       outintCol.put(i, resint);
@@ -782 +1205 @@
   setInsitu(false);
   STSelector sel;
-  std::vector<int> scan1, scan2, beams;
+  std::vector<int> scan1, scan2, beams, types;
   std::vector< vector<int> > scanpair;
-  std::vector<string> calstate;
+  //std::vector<string> calstate;
+  std::vector<int> calstate;
   String msg;
 
@@ -831 +1255 @@
   scanpair.push_back(scan1);
   scanpair.push_back(scan2);
-  calstate.push_back("*calon");
-  calstate.push_back("*[^calon]");
+  //calstate.push_back("*calon");
+  //calstate.push_back("*[^calon]");
+  calstate.push_back(SrcType::NODCAL);
+  calstate.push_back(SrcType::NOD);
   CountedPtr< Scantable > ws = getScantable(s, false);
   uInt l=0;
@@ -841 +1267 @@
           sel.reset();
           sel.setScans(scanpair[i]);
-          sel.setName(calstate[k]);
+          //sel.setName(calstate[k]);
+          types.clear();
+          types.push_back(calstate[k]);
+          sel.setTypes(types);
           beams.clear();
           beams.push_back(j);
@@ -926 +1355 @@
       //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
       // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
+      // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
+      // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
       tsys1[0] = sqrt(tsyssq1 + tsyssq2);
       Array<Float> avtsys =  tsys1;
@@ -952 +1383 @@
   CountedPtr< Scantable > ws = getScantable(s, false);
   CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
-  CountedPtr< Scantable > calsig, calref, out;
+  CountedPtr< Scantable > calsig, calref, out, out1, out2;
+  Bool nofold=False;
+  vector<int> types ; 
 
   //split the data
-  sel.setName("*_fs");
+  //sel.setName("*_fs");
+  types.push_back( SrcType::FSON ) ;
+  sel.setTypes( types ) ;
   ws->setSelection(sel);
   sig = getScantable(ws,false);
   sel.reset();
-  sel.setName("*_fs_calon");
+  types.clear() ;
+  //sel.setName("*_fs_calon");
+  types.push_back( SrcType::FONCAL ) ;
+  sel.setTypes( types ) ;
   ws->setSelection(sel);
   sigwcal = getScantable(ws,false);
   sel.reset();
-  sel.setName("*_fsr");
+  types.clear() ;
+  //sel.setName("*_fsr");
+  types.push_back( SrcType::FSOFF ) ;
+  sel.setTypes( types ) ;
   ws->setSelection(sel);
   ref = getScantable(ws,false);
   sel.reset();
-  sel.setName("*_fsr_calon");
+  types.clear() ;
+  //sel.setName("*_fsr_calon");
+  types.push_back( SrcType::FOFFCAL ) ;
+  sel.setTypes( types ) ;
   ws->setSelection(sel);
   refwcal = getScantable(ws,false);
+  sel.reset() ;
+  types.clear() ;
 
   calsig = dototalpower(sigwcal, sig, tcal=tcal);
   calref = dototalpower(refwcal, ref, tcal=tcal);
 
-  out=dosigref(calsig,calref,smoothref,tsysv,tau);
-
+  out1=dosigref(calsig,calref,smoothref,tsysv,tau); 
+  out2=dosigref(calref,calsig,smoothref,tsysv,tau); 
+
+  Table& tabout1=out1->table();
+  Table& tabout2=out2->table();
+  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
+  ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
+  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
+  Vector<Float> spec; specCol.get(0, spec);
+  uInt nchan = spec.nelements();
+  uInt freqid1; freqidCol1.get(0,freqid1);
+  uInt freqid2; freqidCol2.get(0,freqid2);
+  Double rp1, rp2, rv1, rv2, inc1, inc2;
+  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
+  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
+  //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
+  //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
+  //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
+  if (rp1==rp2) {
+    Double foffset = rv1 - rv2;
+    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
+    if (choffset >= nchan) {
+      //cerr<<"out-band frequency switching, no folding"<<endl;
+      LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
+      os<<"out-band frequency switching, no folding"<<LogIO::POST;
+      nofold = True;
+    }
+  }
+
+  if (nofold) {
+    std::vector< CountedPtr< Scantable > > tabs; 
+    tabs.push_back(out1);
+    tabs.push_back(out2);
+    out = merge(tabs);
+  }
+  else {
+    //out = out1;
+    Double choffset = ( rv1 - rv2 ) / inc2 ;
+    out = dofold( out1, out2, choffset ) ;
+  }
+    
   return out;
+}
+
+CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
+                                      const CountedPtr<Scantable> &ref,
+                                      Double choffset,
+                                      Double choffset2 )
+{
+  LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
+  os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
+
+  // output scantable
+  CountedPtr<Scantable> out = getScantable( sig, false ) ;
+
+  // separate choffset to integer part and decimal part
+  Int ioffset = (Int)choffset ;
+  Double doffset = choffset - ioffset ;
+  Int ioffset2 = (Int)choffset2 ;
+  Double doffset2 = choffset2 - ioffset2 ;
+  os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
+  os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ;  
+
+  // get column
+  ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
+  ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
+  ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
+  ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
+  ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
+  ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
+  ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
+  ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
+  ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
+  ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
+
+  // check
+  if ( ioffset == 0 ) {
+    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
+    os << "channel offset is zero, no folding" << LogIO::POST ;
+    return out ;
+  }
+  int nchan = ref->nchan() ;
+  if ( abs(ioffset) >= nchan ) {
+    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
+    os << "out-band frequency switching, no folding" << LogIO::POST ;
+    return out ;
+  }
+
+  // attach column for output scantable
+  ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
+  ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
+  ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
+  ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
+  ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
+  ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
+
+  // for each row
+  // assume that the data order are same between sig and ref
+  RowAccumulator acc( asap::W_TINTSYS ) ; 
+  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
+    // get values
+    Vector<Float> spsig ;
+    specCol1.get( i, spsig ) ;
+    Vector<Float> spref ;
+    specCol2.get( i, spref ) ;
+    Vector<Float> tsyssig ;
+    tsysCol1.get( i, tsyssig ) ;
+    Vector<Float> tsysref ;
+    tsysCol2.get( i, tsysref ) ;
+    Vector<uChar> flagsig ;
+    flagCol1.get( i, flagsig ) ;
+    Vector<uChar> flagref ;
+    flagCol2.get( i, flagref ) ;
+    Double timesig ;
+    mjdCol1.get( i, timesig ) ;
+    Double timeref ;
+    mjdCol2.get( i, timeref ) ;
+    Double intsig ;
+    intervalCol1.get( i, intsig ) ;
+    Double intref ;
+    intervalCol2.get( i, intref ) ;
+
+    // shift reference spectra
+    int refchan = spref.nelements() ;
+    Vector<Float> sspref( spref.nelements() ) ;
+    Vector<Float> stsysref( tsysref.nelements() ) ;
+    Vector<uChar> sflagref( flagref.nelements() ) ;
+    if ( ioffset > 0 ) {
+      // SPECTRA and FLAGTRA
+      for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
+        sspref[j] = spref[j+ioffset] ; 
+        sflagref[j] = flagref[j+ioffset] ;
+      }
+      for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
+        sspref[j] = spref[j-refchan+ioffset] ;
+        sflagref[j] = flagref[j-refchan+ioffset] ;
+      }
+      spref = sspref.copy() ;
+      flagref = sflagref.copy() ;
+      for ( int j = 0 ; j < refchan - 1 ; j++ ) {
+        sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
+        sflagref[j] = flagref[j+1] + flagref[j] ;
+      }
+      sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
+      sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ; 
+
+      // TSYS
+      if ( spref.nelements() == tsysref.nelements() ) {
+        for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
+          stsysref[j] = tsysref[j+ioffset] ;
+        }
+        for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
+          stsysref[j] = tsysref[j-refchan+ioffset] ;
+        }
+        tsysref = stsysref.copy() ;
+        for ( int j = 0 ; j < refchan - 1 ; j++ ) {
+          stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
+        }
+        stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
+      }
+    }
+    else {
+      // SPECTRA and FLAGTRA
+      for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
+        sspref[j] = spref[refchan+ioffset+j] ;
+        sflagref[j] = flagref[refchan+ioffset+j] ;
+      }
+      for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
+        sspref[j] = spref[j+ioffset] ;
+        sflagref[j] = flagref[j+ioffset] ;
+      }
+      spref = sspref.copy() ;
+      flagref = sflagref.copy() ;
+      sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
+      sflagref[0] = flagref[0] + flagref[refchan-1] ; 
+      for ( int j = 1 ; j < refchan ; j++ ) {
+        sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
+        sflagref[j] = flagref[j-1] + flagref[j] ;
+      }
+      // TSYS
+      if ( spref.nelements() == tsysref.nelements() ) {
+        for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
+          stsysref[j] = tsysref[refchan+ioffset+j] ;
+        }
+        for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
+          stsysref[j] = tsysref[j+ioffset] ;
+        }
+        tsysref = stsysref.copy() ;
+        stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
+        for ( int j = 1 ; j < refchan ; j++ ) {
+          stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
+        }
+      }
+    }
+
+    // shift signal spectra if necessary (only for APEX?)
+    if ( choffset2 != 0.0 ) {
+      int sigchan = spsig.nelements() ;
+      Vector<Float> sspsig( spsig.nelements() ) ;
+      Vector<Float> stsyssig( tsyssig.nelements() ) ;
+      Vector<uChar> sflagsig( flagsig.nelements() ) ;
+      if ( ioffset2 > 0 ) {
+        // SPECTRA and FLAGTRA
+        for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
+          sspsig[j] = spsig[j+ioffset2] ; 
+          sflagsig[j] = flagsig[j+ioffset2] ;
+        }
+        for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
+          sspsig[j] = spsig[j-sigchan+ioffset2] ;
+          sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
+        }
+        spsig = sspsig.copy() ;
+        flagsig = sflagsig.copy() ;
+        for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
+          sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
+          sflagsig[j] = flagsig[j+1] || flagsig[j] ;
+        }
+        sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
+        sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ; 
+        // TSTS
+        if ( spsig.nelements() == tsyssig.nelements() ) {
+          for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
+            stsyssig[j] = tsyssig[j+ioffset2] ;
+          }
+          for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
+            stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
+          }
+          tsyssig = stsyssig.copy() ;
+          for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
+            stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
+          }
+          stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
+        }
+      }
+      else {
+        // SPECTRA and FLAGTRA
+        for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
+          sspsig[j] = spsig[sigchan+ioffset2+j] ;
+          sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
+        }
+        for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
+          sspsig[j] = spsig[j+ioffset2] ;
+          sflagsig[j] = flagsig[j+ioffset2] ;
+        }
+        spsig = sspsig.copy() ;
+        flagsig = sflagsig.copy() ;
+        sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
+        sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ; 
+        for ( int j = 1 ; j < sigchan ; j++ ) {
+          sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
+          sflagsig[j] = flagsig[j-1] + flagsig[j] ;
+        }
+        // TSYS
+        if ( spsig.nelements() == tsyssig.nelements() ) {
+          for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
+            stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
+          }
+          for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
+            stsyssig[j] = tsyssig[j+ioffset2] ;
+          }
+          tsyssig = stsyssig.copy() ;
+          stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
+          for ( int j = 1 ; j < sigchan ; j++ ) {
+            stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
+          }
+        }
+      }
+    }
+
+    // folding
+    acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
+    acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
+    
+    // put result
+    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() ) ;
+    intervalColOut.put( i, acc.getInterval() ) ;
+    mjdColOut.put( i, acc.getTime() ) ;
+    // change FREQ_ID to unshifted IF setting (only for APEX?)
+    if ( choffset2 != 0.0 ) {
+      uInt freqid = fidColOut( 0 ) ; // assume single-IF data
+      double refpix, refval, increment ;
+      out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
+      refval -= choffset * increment ;
+      uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
+      Vector<uInt> freqids = fidColOut.getColumn() ;
+      for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
+        if ( freqids[j] == freqid ) 
+          freqids[j] = newfreqid ;
+      }
+      fidColOut.putColumn( freqids ) ;
+    }
+
+    acc.reset() ;
+  }
+
+  return out ;
 }
 
@@ -1051 +1795 @@
     }
     out.push_back(outstat);
+  }
+  return out;
+}
+
+std::vector< int > STMath::minMaxChan( 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<int> out;
+  for (uInt i=0; i < tab.nrow(); ++i ) {
+    Vector<Float> spec; specCol.get(i, spec);
+    Vector<uChar> flag; flagCol.get(i, flag);
+    MaskedArray<Float> ma  = maskedArray(spec, flag);
+    if (ma.ndim() != 1) {
+      throw (ArrayError(
+          "std::vector<int> STMath::minMaxChan("
+          "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
+          " std::string &which)"
+          " - MaskedArray is not 1D"));
+    }
+    IPosition outpos(1,0);
+    if ( spec.nelements() == m.nelements() ) {
+      outpos = mathutil::minMaxPos(which, ma(m));
+    } else {
+      outpos = mathutil::minMaxPos(which, ma);
+    }
+    out.push_back(outpos[0]);
   }
   return out;
@@ -1571 +2347 @@
     if ( ! (*it)->conformant(*out) ) {
       // non conformant.
-      pushLog(String("Warning: Can't merge scantables as header info differs."));
+      //pushLog(String("Warning: Can't merge scantables as header info differs."));
+      LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
+      os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
     }
     out->appendToHistoryTable((*it)->history());
@@ -1593 +2371 @@
           id = out->frequencies().addEntry(rp, rv, inc);
           freqidcol.put(k,id);
-          String name,fname;Double rf;
+          //String name,fname;Double rf;
+          Vector<String> name,fname;Vector<Double> rf;
           (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
           id = out->molecules().addEntry(rf, name, fname);
@@ -1997 +2776 @@
       int fstart = -1;
       int fend = -1;
-      for (int k=0; k < flag.nelements(); ++k ) {
+      for (unsigned int k=0; k < flag.nelements(); ++k ) {
         if (flag[k] > 0) {
           fstart = k;
@@ -2051 +2830 @@
   return out;
 }
+
+// Averaging spectra with different channel/resolution
+CountedPtr<Scantable>
+STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
+                     const bool& compel,
+                     const std::vector<bool>& mask,
+                     const std::string& weight,
+                     const std::string& avmode )
+  throw ( casa::AipsError )
+{
+  LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
+  if ( avmode == "SCAN" && in.size() != 1 )
+    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
+                    "Use merge first."));
+  
+  // check if OTF observation
+  String obstype = in[0]->getHeader().obstype ;
+  Double tol = 0.0 ;
+  if ( obstype.find( "OTF" ) != String::npos ) {
+    tol = TOL_OTF ;
+  }
+  else {
+    tol = TOL_POINT ;
+  }
+
+  CountedPtr<Scantable> out ;     // processed result 
+  if ( compel ) {
+    std::vector< CountedPtr<Scantable> > newin ; // input for average process
+    uInt insize = in.size() ;    // number of input scantables 
+
+    // TEST: do normal average in each table before IF grouping
+    os << "Do preliminary averaging" << LogIO::POST ;
+    vector< CountedPtr<Scantable> > tmpin( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
+      tmpin[itable] = average( v, mask, weight, avmode ) ;
+    }
+
+    // warning
+    os << "Average spectra with different spectral resolution" << LogIO::POST ;
+
+    // temporarily set coordinfo
+    vector<string> oldinfo( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = in[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // columns
+    ScalarColumn<uInt> freqIDCol ;
+    ScalarColumn<uInt> ifnoCol ;
+    ScalarColumn<uInt> scannoCol ;
+
+
+    // check IF frequency coverage
+    // freqid: list of FREQ_ID, which is used, in each table  
+    // iffreq: list of minimum and maximum frequency for each FREQ_ID in 
+    //         each table
+    // freqid[insize][numIF]
+    // freqid: [[id00, id01, ...],
+    //          [id10, id11, ...],
+    //          ...
+    //          [idn0, idn1, ...]]
+    // iffreq[insize][numIF*2]
+    // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
+    //          [min_id10, max_id10, min_id11, max_id11, ...],
+    //          ...
+    //          [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
+    //os << "Check IF settings in each table" << LogIO::POST ;
+    vector< vector<uInt> > freqid( insize );
+    vector< vector<double> > iffreq( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = tmpin[itable]->nrow() ;
+      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        if ( freqid[itable].size() == freqnrows ) {
+          break ;
+        }
+        else {
+          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
+          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
+          uInt id = freqIDCol( irow ) ;
+          if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
+            //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
+            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
+            freqid[itable].push_back( id ) ;
+            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
+            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
+          }
+        }
+      }
+    }
+
+    // debug
+    //os << "IF settings summary:" << endl ;
+    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
+    //os << "   Table" << i << endl ;
+    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
+    //os << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
+    //}
+    //}
+    //os << endl ;
+    //os.post() ;
+
+    // IF grouping based on their frequency coverage
+    // ifgrp: list of table index and FREQ_ID for all members in each IF group
+    // ifgfreq: list of minimum and maximum frequency in each IF group
+    // ifgrp[numgrp][nummember*2]
+    // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
+    //         [table10, freqrow10, table11, freqrow11, ...],
+    //         ...
+    //         [tablen0, freqrown0, tablen1, freqrown1, ...]]
+    // ifgfreq[numgrp*2] 
+    // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
+    //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
+    vector< vector<uInt> > ifgrp ;
+    vector<double> ifgfreq ;
+
+    // parameter for IF grouping
+    // groupmode = OR    retrieve all region 
+    //             AND   only retrieve overlaped region
+    //string groupmode = "AND" ;
+    string groupmode = "OR" ;
+    uInt sizecr = 0 ;
+    if ( groupmode == "AND" ) 
+      sizecr = 2 ;
+    else if ( groupmode == "OR" ) 
+      sizecr = 0 ;
+
+    vector<double> sortedfreq ;
+    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
+      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
+        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
+          sortedfreq.push_back( iffreq[i][j] ) ;
+      }
+    }
+    sort( sortedfreq.begin(), sortedfreq.end() ) ;
+    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
+      ifgfreq.push_back( *i ) ;
+      ifgfreq.push_back( *(i+1) ) ;
+    }
+    ifgrp.resize( ifgfreq.size()/2 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
+        double range0 = iffreq[itable][2*iif] ;
+        double range1 = iffreq[itable][2*iif+1] ;
+        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
+          double fmin = max( range0, ifgfreq[2*j] ) ;
+          double fmax = min( range1, ifgfreq[2*j+1] ) ;
+          if ( fmin < fmax ) {
+            ifgrp[j].push_back( itable ) ;
+            ifgrp[j].push_back( freqid[itable][iif] ) ;
+          }
+        }
+      }
+    }
+    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ; 
+    vector<double>::iterator giter = ifgfreq.begin() ;
+    while( fiter != ifgrp.end() ) {
+      if ( fiter->size() <= sizecr ) {
+        fiter = ifgrp.erase( fiter ) ;
+        giter = ifgfreq.erase( giter ) ;
+        giter = ifgfreq.erase( giter ) ;
+      }
+      else {
+        fiter++ ;
+        advance( giter, 2 ) ;
+      }
+    }
+
+    // Grouping continuous IF groups (without frequency gap)
+    // freqgrp: list of IF group indexes in each frequency group
+    // freqrange: list of minimum and maximum frequency in each frequency group
+    // freqgrp[numgrp][nummember]
+    // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
+    //           [ifgrp10, ifgrp11, ifgrp12, ...],
+    //           ...
+    //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
+    // freqrange[numgrp*2]
+    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
+    vector< vector<uInt> > freqgrp ;
+    double freqrange = 0.0 ;
+    uInt grpnum = 0 ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      // Assumed that ifgfreq was sorted
+      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
+        freqgrp[grpnum-1].push_back( i ) ;
+      }
+      else {
+        vector<uInt> grp0( 1, i ) ;
+        freqgrp.push_back( grp0 ) ;
+        grpnum++ ;
+      }
+      freqrange = ifgfreq[2*i+1] ;
+    }
+        
+
+    // print IF groups
+    ostringstream oss ;
+    oss << "IF Group summary: " << endl ;
+    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+        oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+      }
+      oss << endl ;
+    }
+    oss << endl ;
+    os << oss.str() << LogIO::POST ;
+    
+    // print frequency group
+    oss.str("") ;
+    oss << "Frequency Group summary: " << endl ;
+    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
+    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
+      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
+        oss << freqgrp[i][j] << " " ;
+      }
+      oss << endl ;
+    }
+    oss << endl ;
+    os << oss.str() << LogIO::POST ;
+
+    // membership check
+    // groups: list of IF group indexes whose frequency range overlaps with 
+    //         that of each table and IF
+    // groups[numtable][numIF][nummembership]
+    // groups: [[[grp, grp,...], [grp, grp,...],...],
+    //          [[grp, grp,...], [grp, grp,...],...],
+    //          ...
+    //          [[grp, grp,...], [grp, grp,...],...]] 
+    vector< vector< vector<uInt> > > groups( insize ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      groups[i].resize( freqid[i].size() ) ;
+    }
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+        uInt tableid = ifgrp[igrp][2*imem] ;
+        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
+        if ( iter != freqid[tableid].end() ) {
+          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
+          groups[tableid][rowid].push_back( igrp ) ;
+        }
+      }
+    }
+
+    // print membership
+    //oss.str("") ;
+    //for ( uInt i = 0 ; i < insize ; i++ ) {
+    //oss << "Table " << i << endl ;
+    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
+    //oss << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
+    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
+    //oss << setw( 2 ) << groups[i][j][k] << " " ;
+    //}
+    //oss << endl ;
+    //}
+    //}
+    //os << oss.str() << LogIO::POST ;
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // Create additional table if needed
+    bool oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    vector< vector<uInt> > addrow( insize ) ;
+    vector<uInt> addtable( insize, 0 ) ;
+    vector<uInt> newtableids( insize ) ;
+    vector<uInt> newifids( insize, 0 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      //os << "Table " << itable << ": " ;
+      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
+        //os << addrow[itable][ifrow] << " " ;
+      }
+      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
+      //os << "(" << addtable[itable] << ")" << LogIO::POST ;
+    }
+    newin.resize( insize ) ;
+    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      newtableids[i] = i ;
+    }
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
+        vector<int> freqidlist ;
+        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
+          if ( groups[itable][i].size() > iadd + 1 ) {
+            freqidlist.push_back( freqid[itable][i] ) ;
+          }
+        }
+        stringstream taqlstream ;
+        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
+        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
+          taqlstream << i ;
+          if ( i < freqidlist.size() - 1 )
+            taqlstream << "," ;
+          else
+            taqlstream << "]" ;
+        }
+        string taql = taqlstream.str() ;
+        //os << "taql = " << taql << LogIO::POST ;
+        STSelector selector = STSelector() ;
+        selector.setTaQL( taql ) ;
+        add->setSelection( selector ) ;
+        newin.push_back( add ) ;
+        newtableids.push_back( itable ) ;
+        newifids.push_back( iadd + 1 ) ;
+      }
+    }
+
+    // udpate ifgrp
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+          if ( groups[itable][ifrow].size() > iadd + 1 ) {
+            uInt igrp = groups[itable][ifrow][iadd+1] ;
+            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
+                ifgrp[igrp][2*imem] = insize + iadd ;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // print IF groups again for debug
+    //oss.str( "" ) ;
+    //oss << "IF Group summary: " << endl ;
+    //oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+    //oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+    //oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+    //}
+    //oss << endl ;
+    //}
+    //oss << endl ;
+    //os << oss.str() << LogIO::POST ;
+
+    // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation 
+    os << "All scan number is set to 0" << LogIO::POST ;
+    //os << "All IF number is set to IF group index" << LogIO::POST ;
+    insize = newin.size() ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      Table &tmpt = newin[itable]->table() ;
+      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
+      scannoCol.attach( tmpt, "SCANNO" ) ;
+      ifnoCol.attach( tmpt, "IFNO" ) ;
+      for ( uInt irow=0 ; irow < rows ; irow++ ) {
+        scannoCol.put( irow, 0 ) ;
+        uInt freqID = freqIDCol( irow ) ;
+        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
+        if ( iter != freqid[newtableids[itable]].end() ) {
+          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
+          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
+        }
+        else {
+          throw(AipsError("IF grouping was wrong in additional tables.")) ;
+        } 
+      }
+    }
+    oldinfo.resize( insize ) ;
+    setInsitu( oldInsitu ) ;
+
+    // temporarily set coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // save column values in the vector
+    vector< vector<uInt> > freqTableIdVec( insize ) ;
+    vector< vector<uInt> > freqIdVec( insize ) ;
+    vector< vector<uInt> > ifNoVec( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      ScalarColumn<uInt> freqIDs ;
+      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
+      ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
+      freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
+      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
+        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
+      }
+      for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
+        freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
+        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
+      }
+    }
+
+    // reset spectra and flagtra: pick up common part of frequency coverage
+    //os << "Pick common frequency range and align resolution" << LogIO::POST ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      int nminchan = -1 ;
+      int nmaxchan = -1 ;
+      vector<uInt> freqIdUpdate ; 
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
+        double minfreq = ifgfreq[2*ifno] ;
+        double maxfreq = ifgfreq[2*ifno+1] ;
+        //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ;
+        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
+        int nchan = abcissa.size() ;
+        double resol = abcissa[1] - abcissa[0] ;
+        //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
+        if ( minfreq <= abcissa[0] ) 
+          nminchan = 0 ;
+        else {
+          //double cfreq = ( minfreq - abcissa[0] ) / resol ;
+          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
+          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
+        }
+        if ( maxfreq >= abcissa[abcissa.size()-1] )
+          nmaxchan = abcissa.size() - 1 ;
+        else {
+          //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
+          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
+          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
+        }
+        //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
+        if ( nmaxchan > nminchan ) {
+          newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
+          int newchan = nmaxchan - nminchan + 1 ;
+          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
+            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
+        }
+        else {
+          throw(AipsError("Failed to pick up common part of frequency range.")) ;
+        }
+      }
+      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
+        uInt freqId = freqIdUpdate[i] ;
+        Double refpix ;
+        Double refval ;
+        Double increment ;
+        
+        // update row
+        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
+        refval = refval - ( refpix - nminchan ) * increment ;
+        refpix = 0 ;
+        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
+      }    
+    }
+
+    
+    // reset spectra and flagtra: align spectral resolution
+    //os << "Align spectral resolution" << LogIO::POST ;
+    // gmaxdnu: the coarsest frequency resolution in the frequency group
+    // gmemid: member index that have a resolution equal to gmaxdnu
+    // gmaxdnu[numfreqgrp]
+    // gmaxdnu: [dnu0, dnu1, ...]
+    // gmemid[numfreqgrp]
+    // gmemid: [id0, id1, ...]
+    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
+    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution 
+      int minchan = INT_MAX ;     // minimum channel number
+      Double refpixref = -1 ;     // reference of 'reference pixel'
+      Double refvalref = -1 ;     // reference of 'reference frequency'
+      Double refinc = -1 ;        // reference frequency resolution 
+      uInt refreqid ;
+      uInt reftable = INT_MAX;
+      // process only if group member > 1
+      if ( ifgrp[igrp].size() > 2 ) {
+        // find minchan and maxdnu in each group
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+          if ( iter != freqIdVec[tableid].end() ) {
+            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+            int nchan = abcissa.size() ;
+            double dnu = abcissa[1] - abcissa[0] ;
+            //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ; 
+            if ( nchan < minchan ) {
+              minchan = nchan ;
+              maxdnu = dnu ;
+              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
+              refreqid = rowid ;
+              reftable = tableid ;
+            }
+          }
+        }
+        // regrid spectra in each group
+        os << "GROUP " << igrp << endl ;
+        os << "   Channel number is adjusted to " << minchan << endl ;
+        os << "   Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ;
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
+          //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ;
+          //os << "   regridChannel applied to " ;
+          if ( tableid != reftable ) 
+            refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
+          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
+            uInt tfreqid = freqIdVec[tableid][irow] ;
+            if ( tfreqid == rowid ) {     
+              //os << irow << " " ;
+              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
+              freqIDCol.put( irow, refreqid ) ;
+              freqIdVec[tableid][irow] = refreqid ;
+            }
+          }
+          //os << LogIO::POST ;
+        } 
+      }
+      else {
+        uInt tableid = ifgrp[igrp][0] ;
+        uInt rowid = ifgrp[igrp][1] ;
+        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+        if ( iter != freqIdVec[tableid].end() ) {
+          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+          minchan = abcissa.size() ;
+          maxdnu = abcissa[1] - abcissa[0] ;
+        }
+      }
+      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
+          if ( maxdnu > gmaxdnu[i] ) {
+            gmaxdnu[i] = maxdnu ;
+            gmemid[i] = igrp ;
+          }
+          break ;
+        }
+      }
+    }
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }      
+
+    // accumulate all rows into the first table 
+    // NOTE: assumed in.size() = 1
+    vector< CountedPtr<Scantable> > tmp( 1 ) ;
+    if ( newin.size() == 1 )
+      tmp[0] = newin[0] ;
+    else 
+      tmp[0] = merge( newin ) ;
+
+    //return tmp[0] ;
+
+    // average
+    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
+
+    //return tmpout ;
+
+    // combine frequency group
+    os << "Combine spectra based on frequency grouping" << LogIO::POST ;
+    os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
+    vector<string> coordinfo = tmpout->getCoordInfo() ;
+    oldinfo[0] = coordinfo[0] ;
+    coordinfo[0] = "Hz" ;
+    tmpout->setCoordInfo( coordinfo ) ;
+    // create proformas of output table
+    stringstream taqlstream ;
+    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
+    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
+      taqlstream << gmemid[i] ;
+      if ( i < gmemid.size() - 1 )
+        taqlstream << "," ;
+      else
+        taqlstream << "]" ;
+    }
+    string taql = taqlstream.str() ;
+    //os << "taql = " << taql << LogIO::POST ;
+    STSelector selector = STSelector() ;
+    selector.setTaQL( taql ) ;
+    oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    out = getScantable( tmpout, false ) ;
+    setInsitu( oldInsitu ) ;
+    out->setSelection( selector ) ;
+    // regrid rows
+    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
+    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
+      uInt ifno = ifnoCol( irow ) ;
+      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
+          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
+          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
+          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
+          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
+          break ;
+        }
+      }
+    }
+    // combine spectra
+    ArrayColumn<Float> specColOut ;
+    specColOut.attach( out->table(), "SPECTRA" ) ;
+    ArrayColumn<uChar> flagColOut ;
+    flagColOut.attach( out->table(), "FLAGTRA" ) ;
+    ScalarColumn<uInt> ifnoColOut ;
+    ifnoColOut.attach( out->table(), "IFNO" ) ;
+    ScalarColumn<uInt> polnoColOut ;
+    polnoColOut.attach( out->table(), "POLNO" ) ;
+    ScalarColumn<uInt> freqidColOut ;
+    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
+    MDirection::ScalarColumn dirColOut ;
+    dirColOut.attach( out->table(), "DIRECTION" ) ;
+    Table &tab = tmpout->table() ;
+    Block<String> cols(1);
+    cols[0] = String("POLNO") ;
+    TableIterator iter( tab, cols ) ;
+    bool done = false ;
+    vector< vector<uInt> > sizes( freqgrp.size() ) ;
+    while( !iter.pastEnd() ) {
+      vector< vector<Float> > specout( freqgrp.size() ) ;
+      vector< vector<uChar> > flagout( freqgrp.size() ) ;
+      ArrayColumn<Float> specCols ;
+      specCols.attach( iter.table(), "SPECTRA" ) ;
+      ArrayColumn<uChar> flagCols ;
+      flagCols.attach( iter.table(), "FLAGTRA" ) ;
+      ifnoCol.attach( iter.table(), "IFNO" ) ;
+      ScalarColumn<uInt> polnos ;
+      polnos.attach( iter.table(), "POLNO" ) ;
+      MDirection::ScalarColumn dircol ;
+      dircol.attach( iter.table(), "DIRECTION" ) ;
+      uInt polno = polnos( 0 ) ;
+      //os << "POLNO iteration: " << polno << LogIO::POST ;
+//       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+//      sizes[igrp].resize( freqgrp[igrp].size() ) ;
+//      for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
+//        for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
+//          uInt ifno = ifnoCol( irow ) ;
+//          if ( ifno == freqgrp[igrp][imem] ) {
+//            Vector<Float> spec = specCols( irow ) ;
+//            Vector<uChar> flag = flagCols( irow ) ;
+//            vector<Float> svec ;
+//            spec.tovector( svec ) ;
+//            vector<uChar> fvec ;
+//            flag.tovector( fvec ) ;
+//            //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
+//            specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
+//            flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
+//            //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
+//            sizes[igrp][imem] = spec.nelements() ;
+//          }
+//        }
+//      }
+//      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
+//        uInt ifout = ifnoColOut( irow ) ;
+//        uInt polout = polnoColOut( irow ) ;
+//        if ( ifout == gmemid[igrp] && polout == polno ) {
+//          // set SPECTRA and FRAGTRA
+//          Vector<Float> newspec( specout[igrp] ) ;
+//          Vector<uChar> newflag( flagout[igrp] ) ;
+//          specColOut.put( irow, newspec ) ;
+//          flagColOut.put( irow, newflag ) ;
+//          // IFNO renumbering
+//          ifnoColOut.put( irow, igrp ) ;
+//        }
+//      } 
+//       }
+      // get a list of number of channels for each frequency group member
+      if ( !done ) {
+        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+          sizes[igrp].resize( freqgrp[igrp].size() ) ;
+          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
+            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
+              uInt ifno = ifnoCol( irow ) ;
+              if ( ifno == freqgrp[igrp][imem] ) {
+                Vector<Float> spec = specCols( irow ) ;
+                sizes[igrp][imem] = spec.nelements() ;
+                break ;
+              }                
+            }
+          }
+        }
+        done = true ;
+      }
+      // combine spectra
+      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
+        uInt polout = polnoColOut( irow ) ;
+        if ( polout == polno ) {
+          uInt ifout = ifnoColOut( irow ) ;
+          Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
+          uInt igrp ;
+          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
+            if ( ifout == gmemid[jgrp] ) {
+              igrp = jgrp ;
+              break ;
+            }
+          }
+          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
+            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
+              uInt ifno = ifnoCol( jrow ) ;
+              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
+              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
+              Double dx = tdir[0] - direction[0] ;
+              Double dy = tdir[1] - direction[1] ;
+              Double dd = sqrt( dx * dx + dy * dy ) ;
+              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
+              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
+                Vector<Float> spec = specCols( jrow ) ;
+                Vector<uChar> flag = flagCols( jrow ) ;
+                vector<Float> svec ;
+                spec.tovector( svec ) ;
+                vector<uChar> fvec ;
+                flag.tovector( fvec ) ;
+                //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
+                specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
+                flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
+                //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
+              }
+            }
+          }
+          // set SPECTRA and FRAGTRA
+          Vector<Float> newspec( specout[igrp] ) ;
+          Vector<uChar> newflag( flagout[igrp] ) ;
+          specColOut.put( irow, newspec ) ;
+          flagColOut.put( irow, newflag ) ;
+          // IFNO renumbering
+          ifnoColOut.put( irow, igrp ) ; 
+        }
+      }
+      iter++ ;
+    }
+    // update FREQUENCIES subtable
+    vector<bool> updated( freqgrp.size(), false ) ;
+    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+      uInt index = 0 ;
+      uInt pixShift = 0 ;
+      while ( freqgrp[igrp][index] != gmemid[igrp] ) {
+        pixShift += sizes[igrp][index++] ;
+      }
+      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
+        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
+          uInt freqidOut = freqidColOut( irow ) ;
+          //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
+          double refpix ;
+          double refval ;
+          double increm ;
+          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
+          refpix += pixShift ;
+          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
+          updated[igrp] = true ;
+        }
+      }
+    }
+
+    //out = tmpout ;
+
+    coordinfo = tmpout->getCoordInfo() ;
+    coordinfo[0] = oldinfo[0] ;
+    tmpout->setCoordInfo( coordinfo ) ;
+  }
+  else {
+    // simple average
+    out =  average( in, mask, weight, avmode ) ;
+  }
+  
+  return out ;
+}
+
+CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
+                                     const String calmode, 
+                                     const String antname )
+{
+  // frequency switch
+  if ( calmode == "fs" ) {
+    return cwcalfs( s, antname ) ;
+  }
+  else {
+    vector<bool> masks = s->getMask( 0 ) ;
+    vector<int> types ;
+
+    // sky scan
+    STSelector sel = STSelector() ;
+    types.push_back( SrcType::SKY ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
+    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+
+    // hot scan
+    types.push_back( SrcType::HOT ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // cold scan
+    CountedPtr<Scantable> acold ;
+//     types.push_back( SrcType::COLD ) ;
+//     sel.setTypes( types ) ;
+//     s->setSelection( sel ) ;
+//     tmp.clear() ;
+//     tmp.push_back( getScantable( s, false ) ) ;
+//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
+//     s->unsetSelection() ;
+//     sel.reset() ;
+//     types.clear() ;
+
+    // off scan
+    types.push_back( SrcType::PSOFF ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // on scan
+    bool insitu = insitu_ ;
+    insitu_ = false ;
+    CountedPtr<Scantable> out = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::PSON ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( out->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // process each on scan
+    ArrayColumn<Float> tsysCol ;
+    tsysCol.attach( out->table(), "TSYS" ) ;
+    for ( int i = 0 ; i < out->nrow() ; i++ ) {
+      vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
+      out->setSpectrum( sp, i ) ;
+      string reftime = out->getTime( i ) ;
+      vector<int> ii( 1, out->getIF( i ) ) ;
+      vector<int> ib( 1, out->getBeam( i ) ) ;
+      vector<int> ip( 1, out->getPol( i ) ) ;
+      sel.setIFs( ii ) ;
+      sel.setBeams( ib ) ;
+      sel.setPolarizations( ip ) ;
+      asky->setSelection( sel ) ;    
+      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
+      const Vector<Float> Vtsys( sptsys ) ;
+      tsysCol.put( i, Vtsys ) ;
+      asky->unsetSelection() ;
+      sel.reset() ;
+    }
+
+    // flux unit
+    out->setFluxUnit( "K" ) ;
+
+    return out ;
+  }
+}
+  
+CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
+                                       const String calmode )
+{
+  // frequency switch
+  if ( calmode == "fs" ) {
+    return almacalfs( s ) ;
+  }
+  else {
+    vector<bool> masks = s->getMask( 0 ) ;
+    
+    // off scan
+    STSelector sel = STSelector() ;
+    vector<int> types ;
+    types.push_back( SrcType::PSOFF ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    // TODO 2010/01/08 TN
+    // Grouping by time should be needed before averaging.
+    // Each group must have own unique SCANNO (should be renumbered).
+    // See PIPELINE/SDCalibration.py
+    CountedPtr<Scantable> soff = getScantable( s, false ) ;
+    Table ttab = soff->table() ;
+    ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
+    uInt nrow = timeCol.nrow() ;
+    Vector<Double> timeSep( nrow - 1 ) ;
+    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
+      timeSep[i] = timeCol(i+1) - timeCol(i) ;
+    }
+    ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
+    Vector<Double> interval = intervalCol.getColumn() ;
+    interval /= 86400.0 ;
+    ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
+    vector<uInt> glist ;
+    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
+      double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
+      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
+      if ( gap > 1.1 ) {
+        glist.push_back( i ) ;
+      }
+    }
+    Vector<uInt> gaplist( glist ) ;
+    //cout << "gaplist = " << gaplist << endl ;
+    uInt newid = 0 ;
+    for ( uInt i = 0 ; i < nrow ; i++ ) {
+      scanCol.put( i, newid ) ;
+      if ( i == gaplist[newid] ) {
+        newid++ ;
+      }
+    }
+    //cout << "new scancol = " << scanCol.getColumn() << endl ;
+    vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
+    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
+    //cout << "aoff.nrow = " << aoff->nrow() << endl ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // on scan
+    bool insitu = insitu_ ;
+    insitu_ = false ;
+    CountedPtr<Scantable> out = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::PSON ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( out->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // process each on scan
+    ArrayColumn<Float> tsysCol ;
+    tsysCol.attach( out->table(), "TSYS" ) ;
+    for ( int i = 0 ; i < out->nrow() ; i++ ) {
+      vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
+      out->setSpectrum( sp, i ) ;
+    }
+
+    // flux unit
+    out->setFluxUnit( "K" ) ;
+
+    return out ;
+  }
+}
+
+CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
+                                       const String antname )
+{
+  vector<int> types ;
+
+  // APEX calibration mode
+  int apexcalmode = 1 ;
+ 
+  if ( antname.find( "APEX" ) != string::npos ) {
+    // check if off scan exists or not
+    STSelector sel = STSelector() ;
+    //sel.setName( offstr1 ) ;
+    types.push_back( SrcType::FLOOFF ) ;
+    sel.setTypes( types ) ;
+    try {
+      s->setSelection( sel ) ;
+    }
+    catch ( AipsError &e ) {
+      apexcalmode = 0 ;
+    }
+    sel.reset() ;
+  }
+  s->unsetSelection() ;
+  types.clear() ;
+
+  vector<bool> masks = s->getMask( 0 ) ;
+  CountedPtr<Scantable> ssig, sref ;
+  CountedPtr<Scantable> out ;
+
+  if ( antname.find( "APEX" ) != string::npos ) {
+    // APEX calibration
+    // sky scan
+    STSelector sel = STSelector() ;
+    types.push_back( SrcType::FLOSKY ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
+    CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    types.push_back( SrcType::FHISKY ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // hot scan
+    types.push_back( SrcType::FLOHOT ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    types.push_back( SrcType::FHIHOT ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // cold scan
+    CountedPtr<Scantable> acoldlo, acoldhi ;
+//     types.push_back( SrcType::FLOCOLD ) ;
+//     sel.setTypes( types ) ;
+//     s->setSelection( sel ) ;
+//     tmp.clear() ;
+//     tmp.push_back( getScantable( s, false ) ) ;
+//     CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
+//     s->unsetSelection() ;
+//     sel.reset() ;
+//     types.clear() ;
+//     types.push_back( SrcType::FHICOLD ) ;
+//     sel.setTypes( types ) ;
+//     s->setSelection( sel ) ;
+//     tmp.clear() ;
+//     tmp.push_back( getScantable( s, false ) ) ;
+//     CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
+//     s->unsetSelection() ;
+//     sel.reset() ;
+//     types.clear() ;
+
+    // ref scan
+    bool insitu = insitu_ ;
+    insitu_ = false ;
+    sref = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::FSLO ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( sref->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // sig scan
+    insitu_ = false ;
+    ssig = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::FSHI ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( ssig->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;  
+    types.clear() ;
+          
+    if ( apexcalmode == 0 ) {
+      // APEX fs data without off scan
+      // process each sig and ref scan
+      ArrayColumn<Float> tsysCollo ;
+      tsysCollo.attach( ssig->table(), "TSYS" ) ;
+      ArrayColumn<Float> tsysColhi ;
+      tsysColhi.attach( sref->table(), "TSYS" ) ;
+      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
+        vector< CountedPtr<Scantable> > sky( 2 ) ;
+        sky[0] = askylo ;
+        sky[1] = askyhi ;
+        vector< CountedPtr<Scantable> > hot( 2 ) ;
+        hot[0] = ahotlo ;
+        hot[1] = ahothi ;
+        vector< CountedPtr<Scantable> > cold( 2 ) ;
+        //cold[0] = acoldlo ;
+        //cold[1] = acoldhi ;
+        vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
+        ssig->setSpectrum( sp, i ) ;
+        string reftime = ssig->getTime( i ) ;
+        vector<int> ii( 1, ssig->getIF( i ) ) ;
+        vector<int> ib( 1, ssig->getBeam( i ) ) ;
+        vector<int> ip( 1, ssig->getPol( i ) ) ;
+        sel.setIFs( ii ) ;
+        sel.setBeams( ib ) ;
+        sel.setPolarizations( ip ) ;
+        askylo->setSelection( sel ) ;
+        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
+        const Vector<Float> Vtsyslo( sptsys ) ;
+        tsysCollo.put( i, Vtsyslo ) ;
+        askylo->unsetSelection() ;
+        sel.reset() ;
+        sky[0] = askyhi ;
+        sky[1] = askylo ;
+        hot[0] = ahothi ;
+        hot[1] = ahotlo ;
+        cold[0] = acoldhi ;
+        cold[1] = acoldlo ;
+        sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
+        sref->setSpectrum( sp, i ) ;
+        reftime = sref->getTime( i ) ;
+        ii[0] = sref->getIF( i )  ;
+        ib[0] = sref->getBeam( i ) ;
+        ip[0] = sref->getPol( i ) ;
+        sel.setIFs( ii ) ;
+        sel.setBeams( ib ) ;
+        sel.setPolarizations( ip ) ;
+        askyhi->setSelection( sel ) ;    
+        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
+        const Vector<Float> Vtsyshi( sptsys ) ;
+        tsysColhi.put( i, Vtsyshi ) ;
+        askyhi->unsetSelection() ;
+        sel.reset() ;
+      }
+    }
+    else if ( apexcalmode == 1 ) {
+      // APEX fs data with off scan
+      // off scan
+      types.push_back( SrcType::FLOOFF ) ;
+      sel.setTypes( types ) ;
+      s->setSelection( sel ) ;
+      tmp.clear() ;
+      tmp.push_back( getScantable( s, false ) ) ;
+      CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
+      s->unsetSelection() ;
+      sel.reset() ;
+      types.clear() ;
+      types.push_back( SrcType::FHIOFF ) ;
+      sel.setTypes( types ) ;
+      s->setSelection( sel ) ;
+      tmp.clear() ;
+      tmp.push_back( getScantable( s, false ) ) ;
+      CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
+      s->unsetSelection() ;
+      sel.reset() ;
+      types.clear() ;
+      
+      // process each sig and ref scan
+      ArrayColumn<Float> tsysCollo ;
+      tsysCollo.attach( ssig->table(), "TSYS" ) ;
+      ArrayColumn<Float> tsysColhi ;
+      tsysColhi.attach( sref->table(), "TSYS" ) ;
+      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
+        vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
+        ssig->setSpectrum( sp, i ) ;
+        sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ; 
+        string reftime = ssig->getTime( i ) ;
+        vector<int> ii( 1, ssig->getIF( i ) ) ;
+        vector<int> ib( 1, ssig->getBeam( i ) ) ;
+        vector<int> ip( 1, ssig->getPol( i ) ) ;
+        sel.setIFs( ii ) ;
+        sel.setBeams( ib ) ;
+        sel.setPolarizations( ip ) ;
+        askylo->setSelection( sel ) ;
+        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
+        const Vector<Float> Vtsyslo( sptsys ) ;
+        tsysCollo.put( i, Vtsyslo ) ;
+        askylo->unsetSelection() ;
+        sel.reset() ;
+        sref->setSpectrum( sp, i ) ;
+        reftime = sref->getTime( i ) ;
+        ii[0] = sref->getIF( i )  ;
+        ib[0] = sref->getBeam( i ) ;
+        ip[0] = sref->getPol( i ) ;
+        sel.setIFs( ii ) ;
+        sel.setBeams( ib ) ;
+        sel.setPolarizations( ip ) ;
+        askyhi->setSelection( sel ) ;    
+        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
+        const Vector<Float> Vtsyshi( sptsys ) ;
+        tsysColhi.put( i, Vtsyshi ) ;
+        askyhi->unsetSelection() ;
+        sel.reset() ;
+      }
+    }
+  }
+  else {
+    // non-APEX fs data
+    // sky scan
+    STSelector sel = STSelector() ;
+    types.push_back( SrcType::SKY ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
+    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // hot scan
+    types.push_back( SrcType::HOT ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    tmp.clear() ;
+    tmp.push_back( getScantable( s, false ) ) ;
+    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+
+    // cold scan
+    CountedPtr<Scantable> acold ;
+//     types.push_back( SrcType::COLD ) ;
+//     sel.setTypes( types ) ;
+//     s->setSelection( sel ) ;
+//     tmp.clear() ;
+//     tmp.push_back( getScantable( s, false ) ) ;
+//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
+//     s->unsetSelection() ;
+//     sel.reset() ;
+//     types.clear() ;
+   
+    // ref scan
+    bool insitu = insitu_ ;
+    insitu_ = false ;
+    sref = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::FSOFF ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( sref->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+    
+    // sig scan
+    insitu_ = false ;
+    ssig = getScantable( s, true ) ;
+    insitu_ = insitu ;
+    types.push_back( SrcType::FSON ) ;
+    sel.setTypes( types ) ;
+    s->setSelection( sel ) ;
+    TableCopy::copyRows( ssig->table(), s->table() ) ;
+    s->unsetSelection() ;
+    sel.reset() ;
+    types.clear() ;
+
+    // process each sig and ref scan
+    ArrayColumn<Float> tsysColsig ;
+    tsysColsig.attach( ssig->table(), "TSYS" ) ;
+    ArrayColumn<Float> tsysColref ;
+    tsysColref.attach( ssig->table(), "TSYS" ) ;
+    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
+      vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
+      ssig->setSpectrum( sp, i ) ;
+      string reftime = ssig->getTime( i ) ;
+      vector<int> ii( 1, ssig->getIF( i ) ) ;
+      vector<int> ib( 1, ssig->getBeam( i ) ) ;
+      vector<int> ip( 1, ssig->getPol( i ) ) ;
+      sel.setIFs( ii ) ;
+      sel.setBeams( ib ) ;
+      sel.setPolarizations( ip ) ;
+      asky->setSelection( sel ) ;
+      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
+      const Vector<Float> Vtsys( sptsys ) ;
+      tsysColsig.put( i, Vtsys ) ;
+      asky->unsetSelection() ;
+      sel.reset() ;
+      sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
+      sref->setSpectrum( sp, i ) ;
+      tsysColref.put( i, Vtsys ) ;
+    }
+  }
+
+  // do folding if necessary
+  Table sigtab = ssig->table() ;
+  Table reftab = sref->table() ;
+  ScalarColumn<uInt> sigifnoCol ;
+  ScalarColumn<uInt> refifnoCol ;
+  ScalarColumn<uInt> sigfidCol ;
+  ScalarColumn<uInt> reffidCol ;
+  Int nchan = (Int)ssig->nchan() ;
+  sigifnoCol.attach( sigtab, "IFNO" ) ;
+  refifnoCol.attach( reftab, "IFNO" ) ;
+  sigfidCol.attach( sigtab, "FREQ_ID" ) ;
+  reffidCol.attach( reftab, "FREQ_ID" ) ;
+  Vector<uInt> sfids( sigfidCol.getColumn() ) ;
+  Vector<uInt> rfids( reffidCol.getColumn() ) ;
+  vector<uInt> sfids_unique ;
+  vector<uInt> rfids_unique ;
+  vector<uInt> sifno_unique ;
+  vector<uInt> rifno_unique ;
+  for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
+    if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
+      sfids_unique.push_back( sfids[i] ) ;
+      sifno_unique.push_back( ssig->getIF( i ) ) ;
+    }
+    if ( count( rfids_unique.begin(), rfids_unique.end(),  rfids[i] ) == 0 ) {
+      rfids_unique.push_back( rfids[i] ) ;
+      rifno_unique.push_back( sref->getIF( i ) ) ;
+    }
+  }
+  double refpix_sig, refval_sig, increment_sig ;
+  double refpix_ref, refval_ref, increment_ref ;
+  vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
+  for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
+    ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
+    sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
+    if ( refpix_sig == refpix_ref ) {
+      double foffset = refval_ref - refval_sig ;
+      int choffset = static_cast<int>(foffset/increment_sig) ;
+      double doffset = foffset / increment_sig ;
+      if ( abs(choffset) >= nchan ) {
+        LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
+        os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
+        os << "Just return signal data" << LogIO::POST ;
+        //std::vector< CountedPtr<Scantable> > tabs ;
+        //tabs.push_back( ssig ) ;
+        //tabs.push_back( sref ) ;
+        //out = merge( tabs ) ;
+        tmp[i] = ssig ;
+      }
+      else {
+        STSelector sel = STSelector() ;
+        vector<int> v( 1, sifno_unique[i] ) ; 
+        sel.setIFs( v ) ;
+        ssig->setSelection( sel ) ;
+        sel.reset() ;
+        v[0] = rifno_unique[i] ;
+        sel.setIFs( v ) ;
+        sref->setSelection( sel ) ;
+        sel.reset() ;
+        if ( antname.find( "APEX" ) != string::npos ) {
+          tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
+          //tmp[i] = dofold( ssig, sref, doffset ) ;
+        }
+        else {
+          tmp[i] = dofold( ssig, sref, doffset ) ;
+        }
+        ssig->unsetSelection() ;
+        sref->unsetSelection() ;
+      }
+    }
+  }
+
+  if ( tmp.size() > 1 ) {
+    out = merge( tmp ) ;
+  }
+  else {
+    out = tmp[0] ;
+  }
+
+  // flux unit
+  out->setFluxUnit( "K" ) ;
+
+  return out ;
+}
+
+CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
+{
+  CountedPtr<Scantable> out ;
+
+  return out ;
+}
+
+vector<float> STMath::getSpectrumFromTime( string reftime, 
+                                           CountedPtr<Scantable>& s,
+                                           string mode ) 
+{
+  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
+  vector<float> sp ;
+
+  if ( s->nrow() == 0 ) {
+    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
+    return sp ;
+  }
+  else if ( s->nrow() == 1 ) {
+    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
+    return s->getSpectrum( 0 ) ; 
+  }
+  else {
+    vector<int> idx = getRowIdFromTime( reftime, s ) ;
+    if ( mode == "before" ) {
+      int id = -1 ;
+      if ( idx[0] != -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[1] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;
+    }
+    else if ( mode == "after" ) {
+      int id = -1 ;
+      if ( idx[1] != -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[0] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;
+    }
+    else if ( mode == "nearest" ) {
+      int id = -1 ;
+      if ( idx[0] == -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[1] == -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        id = idx[0] ;
+      }
+      else {
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
+          id = idx[1] ;
+        }
+        else {
+          id = idx[0] ;
+        }
+      }
+      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+      sp = s->getSpectrum( id ) ;      
+    }
+    else if ( mode == "linear" ) {
+      if ( idx[0] == -1 ) {
+        // use after
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
+        int id = idx[1] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else if ( idx[1] == -1 ) {
+        // use before
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        // use before
+        //os << "No need to interporate." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
+        sp = s->getSpectrum( id ) ;
+      }
+      else {
+        // do interpolation
+        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
+        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
+        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
+          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
+          sp.push_back( v ) ;
+        }
+      }
+    }
+    else {
+      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
+    }
+    return sp ;
+  }
+}
+
+double STMath::getMJD( string strtime ) 
+{
+  if ( strtime.find("/") == string::npos ) {
+    // MJD time string
+    return atof( strtime.c_str() ) ;
+  }
+  else {
+    // string in YYYY/MM/DD/HH:MM:SS format
+    uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
+    uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
+    uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
+    uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
+    uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
+    uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
+    Time t( year, month, day, hour, minute, sec ) ;
+    return t.modifiedJulianDay() ;
+  }
+}
+
+vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
+{
+  double reft = getMJD( reftime ) ;
+  double dtmin = 1.0e100 ;
+  double dtmax = -1.0e100 ;
+  vector<double> dt ;
+  int just_before = -1 ;
+  int just_after = -1 ;
+  for ( int i = 0 ; i < s->nrow() ; i++ ) {
+    dt.push_back( getMJD( s->getTime( i ) ) - reft ) ; 
+  }
+  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
+    if ( dt[i] > 0.0 ) {
+      // after reftime
+      if ( dt[i] < dtmin ) {
+        just_after = i ;
+        dtmin = dt[i] ;
+      }
+    }
+    else if ( dt[i] < 0.0 ) {
+      // before reftime
+      if ( dt[i] > dtmax ) {
+        just_before = i ;
+        dtmax = dt[i] ;
+      }
+    }
+    else {
+      // just a reftime
+      just_before = i ;
+      just_after = i ;
+      dtmax = 0 ;
+      dtmin = 0 ;
+      break ;
+    }
+  }
+
+  vector<int> v ;
+  v.push_back( just_before ) ;
+  v.push_back( just_after ) ;
+
+  return v ;
+}
+
+vector<float> STMath::getTcalFromTime( string reftime, 
+                                       CountedPtr<Scantable>& s,
+                                       string mode ) 
+{
+  LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
+  vector<float> tcal ;
+  STTcal tcalTable = s->tcal() ;
+  String time ;
+  Vector<Float> tcalval ;
+  if ( s->nrow() == 0 ) {
+    os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
+    return tcal ;
+  }
+  else if ( s->nrow() == 1 ) {
+    uInt tcalid = s->getTcalId( 0 ) ;
+    //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+    tcalTable.getEntry( time, tcalval, tcalid ) ;
+    tcalval.tovector( tcal ) ;
+    return tcal ;
+  }
+  else {
+    vector<int> idx = getRowIdFromTime( reftime, s ) ;
+    if ( mode == "before" ) {
+      int id = -1 ;
+      if ( idx[0] != -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[1] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      uInt tcalid = s->getTcalId( id ) ;
+      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+      tcalTable.getEntry( time, tcalval, tcalid ) ;
+      tcalval.tovector( tcal ) ;
+    }
+    else if ( mode == "after" ) {
+      int id = -1 ;
+      if ( idx[1] != -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[0] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      uInt tcalid = s->getTcalId( id ) ;
+      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+      tcalTable.getEntry( time, tcalval, tcalid ) ;
+      tcalval.tovector( tcal ) ;
+    }
+    else if ( mode == "nearest" ) {
+      int id = -1 ;
+      if ( idx[0] == -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[1] == -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        id = idx[0] ;
+      }
+      else {
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
+          id = idx[1] ;
+        }
+        else {
+          id = idx[0] ;
+        }
+      }
+      uInt tcalid = s->getTcalId( id ) ;
+      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+      tcalTable.getEntry( time, tcalval, tcalid ) ;
+      tcalval.tovector( tcal ) ;
+    }
+    else if ( mode == "linear" ) {
+      if ( idx[0] == -1 ) {
+        // use after
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
+        int id = idx[1] ;
+        uInt tcalid = s->getTcalId( id ) ;
+        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+        tcalTable.getEntry( time, tcalval, tcalid ) ;
+        tcalval.tovector( tcal ) ;
+      }
+      else if ( idx[1] == -1 ) {
+        // use before
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
+        int id = idx[0] ;
+        uInt tcalid = s->getTcalId( id ) ;
+        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+        tcalTable.getEntry( time, tcalval, tcalid ) ;
+        tcalval.tovector( tcal ) ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        // use before
+        //os << "No need to interporate." << LogIO::POST ;
+        int id = idx[0] ;
+        uInt tcalid = s->getTcalId( id ) ;
+        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
+        tcalTable.getEntry( time, tcalval, tcalid ) ;
+        tcalval.tovector( tcal ) ;
+      }
+      else {
+        // do interpolation
+        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        vector<float> tcal0 ;
+        vector<float> tcal1 ;
+        uInt tcalid0 = s->getTcalId( idx[0] ) ;
+        uInt tcalid1 = s->getTcalId( idx[1] ) ;
+        tcalTable.getEntry( time, tcalval, tcalid0 ) ;
+        tcalval.tovector( tcal0 ) ;
+        tcalTable.getEntry( time, tcalval, tcalid1 ) ;
+        tcalval.tovector( tcal1 ) ;        
+        for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
+          float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
+          tcal.push_back( v ) ;
+        }
+      }
+    }
+    else {
+      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
+    }
+    return tcal ;
+  }
+}
+
+vector<float> STMath::getTsysFromTime( string reftime, 
+                                       CountedPtr<Scantable>& s,
+                                       string mode ) 
+{
+  LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
+  ArrayColumn<Float> tsysCol ;
+  tsysCol.attach( s->table(), "TSYS" ) ;
+  vector<float> tsys ;
+  String time ;
+  Vector<Float> tsysval ;
+  if ( s->nrow() == 0 ) {
+    os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
+    return tsys ;
+  }
+  else if ( s->nrow() == 1 ) {
+    //os << "use row " << 0 << LogIO::POST ;
+    tsysval = tsysCol( 0 ) ;
+    tsysval.tovector( tsys ) ;
+    return tsys ;
+  }
+  else {
+    vector<int> idx = getRowIdFromTime( reftime, s ) ;
+    if ( mode == "before" ) {
+      int id = -1 ;
+      if ( idx[0] != -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[1] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << LogIO::POST ;
+      tsysval = tsysCol( id ) ;
+      tsysval.tovector( tsys ) ;
+    }
+    else if ( mode == "after" ) {
+      int id = -1 ;
+      if ( idx[1] != -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[0] != -1 ) {
+        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
+        id = idx[1] ;
+      }
+      //os << "use row " << id << LogIO::POST ;
+      tsysval = tsysCol( id ) ;
+      tsysval.tovector( tsys ) ;
+    }
+    else if ( mode == "nearest" ) {
+      int id = -1 ;
+      if ( idx[0] == -1 ) {
+        id = idx[1] ;
+      }
+      else if ( idx[1] == -1 ) {
+        id = idx[0] ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        id = idx[0] ;
+      }
+      else {
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
+          id = idx[1] ;
+        }
+        else {
+          id = idx[0] ;
+        }
+      }
+      //os << "use row " << id << LogIO::POST ;
+      tsysval = tsysCol( id ) ;
+      tsysval.tovector( tsys ) ;
+    }
+    else if ( mode == "linear" ) {
+      if ( idx[0] == -1 ) {
+        // use after
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
+        int id = idx[1] ;
+        //os << "use row " << id << LogIO::POST ;
+        tsysval = tsysCol( id ) ;
+        tsysval.tovector( tsys ) ;
+      }
+      else if ( idx[1] == -1 ) {
+        // use before
+        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << LogIO::POST ;
+        tsysval = tsysCol( id ) ;
+        tsysval.tovector( tsys ) ;
+      }
+      else if ( idx[0] == idx[1] ) {
+        // use before
+        //os << "No need to interporate." << LogIO::POST ;
+        int id = idx[0] ;
+        //os << "use row " << id << LogIO::POST ;
+        tsysval = tsysCol( id ) ;
+        tsysval.tovector( tsys ) ;
+      }
+      else {
+        // do interpolation
+        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
+        double t0 = getMJD( s->getTime( idx[0] ) ) ;
+        double t1 = getMJD( s->getTime( idx[1] ) ) ;
+        double tref = getMJD( reftime ) ;
+        vector<float> tsys0 ;
+        vector<float> tsys1 ;
+        tsysval = tsysCol( idx[0] ) ;
+        tsysval.tovector( tsys0 ) ;
+        tsysval = tsysCol( idx[1] ) ;
+        tsysval.tovector( tsys1 ) ;        
+        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
+          float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
+          tsys.push_back( v ) ;
+        }
+      }
+    }
+    else {
+      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
+    }
+    return tsys ;
+  }
+}
+
+vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
+                                            CountedPtr<Scantable>& off,
+                                            CountedPtr<Scantable>& sky,
+                                            CountedPtr<Scantable>& hot,
+                                            CountedPtr<Scantable>& cold,
+                                            int index,
+                                            string antname ) 
+{
+  string reftime = on->getTime( index ) ;
+  vector<int> ii( 1, on->getIF( index ) ) ;
+  vector<int> ib( 1, on->getBeam( index ) ) ;
+  vector<int> ip( 1, on->getPol( index ) ) ;
+  vector<int> ic( 1, on->getScan( index ) ) ;
+  STSelector sel = STSelector() ;
+  sel.setIFs( ii ) ;
+  sel.setBeams( ib ) ;
+  sel.setPolarizations( ip ) ;
+  sky->setSelection( sel ) ;
+  hot->setSelection( sel ) ;
+  //cold->setSelection( sel ) ;
+  off->setSelection( sel ) ;
+  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
+  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
+  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
+  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
+  vector<float> spec = on->getSpectrum( index ) ;
+  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
+  vector<float> sp( tcal.size() ) ;
+  if ( antname.find( "APEX" ) != string::npos ) {
+    // using gain array
+    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
+      float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
+        * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
+      sp[j] = v ;
+    }
+  }
+  else {
+    // Chopper-Wheel calibration (Ulich & Haas 1976)
+    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
+      float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
+      sp[j] = v ;
+    }
+  }
+  sel.reset() ;
+  sky->unsetSelection() ;
+  hot->unsetSelection() ;
+  //cold->unsetSelection() ;
+  off->unsetSelection() ;
+
+  return sp ;
+}
+
+vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
+                                            CountedPtr<Scantable>& off,
+                                            int index ) 
+{
+  string reftime = on->getTime( index ) ;
+  vector<int> ii( 1, on->getIF( index ) ) ;
+  vector<int> ib( 1, on->getBeam( index ) ) ;
+  vector<int> ip( 1, on->getPol( index ) ) ;
+  vector<int> ic( 1, on->getScan( index ) ) ;
+  STSelector sel = STSelector() ;
+  sel.setIFs( ii ) ;
+  sel.setBeams( ib ) ;
+  sel.setPolarizations( ip ) ;
+  off->setSelection( sel ) ;
+  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
+  vector<float> spec = on->getSpectrum( index ) ;
+  //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
+  //vector<float> tsys = on->getTsysVec( index ) ;
+  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
+  Vector<Float> tsys = tsysCol( index ) ;
+  vector<float> sp( spec.size() ) ;
+  // ALMA Calibration
+  // 
+  // Ta* = Tsys * ( ON - OFF ) / OFF
+  //
+  // 2010/01/07 Takeshi Nakazato
+  unsigned int tsyssize = tsys.nelements() ;
+  unsigned int spsize = sp.size() ;
+  for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
+    float tscale = 0.0 ;
+    if ( tsyssize == spsize ) 
+      tscale = tsys[j] ;
+    else 
+      tscale = tsys[0] ;
+    float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
+    sp[j] = v ;
+  }
+  sel.reset() ;
+  off->unsetSelection() ;
+
+  return sp ;
+}
+
+vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
+                                              CountedPtr<Scantable>& ref,
+                                              CountedPtr<Scantable>& sky,
+                                              CountedPtr<Scantable>& hot,
+                                              CountedPtr<Scantable>& cold,
+                                              int index ) 
+{
+  string reftime = sig->getTime( index ) ;
+  vector<int> ii( 1, sig->getIF( index ) ) ;
+  vector<int> ib( 1, sig->getBeam( index ) ) ;
+  vector<int> ip( 1, sig->getPol( index ) ) ;
+  vector<int> ic( 1, sig->getScan( index ) ) ;
+  STSelector sel = STSelector() ;
+  sel.setIFs( ii ) ;
+  sel.setBeams( ib ) ;
+  sel.setPolarizations( ip ) ;
+  sky->setSelection( sel ) ;
+  hot->setSelection( sel ) ;
+  //cold->setSelection( sel ) ;
+  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
+  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
+  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
+  vector<float> spref = ref->getSpectrum( index ) ;
+  vector<float> spsig = sig->getSpectrum( index ) ;
+  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
+  vector<float> sp( tcal.size() ) ;
+  for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
+    float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
+    sp[j] = v ;
+  }
+  sel.reset() ;
+  sky->unsetSelection() ;
+  hot->unsetSelection() ;
+  //cold->unsetSelection() ;
+
+  return sp ;
+}
+
+vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
+                                              CountedPtr<Scantable>& ref,
+                                              vector< CountedPtr<Scantable> >& sky,
+                                              vector< CountedPtr<Scantable> >& hot,
+                                              vector< CountedPtr<Scantable> >& cold,
+                                              int index ) 
+{
+  string reftime = sig->getTime( index ) ;
+  vector<int> ii( 1, sig->getIF( index ) ) ;
+  vector<int> ib( 1, sig->getBeam( index ) ) ;
+  vector<int> ip( 1, sig->getPol( index ) ) ;
+  vector<int> ic( 1, sig->getScan( index ) ) ;
+  STSelector sel = STSelector() ;
+  sel.setIFs( ii ) ;
+  sel.setBeams( ib ) ;
+  sel.setPolarizations( ip ) ;
+  sky[0]->setSelection( sel ) ;
+  hot[0]->setSelection( sel ) ;
+  //cold[0]->setSelection( sel ) ;
+  vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
+  vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
+  //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
+  vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
+  sel.reset() ;
+  ii[0] = ref->getIF( index ) ;
+  sel.setIFs( ii ) ;
+  sel.setBeams( ib ) ;
+  sel.setPolarizations( ip ) ;
+  sky[1]->setSelection( sel ) ;
+  hot[1]->setSelection( sel ) ;
+  //cold[1]->setSelection( sel ) ;
+  vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
+  vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
+  //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
+  vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ;  
+  vector<float> spref = ref->getSpectrum( index ) ;
+  vector<float> spsig = sig->getSpectrum( index ) ;
+  vector<float> sp( tcals.size() ) ;
+  for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
+    float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
+    sp[j] = v ;
+  }
+  sel.reset() ;
+  sky[0]->unsetSelection() ;
+  hot[0]->unsetSelection() ;
+  //cold[0]->unsetSelection() ;
+  sky[1]->unsetSelection() ;
+  hot[1]->unsetSelection() ;
+  //cold[1]->unsetSelection() ;
+
+  return sp ;
+}

trunk/src/STMath.h

@@ -119 +119 @@
                   const std::string& mode, bool tsys=false );
 
+  // array operation
+  casa::CountedPtr<Scantable>
+    arrayOperate( const casa::CountedPtr<Scantable>& in,
+                  const std::vector<float> val,
+                  const std::string& mode,
+                  const std::string& opmode="channel",  
+                  bool tsys=false );
+
+  // channel operation
+  casa::CountedPtr<Scantable>
+    arrayOperateChannel( const casa::CountedPtr<Scantable>& in, 
+                         const std::vector<float> val,
+                         const std::string& mode, bool tsys=false );
+
+  // row operation
+  casa::CountedPtr<Scantable>
+    arrayOperateRow( const casa::CountedPtr<Scantable>& in, 
+                     const std::vector<float> val,
+                     const std::string& mode, bool tsys=false );
+
+  // 2d array operation
+  casa::CountedPtr<Scantable>
+    array2dOperate( const casa::CountedPtr<Scantable>& in,
+                  const std::vector< std::vector<float> > val,
+                  const std::string& mode, bool tsys=false );
+
   casa::CountedPtr<Scantable>
     binaryOperate( const casa::CountedPtr<Scantable>& left,
@@ -134 +160 @@
   /**
     * Calibrate total power scans (translated from GBTIDL)
-    * @param calon uncalibrated Scantable with CAL noise signal
+    * @param calon uncalibrated Scantable with CAL noise signal 
     * @param caloff uncalibrated Scantable with no CAL signal
     * @param tcal optional scalar Tcal, CAL temperature (K)
-    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
+    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 
     * (spectrum - average of the two CAL on and off spectra;
     * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
-    */
+    */              
   casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
                                             const casa::CountedPtr<Scantable>& caloff,
@@ -151 +177 @@
     * @param smoothref optional Boxcar smooth width of the reference scans
     * default: no smoothing (=1)
-    * @param tsysv optional scalar Tsys value at the zenith, required to
-    * set tau, as well
+    * @param tsysv optional scalar Tsys value at the zenith, required to 
+    * set tau, as well 
     * @param tau optional scalar Tau value
     * @return casa::CountedPtr<Scantable> which holds combined scans
@@ -163 +189 @@
                                         casa::Float tau=0.0 );
 
- /**
+  /**
     * Calibrate GBT Nod scan pairs (translated from GBTIDL)
     * @param s Scantable which contains Nod scans
@@ -170 +196 @@
     * @param tsysv optional scalar Tsys value at the zenith, required to
     * set tau, as well
-    * @param tau optional scalar Tau value
+    * @param tau optional scalar Tau value 
     * @param tcal optional scalar Tcal, CAL temperature (K)
     * @return casa::CountedPtr<Scantable> which holds calibrated scans
@@ -199 +225 @@
                                     casa::Float tcal=0.0 );
 
+  /**
+   * Calibrate data with Chopper-Wheel like calibration method 
+   * which adopts position switching by antenna motion, 
+   * wobbler (nutator) switching and On-The-Fly observation.
+   * 
+   * The method is applicable to APEX, and other telescopes other than GBT.
+   *
+   * @param a Scantable which contains ON and OFF scans
+   * @param a string that indicates calibration mode 
+   * @param a string that indicates antenna name
+   **/
+  casa::CountedPtr<Scantable> cwcal( const casa::CountedPtr<Scantable>& s,
+                                       const casa::String calmode, 
+                                       const casa::String antname );
+
+  /**
+   * Calibrate frequency switched scans with Chopper-Wheel like 
+   * calibration method.
+   *
+   * The method is applicable to APEX, and other telescopes other than GBT.
+   * 
+   * @param a Scantable which contains ON and OFF scans
+   * @param a string that indicates antenna name
+   **/
+  casa::CountedPtr<Scantable> cwcalfs( const casa::CountedPtr<Scantable>& s,
+                                       const casa::String antname );
+
+
+  /**
+   * Folding frequency-switch data
+   * @param sig
+   * @param ref
+   * @param choffset
+   **/
+  casa::CountedPtr<Scantable> dofold( const casa::CountedPtr<Scantable> &sig,
+                                      const casa::CountedPtr<Scantable> &ref,
+                                      casa::Double choffset,
+                                      casa::Double choffset = 0.0 );
+
+  /**
+   * ALMA calibration
+   **/
+  casa::CountedPtr<Scantable> almacal( const casa::CountedPtr<Scantable>& s, 
+                                       const casa::String calmode ) ;
+  casa::CountedPtr<Scantable> almacalfs( const casa::CountedPtr<Scantable>& s ) ;
 
   casa::CountedPtr<Scantable>
@@ -206 +277 @@
                                const std::vector<bool>& mask,
                                const std::string& which);
+
+  std::vector< int > minMaxChan(const casa::CountedPtr<Scantable>& in,
+                                const std::vector<bool>& mask,
+                                const std::string& which);
 
   casa::CountedPtr<Scantable> bin( const casa::CountedPtr<Scantable>& in,
@@ -266 +341 @@
              double end, const std::string& mode="frequency");
 
+  // test for average spectra with different channel/resolution
+  casa::CountedPtr<Scantable>
+    new_average( const std::vector<casa::CountedPtr<Scantable> >& in,
+                 const bool& compel, 
+                 const std::vector<bool>& mask = std::vector<bool>(),
+                 const std::string& weight = "NONE",
+                 const std::string& avmode = "SCAN" )
+    throw (casa::AipsError) ;
+
 private:
   casa::CountedPtr<Scantable>  applyToPol( const casa::CountedPtr<Scantable>& in,
@@ -301 +385 @@
     maskedArray( const casa::Vector<casa::Float>& s,
                  const casa::Vector<casa::uChar>& f );
+  casa::MaskedArray<casa::Double>
+    maskedArray( const casa::Vector<casa::Double>& s,
+                 const casa::Vector<casa::uChar>& f );
   casa::Vector<casa::uChar>
     flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
 
+  vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
+  vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
+  vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
+  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
+
+  // Chopper-Wheel type calibration
+  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
+                                      casa::CountedPtr<Scantable>& off,
+                                      casa::CountedPtr<Scantable>& sky,
+                                      casa::CountedPtr<Scantable>& hot,
+                                      casa::CountedPtr<Scantable>& cold,
+                                      int index,
+                                      string antname ) ;
+  // Tsys * (ON-OFF)/OFF
+  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
+                                      casa::CountedPtr<Scantable>& off,
+                                      int index ) ;
+  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
+                                        casa::CountedPtr<Scantable>& ref,
+                                        casa::CountedPtr<Scantable>& sky,
+                                        casa::CountedPtr<Scantable>& hot,
+                                        casa::CountedPtr<Scantable>& cold,
+                                        int index ) ;
+  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
+                                        casa::CountedPtr<Scantable>& ref,
+                                        vector< casa::CountedPtr<Scantable> >& sky,
+                                        vector< casa::CountedPtr<Scantable> >& hot,
+                                        vector< casa::CountedPtr<Scantable> >& cold,
+                                        int index ) ;
+  double getMJD( string strtime ) ;
+
   bool insitu_;
 };

trunk/src/STMathWrapper.h

@@ -73 +73 @@
   { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); }
 
+  ScantableWrapper arrayOperate( const ScantableWrapper& in, 
+                                 const std::vector<float> val,
+                                 const std::string& mode,
+                                 bool tsys=false )
+  { return ScantableWrapper(STMath::arrayOperateChannel(in.getCP(), val, mode, tsys)); }
+
+  ScantableWrapper array2dOperate( const ScantableWrapper& in, 
+                                   const std::vector< std::vector<float> > val,
+                                   const std::string& mode, bool tsys=false )
+  { return ScantableWrapper(STMath::array2dOperate(in.getCP(), val, mode, tsys)); }
+
   ScantableWrapper binaryOperate( const ScantableWrapper& left,
                                   const ScantableWrapper& right,
@@ -121 +132 @@
   { return STMath::statistic(in.getCP(), mask, which); }
 
+  std::vector<int> minMaxChan(const ScantableWrapper& in,
+                               const std::vector<bool>& mask,
+                               const std::string& which)
+  { return STMath::minMaxChan(in.getCP(), mask, which); }
+
   ScantableWrapper bin( const ScantableWrapper& in, int width=5)
   { return ScantableWrapper(STMath::bin(in.getCP(), width)); }
@@ -175 +191 @@
                                    const std::string& refTime,
                                    const std::string& method  )
-  { return ScantableWrapper(STMath::frequencyAlign(in.getCP(), refTime, method)); }
+  { return ScantableWrapper(STMath::frequencyAlign(in.getCP())); }
 
   ScantableWrapper convertPolarisation( const ScantableWrapper& in,
@@ -191 +207 @@
                                             mode)); }
 
+  // test for average spectra with different channel/resolution
+  ScantableWrapper
+    new_average( const std::vector<ScantableWrapper>& in,
+                 const bool& compel,
+                 const std::vector<bool>& mask,
+                 const std::string& weight,
+                 const std::string& avmode )
+  {
+    std::vector<casa::CountedPtr<Scantable> > sts;
+    for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP());
+    return ScantableWrapper(STMath::new_average(sts, compel, mask, weight, avmode));
+  }
+
+  // cwcal
+  ScantableWrapper cwcal( const ScantableWrapper &in,
+                          const std::string calmode,
+                          const std::string antname )
+  {
+    casa::CountedPtr<Scantable> tab = in.getCP() ;
+    casa::String mode( calmode ) ;
+    casa::String name( antname ) ;
+    return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ;
+  }
+  // almacal
+  ScantableWrapper almacal( const ScantableWrapper &in,
+                          const std::string calmode )
+  {
+    casa::CountedPtr<Scantable> tab = in.getCP() ;
+    casa::String mode( calmode ) ;
+    return ScantableWrapper( STMath::almacal( tab, mode ) ) ;
+  }
 };
 

trunk/src/STMolecules.cpp

@@ -14 +14 @@
 #include <tables/Tables/SetupNewTab.h>
 #include <tables/Tables/ScaColDesc.h>
+#include <tables/Tables/ArrColDesc.h>
 #include <tables/Tables/TableRecord.h>
 #include <tables/Tables/TableParse.h>
 #include <tables/Tables/TableRow.h>
 #include <casa/Containers/RecordField.h>
+
+#include <tables/Tables/TableProxy.h>
 
 #include "STMolecules.h"
@@ -59 +62 @@
 {
   // add to base class table
-  table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
-  table_.addColumn(ScalarColumnDesc<String>("NAME"));
-  table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
+  //table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
+  table_.addColumn(ArrayColumnDesc<Double>("RESTFREQUENCY"));
+  //table_.addColumn(ScalarColumnDesc<String>("NAME"));
+  table_.addColumn(ArrayColumnDesc<String>("NAME"));
+  //table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
+  table_.addColumn(ArrayColumnDesc<String>("FORMATTEDNAME"));
   table_.rwKeywordSet().define("UNIT", String("Hz"));
   // new cached columns
@@ -69 +75 @@
 }
 
+/***
 uInt STMolecules::addEntry( Double restfreq, const String& name,
                             const String& formattedname )
@@ -94 +101 @@
   return resultid;
 }
-
+***/
+uInt STMolecules::addEntry( Vector<Double> restfreq, const Vector<String>& name,
+                            const Vector<String>& formattedname )
+{
+// How to handle this...?
+  Table result =
+    table_( nelements(table_.col("RESTFREQUENCY")) == uInt (restfreq.size()) && 
+            all(table_.col("RESTFREQUENCY")== restfreq) );
+  uInt resultid = 0;
+  if ( result.nrow() > 0) {
+    ROScalarColumn<uInt> c(result, "ID");
+    c.get(0, resultid);
+  } else {
+    uInt rno = table_.nrow();
+    table_.addRow();
+    // get last assigned _id and increment
+    if ( rno > 0 ) {
+      idCol_.get(rno-1, resultid);
+      resultid++;
+    }
+    restfreqCol_.put(rno, restfreq);
+    nameCol_.put(rno, name);
+    formattednameCol_.put(rno, formattedname);
+    idCol_.put(rno, resultid);
+  }
+  return resultid;
+}
+
+
+
+
+/***
 void STMolecules::getEntry( Double& restfreq, String& name,
                             String& formattedname, uInt id ) const
@@ -104 +142 @@
   ROTableRow row(t);
   // get first row - there should only be one matching id
+
   const TableRecord& rec = row.get(0);
   restfreq = rec.asDouble("RESTFREQUENCY");
@@ -109 +148 @@
   formattedname = rec.asString("FORMATTEDNAME");
 }
+***/
+void STMolecules::getEntry( Vector<Double>& restfreq, Vector<String>& name,
+                            Vector<String>& formattedname, uInt id ) const
+{
+  Table t = table_(table_.col("ID") == Int(id) );
+  if (t.nrow() == 0 ) {
+    throw(AipsError("STMolecules::getEntry - id out of range"));
+  }
+  ROTableRow row(t);
+  // get first row - there should only be one matching id
+
+  const TableRecord& rec = row.get(0);
+  //restfreq = rec.asDouble("RESTFREQUENCY");
+  restfreq = rec.asArrayDouble("RESTFREQUENCY");
+  //name = rec.asString("NAME");
+  name = rec.asArrayString("NAME");
+  //formattedname = rec.asString("FORMATTEDNAME");
+  formattedname = rec.asArrayString("FORMATTEDNAME");
+}
 
 std::vector< double > asap::STMolecules::getRestFrequencies( ) const
 {
   std::vector<double> out;
+  //TableProxy itsTable(table_);
+  //Record rec;
   Vector<Double> rfs = restfreqCol_.getColumn();
   rfs.tovector(out);
+  //rec = itsTable.getVarColumn("RESTFREQUENCY", 0, 1, 1);
   return out;
 }
 
-double asap::STMolecules::getRestFrequency( uInt id ) const
-{
+std::vector< double > asap::STMolecules::getRestFrequency( uInt id ) const
+{
+  std::vector<double> out;
   Table t = table_(table_.col("ID") == Int(id) );
   if (t.nrow() == 0 ) {
@@ -126 +188 @@
   ROTableRow row(t);
   const TableRecord& rec = row.get(0);
-  return double(rec.asDouble("RESTFREQUENCY"));
-}
-
+  //return double(rec.asDouble("RESTFREQUENCY"));
+  Vector<Double> rfs = rec.asArrayDouble("RESTFREQUENCY");
+  rfs.tovector(out);
+  return out;
+}
+
+int asap::STMolecules::nrow() const
+{
+  return int(table_.nrow());
+}
 
 }//namespace

trunk/src/STMolecules.h

@@ -17 +17 @@
 #include <tables/Tables/Table.h>
 #include <tables/Tables/ScalarColumn.h>
+#include <tables/Tables/ArrayColumn.h>
+#include <casa/Arrays/Array.h>
 
 #include "STSubTable.h"
@@ -37 +39 @@
   STMolecules& operator=(const STMolecules& other);
 
+/***
   casa::uInt addEntry( casa::Double restfreq, const casa::String& name="",
                        const casa::String& formattedname="");
+***/
 
+  casa::uInt addEntry( casa::Vector<casa::Double> restfreq, const casa::Vector<casa::String>& name=casa::Vector<casa::String>(0),
+                       const casa::Vector<casa::String>& formattedname=casa::Vector<casa::String>(0));
+
+/***
   void getEntry( casa::Double& restfreq, casa::String& name,
                  casa::String& formattedname, casa::uInt id) const;
+***/
+  void getEntry( casa::Vector<casa::Double>& restfreq, casa::Vector<casa::String>& name,
+                 casa::Vector<casa::String>& formattedname, casa::uInt id) const;
 
   std::vector<double> getRestFrequencies() const;
-  double getRestFrequency( casa::uInt id ) const;
+  std::vector<double> getRestFrequency( casa::uInt id ) const;
   const casa::String& name() const { return name_; }
+  int nrow() const;
 
 private:
@@ -52 +64 @@
   //casa::Table table_;
   //casa::ScalarColumn<casa::uInt> freqidCol_;
-  casa::ScalarColumn<casa::Double> restfreqCol_;
-  casa::ScalarColumn<casa::String> nameCol_;
-  casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
+  //casa::ScalarColumn<casa::Double> restfreqCol_;
+  casa::ArrayColumn<casa::Double> restfreqCol_;
+  //casa::ScalarColumn<casa::String> nameCol_;
+  casa::ArrayColumn<casa::String> nameCol_;
+  //casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
+  casa::ArrayColumn<casa::String> formattednameCol_; // e.g. latex
 
 };

trunk/src/STSelector.cpp

@@ -87 +87 @@
 }
 
+void STSelector::setTypes( const std::vector< int >& types )
+{
+  setint("SRCTYPE", types);
+}
+
 void STSelector::setint(const std::string& key, const std::vector< int >& val)
 {
@@ -116 +121 @@
 }
 
+void STSelector::setRows( const std::vector< int >& rows )
+{
+  rowselection_ = rows;
+}
+
+// Table STSelector::apply( const Table& tab )
+// {
+//   if ( empty() ) {
+//     return sort(tab);
+//   }
+//   TableExprNode query;
+//   intidmap::const_iterator it;
+//   for (it = intselections_.begin(); it != intselections_.end(); ++it) {
+//     TableExprNode theset(Vector<Int>( (*it).second ));
+//     if ( query.isNull() ) {
+//       query = tab.col((*it).first).in(theset);
+//     } else {
+//       query = tab.col((*it).first).in(theset) && query;
+//     }
+//   }
+//   stringidmap::const_iterator it1;
+//   for (it1 = stringselections_.begin(); it1 != stringselections_.end(); ++it1) {
+//     TableExprNode theset(mathutil::toVectorString( (*it1).second ));
+//     if ( query.isNull() ) {
+//       query = tab.col((*it1).first).in(theset);
+//     } else {
+//       query = tab.col((*it1).first).in(theset) && query;
+//     }
+//   }
+//   // add taql query
+//   if ( taql_.size() > 0 ) {
+//     Table tmpt = tab;
+//     std::string pytaql = "USING STYLE PYTHON " + taql_;
+
+//     if ( !query.isNull() ) { // taql and selection
+//       tmpt = tableCommand(pytaql, tab(query));
+//     } else { // taql only
+//       tmpt = tableCommand(pytaql, tab);
+//     }
+//     return sort(tmpt);
+//   } else {
+//     if ( query.isNull() ) {
+//       return sort(tab);
+//     } else {
+//       return sort(tab(query));
+//     }
+//   }
+// }
 Table STSelector::apply( const Table& tab )
 {
   if ( empty() ) {
     return sort(tab);
+  }
+  Table basetab = tab;
+  // Important!! Be sure to apply row selection first. 
+  if (rowselection_.size() > 0){
+    //Vector<Int> intrownrs(rowselection_);
+    Vector<uInt> rownrs( rowselection_.size() );
+    convertArray(rownrs, Vector<Int> ( rowselection_ ));
+    basetab = tab( rownrs );
+    ///TableExprNode theset(Vector<Int>( rowselection_ ));
+    ///query = tab.nodeRownr().in(theset);
   }
   TableExprNode query;
@@ -126 +189 @@
     TableExprNode theset(Vector<Int>( (*it).second ));
     if ( query.isNull() ) {
-      query = tab.col((*it).first).in(theset);
+      //query = tab.col((*it).first).in(theset);
+      query = basetab.col((*it).first).in(theset);
     } else {
-      query = tab.col((*it).first).in(theset) && query;
+      //query = tab.col((*it).first).in(theset) && query;
+      query = basetab.col((*it).first).in(theset) && query;
     }
   }
@@ -135 +200 @@
     TableExprNode theset(mathutil::toVectorString( (*it1).second ));
     if ( query.isNull() ) {
-      query = tab.col((*it1).first).in(theset);
+      //query = tab.col((*it1).first).in(theset);
+      query = basetab.col((*it1).first).in(theset);
     } else {
-      query = tab.col((*it1).first).in(theset) && query;
+      //query = tab.col((*it1).first).in(theset) && query;
+      query = basetab.col((*it1).first).in(theset) && query;
     }
   }
   // add taql query
   if ( taql_.size() > 0 ) {
-    Table tmpt = tab;
+    //Table tmpt = tab;
+    Table tmpt = basetab;
     std::string pytaql = "USING STYLE PYTHON " + taql_;
 
     if ( !query.isNull() ) { // taql and selection
-      tmpt = tableCommand(pytaql, tab(query));
+      //tmpt = tableCommand(pytaql, tab(query));
+      tmpt = tableCommand(pytaql, basetab(query));
     } else { // taql only
-      tmpt = tableCommand(pytaql, tab);
+      //tmpt = tableCommand(pytaql, tab);
+      tmpt = tableCommand(pytaql, basetab);
     }
     return sort(tmpt);
   } else {
     if ( query.isNull() ) {
-      return sort(tab);
+      //return sort(tab);
+      return sort(basetab);
     } else {
-      return sort(tab(query));
+      //return sort(tab(query));
+      return sort(basetab(query));
     }
   }
@@ -191 +263 @@
 {
   return getint("CYCLENO");
+}
+
+std::vector< int > asap::STSelector::getTypes( ) const
+{
+  return getint("SRCTYPE") ;
 }
 
@@ -227 +304 @@
 bool asap::STSelector::empty( ) const
 {
-  return (intselections_.empty() && taql_.size() == 0 );
+  //return (intselections_.empty() && taql_.size() == 0 );
+  return (intselections_.empty() && taql_.size() == 0 && rowselection_.size() == 0);
 }
 

trunk/src/STSelector.h

@@ -45 +45 @@
   void setCycles(const std::vector<int>& cycs);
   void setName(const std::string&);
+  void setTypes(const std::vector<int>& types);
   virtual void setTaQL(const std::string& taql);
 
   void setSortOrder(const std::vector<std::string>& order);
+  void setRows(const std::vector<int>& rows);
 
   std::vector<int> getScans() const;
@@ -54 +56 @@
   std::vector<int> getPols() const;
   std::vector<int> getCycles() const;
+  std::vector<int> getTypes() const;
   std::vector<std::string> getPolTypes() const;
   std::string getTaQL() const { return taql_; }
@@ -86 +89 @@
   casa::Block<casa::String> order_;
   std::string taql_;
+  std::vector<int> rowselection_;
 };
 

trunk/src/STWriter.cpp

@@ -42 +42 @@
 #include <atnf/PKSIO/PKSMS2writer.h>
 #include <atnf/PKSIO/PKSSDwriter.h>
+#include <atnf/PKSIO/SrcType.h>
 
 #include <tables/Tables/Table.h>
@@ -154 +155 @@
     havexpol(ifs[i]) = nPol(ifs[i]) > 2;
   }
-  Vector<String> obstypes(2);
-  obstypes(0) = "TR";//on
-  obstypes(1) = "RF TR";//off
+//   Vector<String> obstypes(2);
+//   obstypes(0) = "TR";//on
+//   obstypes(1) = "RF TR";//off
   const Table table = inst->table();
 
@@ -182 +183 @@
     while (!beamit.pastEnd() ) {
       Table btable = beamit.table();
+      //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
+      //pksrec.direction = dirCol(0).getAngle("rad").getValue();
       TableIterator cycit(btable, "CYCLENO");
       ROArrayColumn<Double> srateCol(btable, "SCANRATE");
@@ -188 +191 @@
       Vector<Float> srateflt(sratedbl.nelements());
       convertArray(srateflt, sratedbl);
-      pksrec.scanRate = srateflt;
+      //pksrec.scanRate = srateflt;
+      pksrec.scanRate = sratedbl;
       ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
       spmCol.get(0, pksrec.srcPM);
@@ -200 +204 @@
       while (!cycit.pastEnd() ) {
         Table ctable = cycit.table();
-        TableIterator ifit(ctable, "IFNO");
+        TableIterator ifit(ctable, "IFNO", TableIterator::Ascending, TableIterator::HeapSort);
         MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
         pksrec.direction = dirCol(0).getAngle("rad").getValue();
@@ -213 +217 @@
           uInt nchan = specCol(0).nelements();
           Double crval,crpix;
+          //Vector<Double> restfreq;
           Float tmp0,tmp1,tmp2,tmp3,tmp4;
-          String stmp0,stmp1;
+          String tcalt;
+          Vector<String> stmp0, stmp1;
           inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
                                      rec.asuInt("FREQ_ID"));
@@ -239 +245 @@
           pksrec.fieldName = rec.asString("FIELDNAME");
           pksrec.srcName   = rec.asString("SRCNAME");
-          pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
+          //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
+          pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
           pksrec.bandwidth = nchan * abs(pksrec.freqInc);
           pksrec.azimuth   = rec.asFloat("AZIMUTH");
@@ -253 +260 @@
           pksrec.baseSub   = 0.0f;
           pksrec.xCalFctr  = 0.0;
+          pksrec.flagrow = rec.asuInt("FLAGROW");
 
           status = writer_->write(pksrec);
@@ -348 +356 @@
 }
 
-}
+// get obsType string from SRCTYPE value
+String STWriter::getObsTypes( Int srctype )
+{
+  String obsType ;
+  switch( srctype ) {
+  case Int(SrcType::PSON):
+    obsType = "PSON" ;
+    break ;
+  case Int(SrcType::PSOFF):
+    obsType = "PSOFF" ;
+    break ;
+  case Int(SrcType::NOD):
+    obsType = "NOD" ;
+    break ;
+  case Int(SrcType::FSON):
+    obsType = "FSON" ;
+    break ;
+  case Int(SrcType::FSOFF):
+    obsType = "FSOFF" ;
+    break ;
+  case Int(SrcType::SKY):
+    obsType = "SKY" ;
+    break ;
+  case Int(SrcType::HOT):
+    obsType = "HOT" ;
+    break ;
+  case Int(SrcType::WARM):
+    obsType = "WARM" ;
+    break ;
+  case Int(SrcType::COLD):
+    obsType = "COLD" ;
+    break ;
+  case Int(SrcType::PONCAL):
+    obsType = "PSON:CALON" ;
+    break ;
+  case Int(SrcType::POFFCAL):
+    obsType = "PSOFF:CALON" ;
+    break ;
+  case Int(SrcType::NODCAL):
+    obsType = "NOD:CALON" ;
+    break ;
+  case Int(SrcType::FONCAL):
+    obsType = "FSON:CALON" ;
+    break ;
+  case Int(SrcType::FOFFCAL):
+    obsType = "FSOFF:CALOFF" ;
+    break ;
+  case Int(SrcType::FSLO):
+    obsType = "FSLO" ;
+    break ;
+  case Int(SrcType::FLOOFF):
+    obsType = "FS:LOWER:OFF" ;
+    break ;
+  case Int(SrcType::FLOSKY):
+    obsType = "FS:LOWER:SKY" ;
+    break ;
+  case Int(SrcType::FLOHOT):
+    obsType = "FS:LOWER:HOT" ;
+    break ;
+  case Int(SrcType::FLOWARM):
+    obsType = "FS:LOWER:WARM" ;
+    break ;
+  case Int(SrcType::FLOCOLD):
+    obsType = "FS:LOWER:COLD" ;
+    break ;
+  case Int(SrcType::FSHI):
+    obsType = "FSHI" ;
+    break ;
+  case Int(SrcType::FHIOFF):
+    obsType = "FS:HIGHER:OFF" ;
+    break ;
+  case Int(SrcType::FHISKY):
+    obsType = "FS:HIGHER:SKY" ;
+    break ;
+  case Int(SrcType::FHIHOT):
+    obsType = "FS:HIGHER:HOT" ;
+    break ;
+  case Int(SrcType::FHIWARM):
+    obsType = "FS:HIGHER:WARM" ;
+    break ;
+  case Int(SrcType::FHICOLD):
+    obsType = "FS:HIGHER:COLD" ;
+    break ;
+  default:
+    obsType = "NOTYPE" ;
+  }
+
+  return obsType ;
+}
+
+}

trunk/src/STWriter.h

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

trunk/src/Scantable.cpp

@@ -11 +11 @@
 //
 #include <map>
+#include <fstream>
 
 #include <casa/aips.h>
@@ -24 +25 @@
 #include <casa/Arrays/Vector.h>
 #include <casa/Arrays/VectorSTLIterator.h>
+#include <casa/Arrays/Slice.h>
 #include <casa/BasicMath/Math.h>
 #include <casa/BasicSL/Constants.h>
@@ -29 +31 @@
 #include <casa/Containers/RecordField.h>
 #include <casa/Utilities/GenSort.h>
+#include <casa/Logging/LogIO.h>
 
 #include <tables/Tables/TableParse.h>
@@ -106 +109 @@
 {
   initFactories();
+
   Table tab(name, Table::Update);
   uInt version = tab.keywordSet().asuInt("VERSION");
@@ -121 +125 @@
   attach();
 }
+/*
+Scantable::Scantable(const std::string& name, Table::TableType ttype) :
+  type_(ttype)
+{
+  initFactories();
+  Table tab(name, Table::Update);
+  uInt version = tab.keywordSet().asuInt("VERSION");
+  if (version != version_) {
+    throw(AipsError("Unsupported version of ASAP file."));
+  }
+  if ( type_ == Table::Memory ) {
+    table_ = tab.copyToMemoryTable(generateName());
+  } else {
+    table_ = tab;
+  }
+
+  attachSubtables();
+  originalTable_ = table_;
+  attach();
+}
+*/
 
 Scantable::Scantable( const Scantable& other, bool clear )
@@ -199 +224 @@
   td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
   td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
-  td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
+
+  ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
+  refbeamnoColumn.setDefault(Int(-1));
+  td.addColumn(refbeamnoColumn);
+
+  ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
+  flagrowColumn.setDefault(uInt(0));
+  td.addColumn(flagrowColumn);
 
   td.addColumn(ScalarColumnDesc<Double>("TIME"));
@@ -257 +289 @@
   originalTable_ = table_;
 }
-
 
 void Scantable::attach()
@@ -285 +316 @@
   mfocusidCol_.attach(table_, "FOCUS_ID");
   mmolidCol_.attach(table_, "MOLECULE_ID");
+
+  //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
+  attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
+
+}
+
+template<class T, class T2>
+void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
+                                   const String& colName,
+                                   const T2& defValue)
+{
+  try {
+    col.attach(table_, colName);
+  } catch (TableError& err) {
+    String errMesg = err.getMesg();
+    if (errMesg == "Table column " + colName + " is unknown") {
+      table_.addColumn(ScalarColumnDesc<T>(colName));
+      col.attach(table_, colName);
+      col.fillColumn(static_cast<T>(defValue));
+    } else {
+      throw;
+    }
+  } catch (...) {
+    throw;
+  }
+}
+
+template<class T, class T2>
+void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
+                                   const String& colName,
+                                   const Array<T2>& defValue)
+{
+  try {
+    col.attach(table_, colName);
+  } catch (TableError& err) {
+    String errMesg = err.getMesg();
+    if (errMesg == "Table column " + colName + " is unknown") {
+      table_.addColumn(ArrayColumnDesc<T>(colName));
+      col.attach(table_, colName);
+
+      int size = 0;
+      ArrayIterator<T2>& it = defValue.begin();
+      while (it != defValue.end()) {
+        ++size;
+        ++it;
+      }
+      IPosition ip(1, size);
+      Array<T>& arr(ip);
+      for (int i = 0; i < size; ++i)
+        arr[i] = static_cast<T>(defValue[i]);
+
+      col.fillColumn(arr);
+    } else {
+      throw;
+    }
+  } catch (...) {
+    throw;
+  }
 }
 
@@ -627 +716 @@
 }
 
+void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
+{
+  for (uInt i=0; i<table_.nrow(); ++i) {
+    Vector<uChar> flgs = flagsCol_(i);
+    srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
+    flagsCol_.put(i, flgs);
+  }
+}
+
+std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
+{
+  Vector<uChar> flags;
+  flagsCol_.get(uInt(whichrow), flags);
+  srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
+  Vector<Bool> bflag(flags.shape());
+  convertArray(bflag, flags);
+  //bflag = !bflag;
+
+  std::vector<bool> mask;
+  bflag.tovector(mask);
+  return mask;
+}
+
+void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
+                                   Vector<uChar> flgs)
+{
+    Vector<Float> spcs = specCol_(whichrow);
+    uInt nchannel = nchan();
+    if (spcs.nelements() != nchannel) {
+      throw(AipsError("Data has incorrect number of channels"));
+    }
+    uChar userflag = 1 << 7;
+    if (unflag) {
+      userflag = 0 << 7;
+    }
+    if (clipoutside) {
+      for (uInt j = 0; j < nchannel; ++j) {
+        Float spc = spcs(j);
+        if ((spc >= uthres) || (spc <= dthres)) {
+          flgs(j) = userflag;
+        }
+      }
+    } else {
+      for (uInt j = 0; j < nchannel; ++j) {
+        Float spc = spcs(j);
+        if ((spc < uthres) && (spc > dthres)) {
+          flgs(j) = userflag;
+        }
+      }
+    }
+}
+
 void Scantable::flag(const std::vector<bool>& msk, bool unflag)
 {
@@ -673 +814 @@
     flagsCol_.put(i, flgs);
   }
+}
+
+void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
+{
+  if ( selector_.empty() && (rows.size() == table_.nrow()) )
+    throw(AipsError("Trying to flag whole scantable."));
+
+  uInt rowflag = (unflag ? 0 : 1);
+  std::vector<uInt>::const_iterator it;
+  for (it = rows.begin(); it != rows.end(); ++it)
+    flagrowCol_.put(*it, rowflag);
 }
 
@@ -793 +945 @@
   table_.keywordSet().get("FluxUnit", tmp);
   oss << setw(15) << "Flux Unit:" << tmp << endl;
-  Vector<Double> vec(moleculeTable_.getRestFrequencies());
+  //Vector<Double> vec(moleculeTable_.getRestFrequencies());
+  int nid = moleculeTable_.nrow();
+  Bool firstline = True;
   oss << setw(15) << "Rest Freqs:";
-  if (vec.nelements() > 0) {
-      oss << setprecision(10) << vec << " [Hz]" << endl;
-  } else {
-      oss << "none" << endl;
+  for (int i=0; i<nid; i++) {
+      Table t = table_(table_.col("MOLECULE_ID") == i);
+      if (t.nrow() >  0) {
+          Vector<Double> vec(moleculeTable_.getRestFrequency(i));
+          if (vec.nelements() > 0) {
+               if (firstline) {
+                   oss << setprecision(10) << vec << " [Hz]" << endl;
+                   firstline=False;
+               }
+               else{
+                   oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
+               }
+          } else {
+              oss << "none" << endl;
+          }
+      }
   }
 
@@ -901 +1067 @@
   const MDirection& md = getDirection(whichrow);
   const MEpoch& me = timeCol_(whichrow);
-  Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
+  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
   return freqTable_.getSpectralCoordinate(md, mp, me, rf,
                                           mfreqidCol_(whichrow));
@@ -975 +1142 @@
   const MDirection& md = getDirection(whichrow);
   const MEpoch& me = timeCol_(whichrow);
-  const Double& rf = mmolidCol_(whichrow);
+  //const Double& rf = mmolidCol_(whichrow);
+  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
   SpectralCoordinate spc =
     freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
@@ -992 +1160 @@
 }
 
-void Scantable::setRestFrequencies( double rf, const std::string& name,
+/**
+void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
                                           const std::string& unit )
+**/
+void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
+                                          const std::string& unit )
+
 {
   ///@todo lookup in line table to fill in name and formattedname
   Unit u(unit);
-  Quantum<Double> urf(rf, u);
-  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
+  //Quantum<Double> urf(rf, u);
+  Quantum<Vector<Double> >urf(rf, u);
+  Vector<String> formattedname(0);
+  //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
+
+  //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
+  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
   TableVector<uInt> tabvec(table_, "MOLECULE_ID");
   tabvec = id;
 }
 
-void Scantable::setRestFrequencies( const std::string& name )
+/**
+void asap::Scantable::setRestFrequencies( const std::string& name )
 {
   throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
+  ///@todo implement
+}
+**/
+void Scantable::setRestFrequencies( const vector<std::string>& name )
+{
+  throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
   ///@todo implement
 }
@@ -1044 +1229 @@
 void Scantable::addFit( const STFitEntry& fit, int row )
 {
-  cout << mfitidCol_(uInt(row)) << endl;
+  //cout << mfitidCol_(uInt(row)) << endl;
+  LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
+  os << mfitidCol_(uInt(row)) << LogIO::POST ;
   uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
   mfitidCol_.put(uInt(row), id);
@@ -1059 +1246 @@
 }
 
-std::string Scantable::getAntennaName() const
+String Scantable::getAntennaName() const
 {
   String out;
@@ -1078 +1265 @@
       Table subt = t( t.col("SCAN") == scanlist[i]+1 );
       if (subt.nrow()==0) {
-        cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
+        //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
+        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+        os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
         ret = 1;
         break;
@@ -1090 +1279 @@
           Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
           if ( subt2.nrow() == 0) {
-            cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
+            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+
+            //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
+            os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
             ret = 1;
             break;
@@ -1100 +1292 @@
           if (scan1seqn == 1 && scan2seqn == 2) {
             if (laston1 == laston2) {
-              cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+              //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
               i +=1;
             }
             else {
-              cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+              //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
             }
           }
           else if (scan1seqn==2 && scan2seqn == 1) {
             if (laston1 == laston2) {
-              cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
+              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+              //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
+              os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
               ret = 1;
               break;
@@ -1115 +1313 @@
           }
           else {
-            cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
+            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+            //cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
+            os<<LogIO::WARN<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
             ret = 1;
             break;
@@ -1122 +1322 @@
       }
       else {
-        cerr<<"The scan does not appear to be standard obsevation."<<endl;
+        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+        //cerr<<"The scan does not appear to be standard obsevation."<<endl;
+        os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
       }
     //if ( i >= nscan ) break;
@@ -1128 +1330 @@
   }
   else {
-    cerr<<"No reference to GBT_GO table."<<endl;
+    LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
+    //cerr<<"No reference to GBT_GO table."<<endl;
+    os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
     ret = 1;
   }
@@ -1142 +1346 @@
 }
 
+void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
+  throw( casa::AipsError )
+{
+  // assumed that all rows have same nChan
+  Vector<Float> arr = specCol_( 0 ) ;
+  int nChan = arr.nelements() ;
+
+  // if nmin < 0 or nmax < 0, nothing to do
+  if (  nmin < 0 ) {
+    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
+    }
+  if (  nmax < 0  ) {
+    throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
+  }
+
+  // if nmin > nmax, exchange values
+  if ( nmin > nmax ) {
+    int tmp = nmax ;
+    nmax = nmin ;
+    nmin = tmp ;
+    LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
+    os << "Swap values. Applied range is ["
+       << nmin << ", " << nmax << "]" << LogIO::POST ;
+  }
+
+  // if nmin exceeds nChan, nothing to do
+  if ( nmin >= nChan ) {
+    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
+  }
+
+  // if nmax exceeds nChan, reset nmax to nChan
+  if ( nmax >= nChan ) {
+    if ( nmin == 0 ) {
+      // nothing to do
+      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
+      os << "Whole range is selected. Nothing to do." << LogIO::POST ;
+      return ;
+    }
+    else {
+      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
+      os << "Specified maximum exceeds nChan. Applied range is ["
+         << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
+      nmax = nChan - 1 ;
+    }
+  }
+
+  // reshape specCol_ and flagCol_
+  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
+    reshapeSpectrum( nmin, nmax, irow ) ;
+  }
+
+  // update FREQUENCIES subtable
+  Double refpix ;
+  Double refval ;
+  Double increment ;
+  int freqnrow = freqTable_.table().nrow() ;
+  Vector<uInt> oldId( freqnrow ) ;
+  Vector<uInt> newId( freqnrow ) ;
+  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
+    freqTable_.getEntry( refpix, refval, increment, irow ) ;
+    /***
+     * need to shift refpix to nmin
+     * note that channel nmin in old index will be channel 0 in new one
+     ***/
+    refval = refval - ( refpix - nmin ) * increment ;
+    refpix = 0 ;
+    freqTable_.setEntry( refpix, refval, increment, irow ) ;
+  }
+
+  // update nchan
+  int newsize = nmax - nmin + 1 ;
+  table_.rwKeywordSet().define( "nChan", newsize ) ;
+
+  // update bandwidth
+  // assumed all spectra in the scantable have same bandwidth
+  table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
+
+  return ;
+}
+
+void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
+{
+  // reshape specCol_ and flagCol_
+  Vector<Float> oldspec = specCol_( irow ) ;
+  Vector<uChar> oldflag = flagsCol_( irow ) ;
+  uInt newsize = nmax - nmin + 1 ;
+  specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
+  flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
+
+  return ;
+}
+
+void asap::Scantable::regridChannel( int nChan, double dnu )
+{
+  LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
+  os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
+  // assumed that all rows have same nChan
+  Vector<Float> arr = specCol_( 0 ) ;
+  int oldsize = arr.nelements() ;
+
+  // if oldsize == nChan, nothing to do
+  if ( oldsize == nChan ) {
+    os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
+    return ;
+  }
+
+  // if oldChan < nChan, unphysical operation
+  if ( oldsize < nChan ) {
+    os << "Unphysical operation. Nothing to do." << LogIO::POST ;
+    return ;
+  }
+
+  // change channel number for specCol_ and flagCol_
+  Vector<Float> newspec( nChan, 0 ) ;
+  Vector<uChar> newflag( nChan, false ) ;
+  vector<string> coordinfo = getCoordInfo() ;
+  string oldinfo = coordinfo[0] ;
+  coordinfo[0] = "Hz" ;
+  setCoordInfo( coordinfo ) ;
+  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
+    regridChannel( nChan, dnu, irow ) ;
+  }
+  coordinfo[0] = oldinfo ;
+  setCoordInfo( coordinfo ) ;
+
+
+  // NOTE: this method does not update metadata such as
+  //       FREQUENCIES subtable, nChan, Bandwidth, etc.
+
+  return ;
+}
+
+void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
+{
+  // logging
+  //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
+  //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
+
+  Vector<Float> oldspec = specCol_( irow ) ;
+  Vector<uChar> oldflag = flagsCol_( irow ) ;
+  Vector<Float> newspec( nChan, 0 ) ;
+  Vector<uChar> newflag( nChan, false ) ;
+
+  // regrid
+  vector<double> abcissa = getAbcissa( irow ) ;
+  int oldsize = abcissa.size() ;
+  double olddnu = abcissa[1] - abcissa[0] ;
+  //int refChan = 0 ;
+  //double frac = 0.0 ;
+  //double wedge = 0.0 ;
+  //double pile = 0.0 ;
+  int ichan = 0 ;
+  double wsum = 0.0 ;
+  Vector<Float> z( nChan ) ;
+  z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
+  for ( int ii = 1 ; ii < nChan ; ii++ )
+    z[ii] = z[ii-1] + dnu ;
+  Vector<Float> zi( nChan+1 ) ;
+  Vector<Float> yi( oldsize + 1 ) ;
+  zi[0] = z[0] - 0.5 * dnu ;
+  zi[1] = z[0] + 0.5 * dnu ;
+  for ( int ii = 2 ; ii < nChan ; ii++ )
+    zi[ii] = zi[ii-1] + dnu ;
+  zi[nChan] = z[nChan-1] + 0.5 * dnu ;
+  yi[0] = abcissa[0] - 0.5 * olddnu ;
+  yi[1] = abcissa[1] + 0.5 * olddnu ;
+  for ( int ii = 2 ; ii < oldsize ; ii++ )
+    yi[ii] = abcissa[ii-1] + olddnu ;
+  yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
+  if ( dnu > 0.0 ) {
+    for ( int ii = 0 ; ii < nChan ; ii++ ) {
+      double zl = zi[ii] ;
+      double zr = zi[ii+1] ;
+      for ( int j = ichan ; j < oldsize ; j++ ) {
+        double yl = yi[j] ;
+        double yr = yi[j+1] ;
+        if ( yl <= zl ) {
+          if ( yr <= zl ) {
+            continue ;
+          }
+          else if ( yr <= zr ) {
+            newspec[ii] += oldspec[j] * ( yr - zl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += ( yr - zl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * dnu ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += dnu ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else if ( yl < zr ) {
+          if ( yr <= zr ) {
+              newspec[ii] += oldspec[j] * ( yr - yl ) ;
+              newflag[ii] = newflag[ii] || oldflag[j] ;
+              wsum += ( yr - yl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * ( zr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += ( zr - yl ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else {
+          ichan = j - 1 ;
+          break ;
+        }
+      }
+      newspec[ii] /= wsum ;
+      wsum = 0.0 ;
+    }
+  }
+  else if ( dnu < 0.0 ) {
+    for ( int ii = 0 ; ii < nChan ; ii++ ) {
+      double zl = zi[ii] ;
+      double zr = zi[ii+1] ;
+      for ( int j = ichan ; j < oldsize ; j++ ) {
+        double yl = yi[j] ;
+        double yr = yi[j+1] ;
+        if ( yl >= zl ) {
+          if ( yr >= zl ) {
+            continue ;
+          }
+          else if ( yr >= zr ) {
+            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( yr - zl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * abs( dnu ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( dnu ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else if ( yl > zr ) {
+          if ( yr >= zr ) {
+            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( yr - yl ) ;
+          }
+          else {
+            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
+            newflag[ii] = newflag[ii] || oldflag[j] ;
+            wsum += abs( zr - yl ) ;
+            ichan = j ;
+            break ;
+          }
+        }
+        else {
+          ichan = j - 1 ;
+          break ;
+        }
+      }
+      newspec[ii] /= wsum ;
+      wsum = 0.0 ;
+    }
+  }
+//    * ichan = 0
+//    ***/
+//   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
+//   pile += dnu ;
+//   wedge = olddnu * ( refChan + 1 ) ;
+//   while ( wedge < pile ) {
+//     newspec[0] += olddnu * oldspec[refChan] ;
+//     newflag[0] = newflag[0] || oldflag[refChan] ;
+//     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
+//     refChan++ ;
+//     wedge += olddnu ;
+//     wsum += olddnu ;
+//     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
+//   }
+//   frac = ( wedge - pile ) / olddnu ;
+//   wsum += ( 1.0 - frac ) * olddnu ;
+//   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
+//   newflag[0] = newflag[0] || oldflag[refChan] ;
+//   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
+//   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
+//   newspec[0] /= wsum ;
+//   //ofs << "newspec[0] = " << newspec[0] << endl ;
+//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+
+//   /***
+//    * ichan = 1 - nChan-2
+//    ***/
+//   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
+//     pile += dnu ;
+//     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
+//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
+//     refChan++ ;
+//     wedge += olddnu ;
+//     wsum = frac * olddnu ;
+//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+//     while ( wedge < pile ) {
+//       newspec[ichan] += olddnu * oldspec[refChan] ;
+//       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+//       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
+//       refChan++ ;
+//       wedge += olddnu ;
+//       wsum += olddnu ;
+//       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+//     }
+//     frac = ( wedge - pile ) / olddnu ;
+//     wsum += ( 1.0 - frac ) * olddnu ;
+//     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
+//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
+//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
+//     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
+//     newspec[ichan] /= wsum ;
+//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
+//   }
+
+//   /***
+//    * ichan = nChan-1
+//    ***/
+//   // NOTE: Assumed that all spectra have the same bandwidth
+//   pile += dnu ;
+//   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
+//   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
+//   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
+//   refChan++ ;
+//   wedge += olddnu ;
+//   wsum = frac * olddnu ;
+//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+//   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
+//     newspec[nChan-1] += olddnu * oldspec[jchan] ;
+//     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
+//     wsum += olddnu ;
+//     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
+//     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+//   }
+//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
+//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
+//   newspec[nChan-1] /= wsum ;
+//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
+
+//   specCol_.put( irow, newspec ) ;
+//   flagsCol_.put( irow, newflag ) ;
+
+//   // ofs.close() ;
+
+
+  return ;
+}
+
 std::vector<float> Scantable::getWeather(int whichrow) const
 {
@@ -1154 +1710 @@
 
 }
- //namespace asap
+//namespace asap

trunk/src/Scantable.h

@@ -29 +29 @@
 
 #include <coordinates/Coordinates/SpectralCoordinate.h>
+
+#include <casa/Arrays/Vector.h>
+#include <casa/Quanta/Quantum.h>
+
+#include <casa/Exceptions/Error.h>
 
 #include "Logger.h"
@@ -226 +231 @@
 
   /**
+   * Flag the data in a row-based manner. (CAS-1433 Wataru Kawasaki)
+   * param[in] rows    list of row numbers to be flagged
+   */
+  void flagRow( const std::vector<casa::uInt>& rows = std::vector<casa::uInt>(), bool unflag=false);
+
+  /**
+   * Get flagRow info at the specified row. If true, the whole data
+   * at the row should be flagged.
+   */
+  bool getFlagRow(int whichrow) const
+    { return (flagrowCol_(whichrow) > 0); }
+
+  /**
+   * Flag the data outside a specified range (in a channel-based manner).
+   * (CAS-1807 Wataru Kawasaki)
+   */
+  void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
+
+  /**
+   * Return a list of booleans with the size of nchan for a specified row, to get info
+   * about which channel is clipped.
+   */
+  std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
+  void srchChannelsToClip(casa::uInt whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag,
+                          casa::Vector<casa::uChar> flgs);
+
+  /**
    * Return a list of row numbers with respect to the original table.
    * @return a list of unsigned ints
@@ -276 +308 @@
   std::vector<uint> getScanNos() const { return getNumbers(scanCol_); }
   int getScan(int whichrow) const { return scanCol_(whichrow); }
+
+  //TT addition
+  std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); }
 
   /**
@@ -300 +335 @@
     { return azCol_(whichrow); }
   float getParAngle(int whichrow) const
-  { return focus().getParAngle(mfocusidCol_(whichrow)); }
+    { return focus().getParAngle(mfocusidCol_(whichrow)); }
+  int getTcalId(int whichrow) const
+    { return mtcalidCol_(whichrow); }
 
   std::string getSourceName(int whichrow) const
@@ -353 +390 @@
   std::vector<double> getRestFrequencies() const
     { return moleculeTable_.getRestFrequencies(); }
-
+  std::vector<double> getRestFrequency(int id) const
+    { return moleculeTable_.getRestFrequency(id); }
+
+  /**
   void setRestFrequencies(double rf, const std::string& name = "",
                           const std::string& = "Hz");
-  void setRestFrequencies(const std::string& name);
+  **/
+  // Modified by Takeshi Nakazato 05/09/2008
+  /***
+  void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "",
+                          const std::string& = "Hz");
+  ***/
+  void setRestFrequencies(vector<double> rf,
+                          const vector<std::string>& name = vector<std::string>(1,""),
+                          const std::string& = "Hz");
+
+  //void setRestFrequencies(const std::string& name);
+  void setRestFrequencies(const vector<std::string>& name);
 
   void shift(int npix);
@@ -390 +441 @@
    * @return antenna name string
    */
-  std::string getAntennaName() const;
+  casa::String getAntennaName() const;
 
   /**
@@ -414 +465 @@
   void parallactify(bool flag)
     { focus().setParallactify(flag); }
+
+  /**
+   * Reshape spectrum
+   * @param[in] nmin, nmax minimum and maximum channel
+   * @param[in] irow       row number
+   *
+   * 30/07/2008 Takeshi Nakazato
+   **/
+  void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError );
+  void reshapeSpectrum( int nmin, int nmax, int irow ) ;
+
+  /**
+   * Change channel number under fixed bandwidth
+   * @param[in] nchan, dnu new channel number and spectral resolution
+   * @param[in] irow       row number
+   *
+   * 27/08/2008 Takeshi Nakazato
+   **/
+  void regridChannel( int nchan, double dnu ) ;
+  void regridChannel( int nchan, double dnu, int irow ) ;
+
 
 private:
@@ -489 +561 @@
   casa::ScalarColumn<casa::Float> elCol_;
   casa::ScalarColumn<casa::String> srcnCol_, fldnCol_;
-  casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_;
+  casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_;
   casa::ScalarColumn<casa::Int> rbeamCol_, srctCol_;
   casa::ArrayColumn<casa::Float> specCol_, tsysCol_;
@@ -510 +582 @@
   void initFactories();
 
+  /**
+   * Add an auxiliary column to the main table and attach it to a
+   * cached column. Use for adding new columns that the original asap2
+   * tables do not have.
+   * @param[in] col      reference to the cached column to be attached
+   * @param[in] colName  column name in asap table
+   * @param[in] defValue default value to fill in the column
+   *
+   * 25/10/2009 Wataru Kawasaki
+   */
+  template<class T, class T2> void attachAuxColumnDef(casa::ScalarColumn<T>&,
+                                                       const casa::String&,
+                                                       const T2&);
+  template<class T, class T2> void attachAuxColumnDef(casa::ArrayColumn<T>&,
+                                                      const casa::String&,
+                                                      const casa::Array<T2>&);
 };
 

trunk/src/ScantableWrapper.h

@@ -109 +109 @@
     { table_->flag(msk, unflag); }
 
+  void flagRow(const std::vector<casa::uInt>& rows=std::vector<casa::uInt>(), bool unflag=false)
+    { table_->flagRow(rows, unflag); }
+
+  bool getFlagRow(int whichrow=0) const
+    { return table_->getFlagRow(whichrow); }
+
+  void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside=true, bool unflag=false)
+    { table_->clip(uthres, dthres, clipoutside, unflag); }
+
+  std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag) const
+    { return table_->getClipMask(whichrow, uthres, dthres, clipoutside, unflag); }
+
   std::string getSourceName(int whichrow=0) const
     { return table_->getSourceName(whichrow); }
@@ -134 +146 @@
   std::vector<uint> getScanNos() { return table_->getScanNos(); }
   int getScan(int whichrow) const {return table_->getScan(whichrow);}
+  std::vector<uint> getMolNos() { return table_->getMolNos();}
 
   STSelector getSelection() const { return table_->getSelection(); }
@@ -158 +171 @@
   { table_->shift(npix); }
 
+/**
+  commented out by TT
   void setRestFrequencies(double rf, const std::string& name,
                           const std::string& unit)
     { table_->setRestFrequencies(rf, name, unit); }
+**/
+  void setRestFrequencies(vector<double> rf, const vector<std::string>& name,
+                          const std::string& unit)
+    { table_->setRestFrequencies(rf, name, unit); }
+
 /*
   void setRestFrequencies(const std::string& name) {
@@ -167 +187 @@
 */
 
+/*
   std::vector<double> getRestFrequencies() const
     { return table_->getRestFrequencies(); }
+*/
+  std::vector<double> getRestFrequency(int id) const
+    { return table_->getRestFrequency(id); }
 
   void setCoordInfo(std::vector<string> theinfo) {
@@ -209 +233 @@
   int checkScanInfo(const vector<int>& scanlist) const
     { return table_->checkScanInfo(scanlist); }
-
+ 
   std::vector<double> getDirectionVector(int whichrow) const
     { return table_->getDirectionVector(whichrow); }
@@ -222 +246 @@
   std::vector<float> getWeather(int whichrow) const
     { return table_->getWeather(whichrow); }
+
+  void reshapeSpectrum( int nmin, int nmax )
+  { table_->reshapeSpectrum( nmin, nmax ); }
 
 private:
@@ -229 +256 @@
 } // namespace
 #endif
+

trunk/src/python_STMath.cpp

@@ -49 +49 @@
         .def("_averagebeams", &STMathWrapper::averageBeams)
         .def("_unaryop", &STMathWrapper::unaryOperate)
+        .def("_arrayop", &STMathWrapper::arrayOperate)
+        //.def("_array2dop", &STMathWrapper::array2dOperate)
         .def("_binaryop", &STMathWrapper::binaryOperate)
         .def("_auto_quotient", &STMathWrapper::autoQuotient)
@@ -57 +59 @@
         .def("_dofs", &STMathWrapper::dofs)
         .def("_stats", &STMathWrapper::statistic)
+        .def("_minmaxchan", &STMathWrapper::minMaxChan)
         .def("_freqswitch", &STMathWrapper::freqSwitch)
         .def("_bin", &STMathWrapper::bin)
@@ -73 +76 @@
         .def("_mx_extract", &STMathWrapper::mxExtract)
         .def("_lag_flag", &STMathWrapper::lagFlag)
+        // testing average spectra with different channel/resolution
+        .def("_new_average", &STMathWrapper::new_average)
+        // cwcal
+        .def("cwcal", &STMathWrapper::cwcal)
+        .def("almacal", &STMathWrapper::almacal)
           ;
     };

trunk/src/python_STSelector.cpp

@@ -29 +29 @@
         .def("_getscans", &STSelector::getScans)
         .def("_getcycles", &STSelector::getCycles)
+        .def("_gettypes", &STSelector::getTypes)
         .def("_gettaql", &STSelector::getTaQL)
         .def("_getorder", &STSelector::getSortOrder)
@@ -41 +42 @@
         .def("_settaql", &STSelector::setTaQL)
         .def("_setorder", &STSelector::setSortOrder)
+        .def("_setrows", &STSelector::setRows)
+        .def("_settypes", &STSelector::setTypes)
         .def("_empty", &STSelector::empty)
       ;

trunk/src/python_Scantable.cpp

@@ -58 +58 @@
     .def("getscannos", &ScantableWrapper::getScanNos)
     .def("getcycle", &ScantableWrapper::getCycle)
+    .def("getmolnos", &ScantableWrapper::getMolNos)
     .def("nif", &ScantableWrapper::nif,
          (boost::python::arg("scanno")=-1) )
@@ -87 +88 @@
     .def("_getmask", &ScantableWrapper::getMask,
          (boost::python::arg("whichrow")=0) )
+    .def("_getclipmask", &ScantableWrapper::getClipMask,
+         (boost::python::arg("whichrow")=0) )
     .def("_gettsys", &ScantableWrapper::getTsys)
     .def("_getsourcename", &ScantableWrapper::getSourceName,
@@ -103 +106 @@
          (boost::python::arg("whichrow")=0) )
     .def("get_antennaname", &ScantableWrapper::getAntennaName)
-    .def("_flag", &ScantableWrapper::flag,
-                  (boost::python::arg("unflag") = false) )
+    .def("_flag", &ScantableWrapper::flag)
+    .def("_flag_row", &ScantableWrapper::flagRow)
+    .def("_getflagrow", &ScantableWrapper::getFlagRow,
+         (boost::python::arg("whichrow")=0) )
+    .def("_clip", &ScantableWrapper::clip,
+         (boost::python::arg("clipoutside")=true,
+          boost::python::arg("unflag")=false) )
     .def("_save",  &ScantableWrapper::makePersistent)
     .def("_summary",  &ScantableWrapper::summary,
          (boost::python::arg("verbose")=true) )
-    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
+    //.def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
+    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequency)
     .def("_setrestfreqs",  &ScantableWrapper::setRestFrequencies)
     .def("shift_refpix", &ScantableWrapper::shift)
@@ -127 +136 @@
     .def("get_coordinate", &ScantableWrapper::getCoordinate)
     .def("_get_weather", &ScantableWrapper::getWeather)
+    .def("_reshape", &ScantableWrapper::reshapeSpectrum, 
+         (boost::python::arg("nmin")=-1, 
+          boost::python::arg("nmax")=-1) )
   ;
 };

trunk/src/python_asap.cpp

@@ -67 +67 @@
   asap::python::python_Scantable();
   asap::python::python_STFiller();
+  asap::python::python_Filler();
   asap::python::python_STSelector();
   asap::python::python_STMath();
@@ -75 +76 @@
   asap::python::python_LineCatalog();
   asap::python::python_Logger();
+  asap::python::python_LogSink();
   asap::python::python_STCoordinate();
   asap::python::python_STAtmosphere();
+  asap::python::python_SrcType();
 
 #ifndef HAVE_LIBPYRAP

trunk/src/python_asap.h

@@ -39 +39 @@
     void python_Scantable();
     void python_STFiller();
+    void python_Filler();
     void python_STSelector();
     void python_STMath();
@@ -47 +48 @@
     void python_LineCatalog();
     void python_Logger();
+    void python_LogSink();
     void python_STCoordinate();
     void python_STAtmosphere();
+    void python_SrcType();
 
   } // python