Ignore:
Timestamp:
07/30/10 16:59:56 (14 years ago)
Author:
Kana Sugimoto
Message:

Merged the all developments from alma branch. copied external-alma.

Location:
branches/mergetest/src
Files:
30 edited
16 copied

Legend:

Unmodified
Added
Removed
  • branches/mergetest/src

  • branches/mergetest/src/MathUtils.cpp

    r1570 r1789  
    3838#include <casa/BasicSL/String.h>
    3939#include <scimath/Mathematics/MedianSlider.h>
     40#include <casa/Exceptions/Error.h>
    4041
    4142#include <scimath/Fitting/LinearFit.h>
     
    5354   String str(which);
    5455   str.upcase();
    55    if (str.contains(String("MIN"))) {
     56   if (str.matches(String("MIN"))) {
    5657      return min(data);
    57    } else if (str.contains(String("MAX"))) {
     58   } else if (str.matches(String("MAX"))) {
    5859      return max(data);
    59    } else if (str.contains(String("SUMSQ"))) {
     60   } else if (str.matches(String("SUMSQ"))) {
    6061      return sumsquares(data);
    61    } else if (str.contains(String("SUM"))) {
     62   } else if (str.matches(String("SUM"))) {
    6263      return sum(data);
    63    } else if (str.contains(String("MEAN"))) {
     64   } else if (str.matches(String("MEAN"))) {
    6465      return mean(data);
    65    } else if (str.contains(String("VAR"))) {
     66   } else if (str.matches(String("VAR"))) {
    6667      return variance(data);
    67    } else if (str.contains(String("STDDEV"))) {
     68      } else if (str.matches(String("STDDEV"))) {
    6869      return stddev(data);
    69    } else if (str.contains(String("AVDEV"))) {
     70   } else if (str.matches(String("AVDEV"))) {
    7071      return avdev(data);
    71    } else if (str.contains(String("RMS"))) {
     72   } else if (str.matches(String("RMS"))) {
    7273      uInt n = data.nelementsValid();
    7374      return sqrt(sumsquares(data)/n);
    74    } else if (str.contains(String("MED"))) {
     75   } else if (str.matches(String("MEDIAN"))) {
    7576      return median(data);
    76    }
     77   } else {
     78      String msg = str + " is not a valid type of statistics";
     79      throw(AipsError(msg));
     80   }
    7781   return 0.0;
    7882}
    7983
     84IPosition mathutil::minMaxPos(const String& which,
     85                           const MaskedArray<Float>& data)
     86{
     87   Float minVal, maxVal;
     88   IPosition minPos(data.ndim(), 0), maxPos(data.ndim(), 0);
     89   minMax(minVal, maxVal, minPos, maxPos, data);
     90   String str(which);
     91   str.upcase();
     92   if (str.contains(String("MIN"))) {
     93     return minPos;
     94   } else if (str.contains(String("MAX"))) {
     95     return maxPos;
     96   } else {
     97      String msg = str + " is not a valid type of statistics";
     98      throw(AipsError(msg));
     99   }
     100   //return 0.0;
     101}
    80102
    81103void mathutil::replaceMaskByZero(Vector<Float>& data, const Vector<Bool>& mask)
  • branches/mergetest/src/MathUtils.h

    r1570 r1789  
    3737#include <casa/Arrays/Vector.h>
    3838#include <casa/BasicSL/String.h>
     39#include <casa/Arrays/IPosition.h>
    3940
    4041namespace mathutil {
     
    7980                   float hwidth, int order);
    8081
    81 
    8282// Generate specified statistic
    8383float statistics(const casa::String& which,
     84                 const casa::MaskedArray<casa::Float>& data);
     85
     86// Return a position of min or max value
     87 casa::IPosition minMaxPos(const casa::String& which,
    8488                 const casa::MaskedArray<casa::Float>& data);
    8589
  • branches/mergetest/src/RowAccumulator.cpp

    r1569 r1789  
    152152  userMask_ = m;
    153153}
     154
     155// Added by TT  check the state of RowAccumulator
     156casa::Bool RowAccumulator::state() const
     157{
     158  return initialized_;
     159}
     160
  • branches/mergetest/src/RowAccumulator.h

    r1569 r1789  
    8585    */
    8686  void reset();
     87  /**
     88    * check the initialization state
     89    */
     90  casa::Bool state() const;
    8791
    8892private:
  • branches/mergetest/src/SConscript

  • branches/mergetest/src/STAsciiWriter.cpp

    r1552 r1789  
     1
    12//#---------------------------------------------------------------------------
    23//# STAsciiWriter.cc: ASAP class to write out single dish spectra as FITS images
     
    8889
    8990   String rootName(fileName);
    90    
     91
    9192  Block<String> cols(4);
    9293  cols[0] = String("SCANNO");
     
    132133    String wcs = stable.frequencies().print(rec.asuInt("FREQ_ID"), True);
    133134    addLine(of, "WCS", wcs);
    134     addLine(of, "Rest Freq.",
    135             stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
     135    std::vector<double> restfreqs= stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID"));
     136    int nf = restfreqs.size();
     137    //addLine(of, "Rest Freq.",
     138    //        stable.molecules().getRestFrequency(rec.asuInt("MOLECULE_ID") ));
     139    addLine(of, "Rest Freq.", restfreqs[0]);
     140    for ( unsigned int i=1; i<nf; ++i) {
     141      addLine(of, " ", restfreqs[i]);
     142    }
     143    ostringstream osflagrow;
     144    for ( unsigned int i=0; i<t.nrow(); ++i) {
     145      osflagrow << "Pol" << i << ":" << ((row.get(i).asuInt("FLAGROW") > 0) ? "True" : "False") << " ";
     146    }
     147    addLine(of, "Row_Flagged", String(osflagrow));
    136148    of << setfill('#') << setw(70) << "" << setfill(' ') << endl;
    137149
  • branches/mergetest/src/STFiller.cpp

    r1725 r1789  
    2525#include <tables/Tables/TableRow.h>
    2626
     27#include <measures/Measures/MDirection.h>
     28#include <measures/Measures/MeasConvert.h>
     29
    2730#include <atnf/PKSIO/PKSrecord.h>
    2831#include <atnf/PKSIO/PKSreader.h>
     
    3033 #include <casa/System/ProgressMeter.h>
    3134#endif
     35#include <casa/System/ProgressMeter.h>
     36#include <atnf/PKSIO/NROReader.h>
     37#include <casa/Logging/LogIO.h>
     38
     39#include <time.h>
     40
    3241
    3342#include "STDefs.h"
     
    4554  header_(0),
    4655  table_(0),
    47   refRx_(".*(e|w|_R)$")
     56  refRx_(".*(e|w|_R)$"),
     57  nreader_(0)
    4858{
    4959}
     
    5363  header_(0),
    5464  table_(stbl),
    55   refRx_(".*(e|w|_R)$")
     65  refRx_(".*(e|w|_R)$"),
     66  nreader_(0)
    5667{
    5768}
     
    6172  header_(0),
    6273  table_(0),
    63   refRx_(".*(e|w|_R)$")
     74  refRx_(".*(e|w|_R)$"),
     75  nreader_(0)
    6476{
    65   open(filename, whichIF, whichBeam);
     77  open(filename, "", whichIF, whichBeam);
    6678}
    6779
     
    7082  close();
    7183}
    72 
    73 void STFiller::open( const std::string& filename, int whichIF, int whichBeam )
     84 
     85void STFiller::open( const std::string& filename, const std::string& antenna, int whichIF, int whichBeam, casa::Bool getPt )
    7486{
    7587  if (table_.null())  {
     
    93105  Vector<Bool> beams, ifs;
    94106  Vector<uInt> nchans,npols;
    95   if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
     107
     108  //
     109  // if isNRO_ is true, try NROReader
     110  //
     111  // 2008/11/11 Takeshi Nakazato
     112  isNRO_ = fileCheck() ;
     113  if ( isNRO_ ) {
     114    if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
     115      throw(AipsError("Creation of NROReader failed")) ;
     116    }
     117    else {
     118      openNRO( whichIF, whichBeam ) ;
     119      return ;
     120    }
     121  }
     122  //
     123
     124  if ( (reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
    96125                              nchans, npols, haveXPol_,haveBase, haveSpectra
    97126                              )) == 0 )  {
     
    118147  header_->npol = max(npols);
    119148  header_->nbeam = nBeam_;
    120  
     149
    121150  Int status = reader_->getHeader(header_->observer, header_->project,
    122151                                  header_->antennaname, header_->antennaposition,
     
    197226  Vector<Int> start(nIF_, 1);
    198227  Vector<Int> end(nIF_, 0);
    199   reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0]);
     228  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
    200229  table_->setHeader(*header_);
    201230  //For MS, add the location of POINTING of the input MS so one get
     
    216245    }
    217246  }
     247  String freqFrame = header_->freqref;
     248  //translate frequency reference frame back to
     249  //MS style (as PKSMS2reader converts the original frame
     250  //in FITS standard style)
     251  if (freqFrame == "TOPOCENT") {
     252    freqFrame = "TOPO";
     253  } else if (freqFrame == "GEOCENER") {
     254    freqFrame = "GEO";
     255  } else if (freqFrame == "BARYCENT") {
     256    freqFrame = "BARY";
     257  } else if (freqFrame == "GALACTOC") {
     258    freqFrame = "GALACTO";
     259  } else if (freqFrame == "LOCALGRP") {
     260    freqFrame = "LGROUP";
     261  } else if (freqFrame == "CMBDIPOL") {
     262    freqFrame = "CMB";
     263  } else if (freqFrame == "SOURCE") {
     264    freqFrame = "REST";
     265  }
     266  // set both "FRAME" and "BASEFRAME"
     267  table_->frequencies().setFrame(freqFrame, false);
     268  table_->frequencies().setFrame(freqFrame,true);
    218269  //table_->focus().setParallactify(true);
    219270}
     
    222273{
    223274  delete reader_;reader_=0;
     275  delete nreader_;nreader_=0;
    224276  delete header_;header_=0;
    225277  table_ = 0;
     
    229281{
    230282  int status = 0;
     283
     284  //
     285  // for NRO data
     286  //
     287  // 2008/11/12 Takeshi Nakazato
     288  if ( isNRO_ ) {
     289    status = readNRO() ;
     290    return status ;
     291  }
     292  //
     293
     294/**
     295  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
     296  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
     297    humidity, parAngle, pressure, temperature, windAz, windSpeed;
     298  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
     299  String          fieldName, srcName, tcalTime, obsType;
     300  Vector<Float>   calFctr, sigma, tcal, tsys;
     301  Matrix<Float>   baseLin, baseSub;
     302  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
     303  Matrix<Float>   spectra;
     304  Matrix<uChar>   flagtra;
     305  Complex         xCalFctr;
     306  Vector<Complex> xPol;
     307**/
    231308
    232309  Double min = 0.0;
     
    236313#endif
    237314  PKSrecord pksrec;
     315  pksrec.srcType=-1;
    238316  int n = 0;
     317  bool isGBTFITS = false ;
     318  if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
     319    FILE *fp = fopen( filename_.c_str(), "r" ) ;
     320    fseek( fp, 640, SEEK_SET ) ;
     321    char buf[81] ;
     322    fread( buf, 80, 1, fp ) ;
     323    buf[80] = '\0' ;
     324    if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
     325      isGBTFITS = true ;
     326    }
     327    fclose( fp ) ;
     328  }
    239329  while ( status == 0 ) {
    240330    status = reader_->read(pksrec);
     
    288378    //*srcnCol = pksrec.srcName;//.before(rx2);
    289379    *srctCol = match;
     380    if ( pksrec.srcType != -1 ) {
     381      *srctCol = pksrec.srcType ;
     382    }
    290383    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
    291384    *beamCol = pksrec.beamNo-beamOffset_-1;
     
    298391    uInt id;
    299392    /// @todo this has to change when nchan isn't global anymore
    300     id = table_->frequencies().addEntry(Double(header_->nchan/2),
    301                                         pksrec.refFreq, pksrec.freqInc);
     393    //id = table_->frequencies().addEntry(Double(header_->nchan/2),
     394    //                                    pksrec.refFreq, pksrec.freqInc);
     395    if ( pksrec.nchan == 1 ) {
     396      id = table_->frequencies().addEntry(Double(0),
     397                                          pksrec.refFreq, pksrec.freqInc);
     398    }
     399    else {
     400      id = table_->frequencies().addEntry(Double(pksrec.nchan/2),
     401                                          pksrec.refFreq, pksrec.freqInc);
     402    }
    302403    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
    303404    *mfreqidCol = id;
     405    //*ifCol = id;
    304406
    305407    id = table_->molecules().addEntry(pksrec.restFreq);
     
    335437    // into 2-4 rows in the scantable
    336438    Vector<Float> tsysvec(1);
    337     // Why is pksrec.spectra.ncolumn() == 3 for haveXPol_ == True
     439    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
    338440    uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
    339441    for ( uInt i=0; i< npol; ++i ) {
    340442      tsysvec = pksrec.tsys(i);
    341443      *tsysCol = tsysvec;
    342       *polnoCol = i;
     444      if (isGBTFITS)
     445        *polnoCol = pksrec.polNo ;
     446      else
     447        *polnoCol = i;
    343448
    344449      *specCol = pksrec.spectra.column(i);
     
    347452      row.put(table_->table().nrow()-1, rec);
    348453    }
     454
     455    RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
     456    *flagrowCol = pksrec.flagrow;
     457
    349458    if ( haveXPol_[0] ) {
    350459      // no tsys given for xpol, so emulate it
     
    381490}
    382491
     492/**
     493 * For NRO data
     494 *
     495 * 2008/11/11 Takeshi Nakazato
     496 **/
     497void STFiller::openNRO( int whichIF, int whichBeam )
     498{
     499  // open file
     500  // DEBUG
     501  time_t t0 ;
     502  time( &t0 ) ;
     503  tm *ttm = localtime( &t0 ) ;
     504  LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
     505//   cout << "STFiller::openNRO()  Start time = " << t0
     506//        << " ("
     507//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     508//        << " "
     509//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     510//        << ")" << endl ;
     511  os << "Start time = " << t0
     512     << " ("
     513     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     514     << " "
     515     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     516     << ")" << LogIO::POST ;
     517
     518  // fill STHeader
     519  header_ = new STHeader() ;
     520
     521  if ( nreader_->getHeaderInfo( header_->nchan,
     522                                header_->npol,
     523                                nIF_,
     524                                nBeam_,
     525                                header_->observer,
     526                                header_->project,
     527                                header_->obstype,
     528                                header_->antennaname,
     529                                header_->antennaposition,
     530                                header_->equinox,
     531                                header_->freqref,
     532                                header_->reffreq,
     533                                header_->bandwidth,
     534                                header_->utc,
     535                                header_->fluxunit,
     536                                header_->epoch,
     537                                header_->poltype ) ) {
     538//     cout << "STFiller::openNRO()  Failed to get header information." << endl ;
     539//     return ;
     540    throw( AipsError("Failed to get header information.") ) ;
     541  }
     542
     543  // set FRAME and BASEFRAME keyword of FREQUENCIES table
     544  if ( header_->freqref != "TOPO" ) {
     545    table_->frequencies().setFrame( header_->freqref, false ) ;
     546    table_->frequencies().setFrame( header_->freqref, true ) ;
     547  }
     548
     549  ifOffset_ = 0;
     550  vector<Bool> ifs = nreader_->getIFs() ;
     551  if ( whichIF >= 0 ) {
     552    if ( whichIF >= 0 && whichIF < nIF_ ) {
     553      for ( int i = 0 ; i < nIF_ ; i++ )
     554        ifs[i] = False ;
     555      ifs[whichIF] = True ;
     556      header_->nif = 1;
     557      nIF_ = 1;
     558      ifOffset_ = whichIF;
     559    } else {
     560      delete reader_;
     561      reader_ = 0;
     562      delete header_;
     563      header_ = 0;
     564      throw(AipsError("Illegal IF selection"));
     565    }
     566  }
     567
     568  beamOffset_ = 0;
     569  vector<Bool> beams = nreader_->getBeams() ;
     570  if (whichBeam>=0) {
     571    if (whichBeam>=0 && whichBeam<nBeam_) {
     572      for ( int i = 0 ; i < nBeam_ ; i++ )
     573        beams[i] = False ;
     574      beams[whichBeam] = True;
     575      header_->nbeam = 1;
     576      nBeam_ = 1;
     577      beamOffset_ = whichBeam;
     578    } else {
     579      delete reader_;
     580      reader_ = 0;
     581      delete header_;
     582      header_ = 0;
     583      throw(AipsError("Illegal Beam selection"));
     584    }
     585  }
     586
     587  header_->nbeam = nBeam_ ;
     588  header_->nif = nIF_ ;
     589
     590  // set header
     591  table_->setHeader( *header_ ) ;
     592
     593  // DEBUG
     594  time_t t1 ;
     595  time( &t1 ) ;
     596  ttm = localtime( &t1 ) ;
     597//   cout << "STFiller::openNRO()  End time = " << t1
     598//        << " ("
     599//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     600//        << " "
     601//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     602//        << ")" << endl ;
     603//   cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
     604  os << "End time = " << t1
     605     << " ("
     606     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     607     << " "
     608     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     609     << ")" << endl ;
     610  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
     611  os.post() ;
     612  //
     613
     614  return ;
     615}
     616
     617int STFiller::readNRO()
     618{
     619  // DEBUG
     620  time_t t0 ;
     621  time( &t0 ) ;
     622  tm *ttm = localtime( &t0 ) ;
     623  LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
     624//   cout << "STFiller::readNRO()  Start time = " << t0
     625//        << " ("
     626//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     627//        << " "
     628//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     629//        << ")" << endl ;
     630  os << "Start time = " << t0
     631     << " ("
     632     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     633     << " "
     634     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     635     << ")" << LogIO::POST ;
     636  //
     637
     638  // fill row
     639  uInt id ;
     640  uInt imax = nreader_->getRowNum() ;
     641  vector< vector<double > > freqs ;
     642  uInt i = 0 ;
     643  int count = 0 ;
     644  uInt scanno ;
     645  uInt cycleno ;
     646  uInt beamno ;
     647  uInt polno ;
     648  vector<double> fqs ;
     649  Vector<Double> restfreq ;
     650  uInt refbeamno ;
     651  Double scantime ;
     652  Double interval ;
     653  String srcname ;
     654  String fieldname ;
     655  Array<Float> spectra ;
     656  Array<uChar> flagtra ;
     657  Array<Float> tsys ;
     658  Array<Double> direction ;
     659  Float azimuth ;
     660  Float elevation ;
     661  Float parangle ;
     662  Float opacity ;
     663  uInt tcalid ;
     664  Int fitid ;
     665  uInt focusid ;
     666  Float temperature ;
     667  Float pressure ;
     668  Float humidity ;
     669  Float windvel ;
     670  Float winddir ;
     671  Double srcvel ;
     672  Array<Double> propermotion ;
     673  Vector<Double> srcdir ;
     674  Array<Double> scanrate ;
     675  for ( i = 0 ; i < imax ; i++ ) {
     676    string scanType = nreader_->getScanType( i ) ;
     677    Int srcType = -1 ;
     678    if ( scanType.compare( 0, 2, "ON") == 0 ) {
     679      // os << "ON srcType: " << i << LogIO::POST ;
     680      srcType = 0 ;
     681    }
     682    else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
     683      //os << "OFF srcType: " << i << LogIO::POST ;
     684      srcType = 1 ;
     685    }
     686    else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
     687      //os << "ZERO srcType: " << i << LogIO::POST ;
     688      srcType = 2 ;
     689    }
     690    else {
     691      //os << "Undefined srcType: " << i << LogIO::POST ;
     692      srcType = 3 ;
     693    }
     694 
     695    // if srcType is 2 (ZERO scan), ignore scan
     696    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
     697      TableRow row( table_->table() ) ;
     698      TableRecord& rec = row.record();
     699
     700      if ( nreader_->getScanInfo( i,
     701                                  scanno,
     702                                  cycleno,
     703                                  beamno,
     704                                  polno,
     705                                  fqs,
     706                                  restfreq,
     707                                  refbeamno,
     708                                  scantime,
     709                                  interval,
     710                                  srcname,
     711                                  fieldname,
     712                                  spectra,
     713                                  flagtra,
     714                                  tsys,
     715                                  direction,
     716                                  azimuth,
     717                                  elevation,
     718                                  parangle,
     719                                  opacity,
     720                                  tcalid,
     721                                  fitid,
     722                                  focusid,
     723                                  temperature,
     724                                  pressure,
     725                                  humidity,
     726                                  windvel,
     727                                  winddir,
     728                                  srcvel,
     729                                  propermotion,
     730                                  srcdir,
     731                                  scanrate ) ) {
     732//         cerr << "STFiller::readNRO()  Failed to get scan information." << endl ;
     733//         return 1 ;
     734        throw( AipsError("Failed to get scan information.") ) ;
     735      }
     736
     737      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
     738      *scannoCol = scanno ;
     739      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
     740      *cyclenoCol = cycleno ;
     741      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
     742      *beamCol = beamno ;
     743      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
     744      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
     745      *polnoCol = polno ;
     746      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
     747      if ( freqs.size() == 0 ) {
     748        id = table_->frequencies().addEntry( Double( fqs[0] ),
     749                                             Double( fqs[1] ),
     750                                             Double( fqs[2] ) ) ;
     751        *mfreqidCol = id ;
     752        *ifCol = id ;
     753        freqs.push_back( fqs ) ;
     754      }
     755      else {
     756        int iadd = -1 ;
     757        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
     758          //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
     759          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
     760          //os << "fdiff = " << fdiff << LogIO::POST ;
     761          if ( fdiff < 1.0e-8 ) {
     762            iadd = iif ;
     763            break ;
     764          }
     765        }
     766        if ( iadd == -1 ) {
     767          id = table_->frequencies().addEntry( Double( fqs[0] ),
     768                                               Double( fqs[1] ),
     769                                               Double( fqs[2] ) ) ;
     770          *mfreqidCol = id ;
     771          *ifCol = id ;
     772          freqs.push_back( fqs ) ;
     773        }
     774        else {
     775          *mfreqidCol = iadd ;
     776          *ifCol = iadd ;
     777        }
     778      }
     779      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
     780      id = table_->molecules().addEntry( restfreq ) ;
     781      *molidCol = id ;
     782      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
     783      *rbCol = refbeamno ;
     784      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
     785      *mjdCol = scantime ;
     786      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
     787      *intervalCol = interval ;
     788      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
     789      *srcnCol = srcname ;
     790      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
     791      *srctCol = srcType ;     
     792      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
     793      *fieldnCol = fieldname ;
     794      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
     795      *specCol = spectra ;
     796      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
     797      *flagCol = flagtra ;
     798      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
     799      *tsysCol = tsys ;
     800      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
     801      *dirCol = direction ;
     802      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
     803      *azCol = azimuth ;
     804      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
     805      *elCol = elevation ;
     806      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
     807      *parCol = parangle ;
     808      RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
     809      *tauCol = opacity ;
     810      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
     811      *mcalidCol = tcalid ;
     812      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
     813      *fitCol = fitid ;
     814      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
     815      *mfocusidCol = focusid ;
     816      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
     817      id = table_->weather().addEntry( temperature,
     818                                       pressure,
     819                                       humidity,
     820                                       windvel,
     821                                       winddir ) ;
     822      *mweatheridCol = id ;         
     823      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
     824      *svelCol = srcvel ;
     825      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
     826      *spmCol = propermotion ;
     827      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
     828      *sdirCol = srcdir ;
     829      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
     830      *srateCol = scanrate ;
     831
     832      table_->table().addRow() ;
     833      row.put(table_->table().nrow()-1, rec) ;
     834    }
     835    else {
     836      count++ ;
     837    }
     838    // DEBUG
     839    //int rownum = nreader_->getRowNum() ;
     840    //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
     841    //
     842  }
     843
     844  // DEBUG
     845  time_t t1 ;
     846  time( &t1 ) ;
     847  ttm = localtime( &t1 ) ;
     848//   cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
     849//   cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
     850//        << count << " \"ZERO\" scans)" << endl ;
     851//   cout << "STFiller::readNRO()  End time = " << t1
     852//        << " ("
     853//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     854//        << " "
     855//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     856//        << ")" << endl ;
     857//   cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
     858  os << "Processed " << i << " rows" << endl ;
     859  os << "Added " << i - count << " rows (ignored "
     860     << count << " \"ZERO\" scans)" << endl ;
     861  os.post() ;
     862  os << "End time = " << t1
     863     << " ("
     864     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
     865     << " "
     866     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
     867     << ")" << endl ;
     868  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
     869  os.post() ;
     870  //
     871
     872  return 0 ;
     873}
     874
     875Bool STFiller::fileCheck()
     876{
     877  bool bval = false ;
     878
     879  // if filename_ is directory, return false
     880  File inFile( filename_ ) ;
     881  if ( inFile.isDirectory() )
     882    return bval ;
     883 
     884  // if beginning of header data is "RW", return true
     885  // otherwise, return false ;
     886  FILE *fp = fopen( filename_.c_str(), "r" ) ;
     887  char buf[9] ;
     888  char buf2[80] ;
     889  fread( buf, 4, 1, fp ) ;
     890  buf[4] = '\0' ;
     891  fseek( fp, 640, SEEK_SET ) ;
     892  fread( buf2, 80, 1, fp ) ;
     893  if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
     894    bval = true ;
     895  }
     896  fclose( fp ) ;
     897  return bval ;
     898}
     899
    383900}//namespace asap
  • branches/mergetest/src/STFiller.h

    r1504 r1789  
    2727
    2828class PKSreader;
     29class NROReader;
    2930
    3031namespace asap {
     
    6162    */
    6263  explicit STFiller( const std::string& filename, int whichIF=-1,
    63                   int whichBeam=-1 );
     64                     int whichBeam=-1 );
    6465
    6566  /**
     
    7576   * @exception AipsError Creation of PKSreader failed
    7677   */
    77   void open( const std::string& filename, int whichIF=-1, int whichBeam=-1 );
     78  void open( const std::string& filename, const std::string& antenna, int whichIF=-1, int whichBeam=-1, casa::Bool getPt=casa::False );
    7879
    7980  /**
     
    9394  casa::CountedPtr<Scantable> getTable() const { return table_;}
    9495
     96  /**
     97   * For NRO data
     98   *
     99   * 2008/11/11 Takeshi Nakazato
     100   *
     101   * openNRO  : NRO version of open(), which performs to open file and
     102   *            read header data.
     103   * 
     104   * readNRO  : NRO version of read(), which performs to read scan
     105   *            records.
     106   *
     107   * fileCheck: Identify a type (NRO data or not) of filename_.
     108   **/
     109  void openNRO( int whichIF=-1, int whichBeam=-1 ) ;
     110  int readNRO() ;
     111  casa::Bool fileCheck() ;
     112
    95113  void setReferenceExpr(const std::string& rx) { refRx_ = rx; }
    96114
     
    105123  casa::Vector<casa::Bool> haveXPol_;
    106124  casa::String refRx_;
     125  NROReader *nreader_ ;
     126  casa::Bool isNRO_ ;
    107127};
    108128
  • branches/mergetest/src/STFitter.cpp

    r1391 r1789  
    3232#include <casa/Arrays/ArrayMath.h>
    3333#include <casa/Arrays/ArrayLogical.h>
     34#include <casa/Logging/LogIO.h>
    3435#include <scimath/Fitting.h>
    3536#include <scimath/Fitting/LinearFit.h>
     
    3738#include <scimath/Functionals/CompoundFunction.h>
    3839#include <scimath/Functionals/Gaussian1D.h>
     40#include "Lorentzian1D.h"
    3941#include <scimath/Functionals/Polynomial.h>
    4042#include <scimath/Mathematics/AutoDiff.h>
     
    146148    funcs_.resize(1);
    147149    funcs_[0] = new Polynomial<Float>(ncomp);
     150  } else if (expr == "lorentz") {
     151    if (ncomp < 1) throw (AipsError("Need at least one lorentzian to fit."));
     152    funcs_.resize(ncomp);
     153    for (Int k=0; k<ncomp; ++k) {
     154      funcs_[k] = new Lorentzian1D<Float>();
     155    }
    148156  } else {
    149     cerr << " compiled functions not yet implemented" << endl;
     157    //cerr << " compiled functions not yet implemented" << endl;
     158    LogIO os( LogOrigin( "Fitter", "setExpression()", WHERE ) ) ;
     159    os << LogIO::WARN << " compiled functions not yet implemented" << LogIO::POST;
    150160    //funcs_.resize(1);
    151161    //funcs_[0] = new CompiledFunction<Float>();
     
    227237            (funcs_[0]->parameters())[i] =  tmppar[i];
    228238        }
     239    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
     240        uInt count = 0;
     241        for (uInt j=0; j < funcs_.nelements(); ++j) {
     242            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     243                (funcs_[j]->parameters())[i] = tmppar[count];
     244                parameters_[count] = tmppar[count];
     245                ++count;
     246            }
     247        }
    229248    }
    230249    // reset
     
    260279            funcs_[0]->mask(i) =  !fixed[i];
    261280        }
     281    } else if (dynamic_cast<Lorentzian1D<Float>* >(funcs_[0]) != 0) {
     282      uInt count = 0;
     283        for (uInt j=0; j < funcs_.nelements(); ++j) {
     284            for (uInt i=0; i < funcs_[j]->nparameters(); ++i) {
     285                funcs_[j]->mask(i) = !fixed[count];
     286                fixedpar_[count] = fixed[count];
     287                ++count;
     288            }
     289        }
    262290    }
    263291    return true;
  • branches/mergetest/src/STFrequencies.cpp

    r1694 r1789  
    128128}
    129129
     130void STFrequencies::setEntry( Double refpix, Double refval, Double inc, uInt id )
     131{
     132  Table t = table_(table_.col("ID") == Int(id) );
     133  if (t.nrow() == 0 ) {
     134    throw(AipsError("STFrequencies::getEntry - freqID out of range"));
     135  }
     136  for ( uInt i = 0 ; i < table_.nrow() ; i++ ) {
     137    uInt fid ;
     138    idCol_.get( i, fid ) ;
     139    if ( fid == id ) {
     140      refpixCol_.put( i, refpix ) ;
     141      refvalCol_.put( i, refval ) ;
     142      incrCol_.put( i, inc ) ;
     143    }
     144  }
     145}
     146
    130147SpectralCoordinate STFrequencies::getSpectralCoordinate( uInt id ) const
    131148{
     
    145162}
    146163
     164/**
    147165SpectralCoordinate
    148166  STFrequencies::getSpectralCoordinate( const MDirection& md,
     
    150168                                        const MEpoch& me,
    151169                                        Double restfreq, uInt id ) const
     170**/
     171SpectralCoordinate
     172  STFrequencies::getSpectralCoordinate( const MDirection& md,
     173                                              const MPosition& mp,
     174                                              const MEpoch& me,
     175                                              Vector<Double> restfreq, uInt id ) const
    152176{
    153177  SpectralCoordinate spc = getSpectralCoordinate(id);
    154   spc.setRestFrequency(restfreq, True);
     178  //spc.setRestFrequency(restfreq, True);
     179  // for now just use the first rest frequency
     180  if (restfreq.nelements()==0 ) {
     181    restfreq.resize(1);
     182    restfreq[0] = 0;
     183  }
     184  spc.setRestFrequency(restfreq[0], True);
    155185  if ( !spc.setReferenceConversion(getFrame(), me, mp, md) ) {
    156186    throw(AipsError("Couldn't convert frequency frame."));
  • branches/mergetest/src/STFrequencies.h

    r1375 r1789  
    5959                 casa::Double& inc, casa::uInt id );
    6060
     61  /***
     62   * Set the frequency values for a specific id via references
     63   * @param refpix the reference pixel
     64   * @param refval the reference value
     65   * @param inc    the increment
     66   * @param id     the identifier
     67   *
     68   * 17/09/2008 Takeshi Nakazato
     69   ***/
     70  void setEntry( casa::Double refpix, casa::Double refval,
     71                 casa::Double inc, casa::uInt id ) ;
     72
    6173
    6274  bool conformant(const STFrequencies& other) const;
     
    6981  casa::SpectralCoordinate getSpectralCoordinate( casa::uInt freqID ) const;
    7082
     83  /**
    7184  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
    7285                                                  const casa::MPosition& mp,
     
    7588                                                  casa::uInt freqID
    7689                                                  ) const;
     90  **/
     91  casa::SpectralCoordinate getSpectralCoordinate( const casa::MDirection& md,
     92                                                  const casa::MPosition& mp,
     93                                                  const casa::MEpoch& me,
     94                                                  casa::Vector<casa::Double> restfreq,
     95                                                  casa::uInt freqID
     96                                                  ) const;
    7797
    7898  /**
  • branches/mergetest/src/STHeader.cpp

    r1439 r1789  
    3737#include <casa/Arrays/IPosition.h>
    3838#include <casa/Quanta/MVTime.h>
     39#include <casa/Logging/LogIO.h>
    3940
    40 
     41#include <sstream>
    4142
    4243#include "STDefs.h"
     
    7980  MVTime mvt(this->utc);
    8081  mvt.setFormat(MVTime::YMD);
    81   cout << "Observer: " << this->observer << endl
    82        << "Project: " << this->project << endl
    83        << "Obstype: " << this->obstype << endl
    84        << "Antenna: " << this->antennaname << endl
    85        << "Ant. Position: " << this->antennaposition << endl
    86        << "Equinox: " << this->equinox << endl
    87        << "Freq. ref.: " << this->freqref << endl
    88        << "Ref. frequency: " << this->reffreq << endl
    89        << "Bandwidth: "  << this->bandwidth << endl
    90        << "Time (utc): "
    91        << mvt
    92        << endl;
     82//   cout << "Observer: " << this->observer << endl
     83//        << "Project: " << this->project << endl
     84//        << "Obstype: " << this->obstype << endl
     85//        << "Antenna: " << this->antennaname << endl
     86//        << "Ant. Position: " << this->antennaposition << endl
     87//        << "Equinox: " << this->equinox << endl
     88//        << "Freq. ref.: " << this->freqref << endl
     89//        << "Ref. frequency: " << this->reffreq << endl
     90//        << "Bandwidth: "  << this->bandwidth << endl
     91//        << "Time (utc): "
     92//        << mvt
     93//        << endl;
     94  LogIO os( LogOrigin( "STHeader", "print()", WHERE ) ) ;
     95  os << "Observer: " << this->observer << endl
     96     << "Project: " << this->project << endl
     97     << "Obstype: " << this->obstype << endl
     98     << "Antenna: " << this->antennaname << endl
     99     << "Ant. Position: " << this->antennaposition << endl
     100     << "Equinox: " << this->equinox << endl
     101     << "Freq. ref.: " << this->freqref << endl
     102     << "Ref. frequency: " << this->reffreq << endl
     103     << "Bandwidth: "  << this->bandwidth << endl
     104     << "Time (utc): "
     105     << mvt
     106     << LogIO::POST ;
    93107  //setprecision(10) << this->utc << endl;
    94108}
     
    129143{
    130144   if (n_>0) {
    131       cerr << "Source    ID" << endl;
    132       for (uInt i=0; i<n_; i++) {
    133          cout << setw(11) << source_(i) << ID_(i) << endl;
    134       }
     145//       cerr << "Source    ID" << endl;
     146//       for (uInt i=0; i<n_; i++) {
     147//          cout << setw(11) << source_(i) << ID_(i) << endl;
     148     LogIO os( LogOrigin( "SDDataDesc", "summary()", WHERE ) ) ;
     149     ostringstream oss ;
     150     oss << "Source    ID" << endl;
     151     for (uInt i=0; i<n_; i++) {
     152       oss << setw(11) << source_(i) << ID_(i) << endl;
     153     }
     154     os << oss.str() << LogIO::POST ;
    135155   }
    136156}
  • branches/mergetest/src/STMath.cpp

    r1689 r1789  
    4444#include <scimath/Functionals/Polynomial.h>
    4545
     46#include <atnf/PKSIO/SrcType.h>
     47
     48#include <casa/Logging/LogIO.h>
     49#include <sstream>
     50
    4651#include "MathUtils.h"
    4752#include "RowAccumulator.h"
     
    5358
    5459using namespace asap;
     60
     61// tolerance for direction comparison (rad)
     62#define TOL_OTF    1.0e-15
     63#define TOL_POINT  2.9088821e-4  // 1 arcmin
    5564
    5665STMath::STMath(bool insitu) :
     
    7079                 const std::string& avmode)
    7180{
     81  LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
    7282  if ( avmode == "SCAN" && in.size() != 1 )
    7383    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
    7484                    "Use merge first."));
    7585  WeightType wtype = stringToWeight(weight);
     86
     87  // check if OTF observation
     88  String obstype = in[0]->getHeader().obstype ;
     89  Double tol = 0.0 ;
     90  if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
     91    tol = TOL_OTF ;
     92  }
     93  else {
     94    tol = TOL_POINT ;
     95  }
    7696
    7797  // output
     
    113133  }
    114134  if ( avmode == "SCAN"  && in.size() == 1) {
    115     cols.resize(4);
    116     cols[3] = String("SCANNO");
     135    //cols.resize(4);
     136    //cols[3] = String("SCANNO");
     137    cols.resize(5);
     138    cols[3] = String("SRCNAME");
     139    cols[4] = String("SCANNO");
    117140  }
    118141  uInt outrowCount = 0;
    119142  TableIterator iter(baset, cols);
     143//   int count = 0 ;
    120144  while (!iter.pastEnd()) {
    121145    Table subt = iter.table();
    122     // copy the first row of this selection into the new table
    123     tout.addRow();
    124     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    125     // re-index to 0
    126     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    127       scanColOut.put(outrowCount, uInt(0));
    128     }
    129     ++outrowCount;
     146//     // copy the first row of this selection into the new table
     147//     tout.addRow();
     148//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     149//     // re-index to 0
     150//     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     151//       scanColOut.put(outrowCount, uInt(0));
     152//     }
     153//     ++outrowCount;
     154    MDirection::ScalarColumn dircol ;
     155    dircol.attach( subt, "DIRECTION" ) ;
     156    Int length = subt.nrow() ;
     157    vector< Vector<Double> > dirs ;
     158    vector<int> indexes ;
     159    for ( Int i = 0 ; i < length ; i++ ) {
     160      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     161      //os << << count++ << ": " ;
     162      //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
     163      bool adddir = true ;
     164      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     165        //if ( allTrue( t == dirs[j] ) ) {
     166        Double dx = t[0] - dirs[j][0] ;
     167        Double dy = t[1] - dirs[j][1] ;
     168        Double dd = sqrt( dx * dx + dy * dy ) ;
     169        //if ( allNearAbs( t, dirs[j], tol ) ) {
     170        if ( dd <= tol ) {
     171          adddir = false ;
     172          break ;
     173        }
     174      }
     175      if ( adddir ) {
     176        dirs.push_back( t ) ;
     177        indexes.push_back( i ) ;
     178      }
     179    }
     180    uInt rowNum = dirs.size() ;
     181    tout.addRow( rowNum ) ;
     182    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     183      TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
     184      // re-index to 0
     185      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     186        scanColOut.put(outrowCount+i, uInt(0));
     187      }       
     188    }
     189    outrowCount += rowNum ;
    130190    ++iter;
    131191  }
    132 
    133192  RowAccumulator acc(wtype);
    134193  Vector<Bool> cmask(mask);
     
    155214        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
    156215      } else if (avmode == "SCAN") {
    157         subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
     216        //subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
     217        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO"))
     218                         && basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
    158219      } else {
    159220        subt = basesubt;
    160221      }
     222
     223      vector<uInt> removeRows ;
     224      uInt nrsubt = subt.nrow() ;
     225      for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
     226        //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
     227        Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
     228        Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
     229        double dx = x0[0] - x1[0] ;
     230        double dy = x0[0] - x1[0] ;
     231        Double dd = sqrt( dx * dx + dy * dy ) ;
     232        //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
     233        if ( dd > tol ) {
     234          removeRows.push_back( irow ) ;
     235        }
     236      }
     237      if ( removeRows.size() != 0 ) {
     238        subt.removeRow( removeRows ) ;
     239      }
     240     
     241      if ( nrsubt == removeRows.size() )
     242        throw(AipsError("Averaging data is empty.")) ;
     243
    161244      specCol.attach(subt,"SPECTRA");
    162245      flagCol.attach(subt,"FLAGTRA");
     
    192275    }
    193276    //write out
    194     Vector<uChar> flg(msk.shape());
    195     convertArray(flg, !msk);
    196     flagColOut.put(i, flg);
    197     specColOut.put(i, acc.getSpectrum());
    198     tsysColOut.put(i, acc.getTsys());
    199     intColOut.put(i, acc.getInterval());
    200     mjdColOut.put(i, acc.getTime());
    201     // we should only have one cycle now -> reset it to be 0
    202     // frequency switched data has different CYCLENO for different IFNO
    203     // which requires resetting this value
    204     cycColOut.put(i, uInt(0));
     277    if (acc.state()) {
     278      Vector<uChar> flg(msk.shape());
     279      convertArray(flg, !msk);
     280      flagColOut.put(i, flg);
     281      specColOut.put(i, acc.getSpectrum());
     282      tsysColOut.put(i, acc.getTsys());
     283      intColOut.put(i, acc.getInterval());
     284      mjdColOut.put(i, acc.getTime());
     285      // we should only have one cycle now -> reset it to be 0
     286      // frequency switched data has different CYCLENO for different IFNO
     287      // which requires resetting this value
     288      cycColOut.put(i, uInt(0));
     289    } else {
     290      ostringstream oss;
     291      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
     292      pushLog(String(oss));
     293    }
    205294    acc.reset();
    206295  }
    207296  if (rowstodelete.nelements() > 0) {
     297    //cout << rowstodelete << endl;
     298    os << rowstodelete << LogIO::POST ;
    208299    tout.removeRow(rowstodelete);
    209300    if (tout.nrow() == 0) {
     
    219310                          const std::string& avmode )
    220311{
     312  // check if OTF observation
     313  String obstype = in->getHeader().obstype ;
     314  Double tol = 0.0 ;
     315  if ( obstype.find( "OTF" ) != String::npos ) {
     316    tol = TOL_OTF ;
     317  }
     318  else {
     319    tol = TOL_POINT ;
     320  }
     321
    221322  // clone as this is non insitu
    222323  bool insitu = insitu_;
     
    250351    flagCol.attach(subt,"FLAGTRA");
    251352    tsysCol.attach(subt,"TSYS");
    252     tout.addRow();
    253     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    254     if ( avmode != "SCAN") {
    255       scanColOut.put(outrowCount, uInt(0));
    256     }
    257     Vector<Float> tmp;
    258     specCol.get(0, tmp);
    259     uInt nchan = tmp.nelements();
    260     // have to do channel by channel here as MaskedArrMath
    261     // doesn't have partialMedians
    262     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
    263     Vector<Float> outspec(nchan);
    264     Vector<uChar> outflag(nchan,0);
    265     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
    266     for (uInt i=0; i<nchan; ++i) {
    267       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
    268       MaskedArray<Float> ma = maskedArray(specs,flags);
    269       outspec[i] = median(ma);
    270       if ( allEQ(ma.getMask(), False) )
    271         outflag[i] = userflag;// flag data
    272     }
    273     outtsys[0] = median(tsysCol.getColumn());
    274     specColOut.put(outrowCount, outspec);
    275     flagColOut.put(outrowCount, outflag);
    276     tsysColOut.put(outrowCount, outtsys);
    277     Double intsum = sum(intCol.getColumn());
    278     intColOut.put(outrowCount, intsum);
    279     ++outrowCount;
     353//     tout.addRow();
     354//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     355//     if ( avmode != "SCAN") {
     356//       scanColOut.put(outrowCount, uInt(0));
     357//     }
     358//     Vector<Float> tmp;
     359//     specCol.get(0, tmp);
     360//     uInt nchan = tmp.nelements();
     361//     // have to do channel by channel here as MaskedArrMath
     362//     // doesn't have partialMedians
     363//     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     364//     Vector<Float> outspec(nchan);
     365//     Vector<uChar> outflag(nchan,0);
     366//     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     367//     for (uInt i=0; i<nchan; ++i) {
     368//       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     369//       MaskedArray<Float> ma = maskedArray(specs,flags);
     370//       outspec[i] = median(ma);
     371//       if ( allEQ(ma.getMask(), False) )
     372//         outflag[i] = userflag;// flag data
     373//     }
     374//     outtsys[0] = median(tsysCol.getColumn());
     375//     specColOut.put(outrowCount, outspec);
     376//     flagColOut.put(outrowCount, outflag);
     377//     tsysColOut.put(outrowCount, outtsys);
     378//     Double intsum = sum(intCol.getColumn());
     379//     intColOut.put(outrowCount, intsum);
     380//     ++outrowCount;
     381//     ++iter;
     382    MDirection::ScalarColumn dircol ;
     383    dircol.attach( subt, "DIRECTION" ) ;
     384    Int length = subt.nrow() ;
     385    vector< Vector<Double> > dirs ;
     386    vector<int> indexes ;
     387    for ( Int i = 0 ; i < length ; i++ ) {
     388      Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     389      bool adddir = true ;
     390      for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     391        //if ( allTrue( t == dirs[j] ) ) {
     392        Double dx = t[0] - dirs[j][0] ;
     393        Double dy = t[1] - dirs[j][1] ;
     394        Double dd = sqrt( dx * dx + dy * dy ) ;
     395        //if ( allNearAbs( t, dirs[j], tol ) ) {
     396        if ( dd <= tol ) {
     397          adddir = false ;
     398          break ;
     399        }
     400      }
     401      if ( adddir ) {
     402        dirs.push_back( t ) ;
     403        indexes.push_back( i ) ;
     404      }
     405    }
     406    uInt rowNum = dirs.size() ;
     407    tout.addRow( rowNum );
     408    for ( uInt i = 0 ; i < rowNum ; i++ ) {
     409      TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
     410      if ( avmode != "SCAN") {
     411        //scanColOut.put(outrowCount+i, uInt(0));
     412      }
     413    }
     414    MDirection::ScalarColumn dircolOut ;
     415    dircolOut.attach( tout, "DIRECTION" ) ;
     416    for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
     417      Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
     418      Vector<Float> tmp;
     419      specCol.get(0, tmp);
     420      uInt nchan = tmp.nelements();
     421      // have to do channel by channel here as MaskedArrMath
     422      // doesn't have partialMedians
     423      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     424      // mask spectra for different DIRECTION
     425      for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
     426        Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     427        //if ( t[0] != direction[0] || t[1] != direction[1] ) {
     428        Double dx = t[0] - direction[0] ;
     429        Double dy = t[1] - direction[1] ;
     430        Double dd = sqrt( dx * dx + dy * dy ) ;
     431        //if ( !allNearAbs( t, direction, tol ) ) {
     432        if ( dd > tol ) {
     433          flags[jrow] = userflag ;
     434        }
     435      }
     436      Vector<Float> outspec(nchan);
     437      Vector<uChar> outflag(nchan,0);
     438      Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     439      for (uInt i=0; i<nchan; ++i) {
     440        Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     441        MaskedArray<Float> ma = maskedArray(specs,flags);
     442        outspec[i] = median(ma);
     443        if ( allEQ(ma.getMask(), False) )
     444          outflag[i] = userflag;// flag data
     445      }
     446      outtsys[0] = median(tsysCol.getColumn());
     447      specColOut.put(outrowCount+irow, outspec);
     448      flagColOut.put(outrowCount+irow, outflag);
     449      tsysColOut.put(outrowCount+irow, outtsys);
     450      Vector<Double> integ = intCol.getColumn() ;
     451      MaskedArray<Double> mi = maskedArray( integ, flags ) ;
     452      Double intsum = sum(mi);
     453      intColOut.put(outrowCount+irow, intsum);
     454    }
     455    outrowCount += rowNum ;
    280456    ++iter;
    281457  }
     
    323499      if ( tsys ) {
    324500        ts += val;
     501        tsysCol.put(i, ts);
     502      }
     503    }
     504  }
     505  return out;
     506}
     507
     508CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
     509                                              const std::vector<float> val,
     510                                              const std::string& mode,
     511                                              const std::string& opmode,
     512                                              bool tsys )
     513{
     514  CountedPtr< Scantable > out ;
     515  if ( opmode == "channel" ) {
     516    out = arrayOperateChannel( in, val, mode, tsys ) ;
     517  }
     518  else if ( opmode == "row" ) {
     519    out = arrayOperateRow( in, val, mode, tsys ) ;
     520  }
     521  else {
     522    throw( AipsError( "Unknown array operation mode." ) ) ;
     523  }
     524  return out ;
     525}
     526
     527CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
     528                                                     const std::vector<float> val,
     529                                                     const std::string& mode,
     530                                                     bool tsys )
     531{
     532  if ( val.size() == 1 ){
     533    return unaryOperate( in, val[0], mode, tsys ) ;
     534  }
     535
     536  // conformity of SPECTRA and TSYS
     537  if ( tsys ) {
     538    TableIterator titer(in->table(), "IFNO");
     539    while ( !titer.pastEnd() ) {
     540      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     541      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     542      Array<Float> spec = specCol.getColumn() ;
     543      Array<Float> ts = tsysCol.getColumn() ;
     544      if ( !spec.conform( ts ) ) {
     545        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     546      }
     547      titer.next() ;
     548    }
     549  }
     550
     551  // check if all spectra in the scantable have the same number of channel
     552  vector<uInt> nchans;
     553  vector<uInt> ifnos = in->getIFNos() ;
     554  for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
     555    nchans.push_back( in->nchan( ifnos[i] ) ) ;
     556  }
     557  Vector<uInt> mchans( nchans ) ;
     558  if ( anyNE( mchans, mchans[0] ) ) {
     559    throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
     560  }
     561
     562  // check if vector size is equal to nchan
     563  Vector<Float> fact( val ) ;
     564  if ( fact.nelements() != mchans[0] ) {
     565    throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
     566  }
     567
     568  // check divided by zero
     569  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
     570    throw( AipsError("Divided by zero is not recommended." ) ) ;
     571  }
     572
     573  CountedPtr< Scantable > out = getScantable(in, false);
     574  Table& tab = out->table();
     575  ArrayColumn<Float> specCol(tab,"SPECTRA");
     576  ArrayColumn<Float> tsysCol(tab,"TSYS");
     577  for (uInt i=0; i<tab.nrow(); ++i) {
     578    Vector<Float> spec;
     579    Vector<Float> ts;
     580    specCol.get(i, spec);
     581    tsysCol.get(i, ts);
     582    if (mode == "MUL" || mode == "DIV") {
     583      if (mode == "DIV") fact = (float)1.0 / fact;
     584      spec *= fact;
     585      specCol.put(i, spec);
     586      if ( tsys ) {
     587        ts *= fact;
     588        tsysCol.put(i, ts);
     589      }
     590    } else if ( mode == "ADD"  || mode == "SUB") {
     591      if (mode == "SUB") fact *= (float)-1.0 ;
     592      spec += fact;
     593      specCol.put(i, spec);
     594      if ( tsys ) {
     595        ts += fact;
     596        tsysCol.put(i, ts);
     597      }
     598    }
     599  }
     600  return out;
     601}
     602
     603CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
     604                                                 const std::vector<float> val,
     605                                                 const std::string& mode,
     606                                                 bool tsys )
     607{
     608  if ( val.size() == 1 ) {
     609    return unaryOperate( in, val[0], mode, tsys ) ;
     610  }
     611
     612  // conformity of SPECTRA and TSYS
     613  if ( tsys ) {
     614    TableIterator titer(in->table(), "IFNO");
     615    while ( !titer.pastEnd() ) {
     616      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     617      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     618      Array<Float> spec = specCol.getColumn() ;
     619      Array<Float> ts = tsysCol.getColumn() ;
     620      if ( !spec.conform( ts ) ) {
     621        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     622      }
     623      titer.next() ;
     624    }
     625  }
     626
     627  // check if vector size is equal to nrow
     628  Vector<Float> fact( val ) ;
     629  if ( fact.nelements() != in->nrow() ) {
     630    throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
     631  }
     632
     633  // check divided by zero
     634  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
     635    throw( AipsError("Divided by zero is not recommended." ) ) ;
     636  }
     637
     638  CountedPtr< Scantable > out = getScantable(in, false);
     639  Table& tab = out->table();
     640  ArrayColumn<Float> specCol(tab,"SPECTRA");
     641  ArrayColumn<Float> tsysCol(tab,"TSYS");
     642  if (mode == "DIV") fact = (float)1.0 / fact;
     643  if (mode == "SUB") fact *= (float)-1.0 ;
     644  for (uInt i=0; i<tab.nrow(); ++i) {
     645    Vector<Float> spec;
     646    Vector<Float> ts;
     647    specCol.get(i, spec);
     648    tsysCol.get(i, ts);
     649    if (mode == "MUL" || mode == "DIV") {
     650      spec *= fact[i];
     651      specCol.put(i, spec);
     652      if ( tsys ) {
     653        ts *= fact[i];
     654        tsysCol.put(i, ts);
     655      }
     656    } else if ( mode == "ADD"  || mode == "SUB") {
     657      spec += fact[i];
     658      specCol.put(i, spec);
     659      if ( tsys ) {
     660        ts += fact[i];
     661        tsysCol.put(i, ts);
     662      }
     663    }
     664  }
     665  return out;
     666}
     667
     668CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
     669                                                const std::vector< std::vector<float> > val,
     670                                                const std::string& mode,
     671                                                bool tsys )
     672{
     673  // conformity of SPECTRA and TSYS
     674  if ( tsys ) {
     675    TableIterator titer(in->table(), "IFNO");
     676    while ( !titer.pastEnd() ) {
     677      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
     678      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
     679      Array<Float> spec = specCol.getColumn() ;
     680      Array<Float> ts = tsysCol.getColumn() ;
     681      if ( !spec.conform( ts ) ) {
     682        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
     683      }
     684      titer.next() ;
     685    }
     686  }
     687
     688  // some checks
     689  vector<uInt> nchans;
     690  for ( uInt i = 0 ; i < in->nrow() ; i++ ) {
     691    nchans.push_back( (in->getSpectrum( i )).size() ) ;
     692  }
     693  //Vector<uInt> mchans( nchans ) ;
     694  vector< Vector<Float> > facts ;
     695  for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
     696    Vector<Float> tmp( val[i] ) ;
     697    // check divided by zero
     698    if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
     699      throw( AipsError("Divided by zero is not recommended." ) ) ;
     700    }
     701    // conformity check
     702    if ( tmp.nelements() != nchans[i] ) {
     703      stringstream ss ;
     704      ss << "Row " << i << ": Vector size must be same as number of channel." ;
     705      throw( AipsError( ss.str() ) ) ;
     706    }
     707    facts.push_back( tmp ) ;
     708  }
     709
     710
     711  CountedPtr< Scantable > out = getScantable(in, false);
     712  Table& tab = out->table();
     713  ArrayColumn<Float> specCol(tab,"SPECTRA");
     714  ArrayColumn<Float> tsysCol(tab,"TSYS");
     715  for (uInt i=0; i<tab.nrow(); ++i) {
     716    Vector<Float> fact = facts[i] ;
     717    Vector<Float> spec;
     718    Vector<Float> ts;
     719    specCol.get(i, spec);
     720    tsysCol.get(i, ts);
     721    if (mode == "MUL" || mode == "DIV") {
     722      if (mode == "DIV") fact = (float)1.0 / fact;
     723      spec *= fact;
     724      specCol.put(i, spec);
     725      if ( tsys ) {
     726        ts *= fact;
     727        tsysCol.put(i, ts);
     728      }
     729    } else if ( mode == "ADD"  || mode == "SUB") {
     730      if (mode == "SUB") fact *= (float)-1.0 ;
     731      spec += fact;
     732      specCol.put(i, spec);
     733      if ( tsys ) {
     734        ts += fact;
    325735        tsysCol.put(i, ts);
    326736      }
     
    386796}
    387797
     798MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
     799                                         const Vector<uChar>& f)
     800{
     801  Vector<Bool> mask;
     802  mask.resize(f.shape());
     803  convertArray(mask, f);
     804  return MaskedArray<Double>(s,!mask);
     805}
     806
    388807Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
    389808{
     
    402821  // make this operation non insitu
    403822  const Table& tin = in->table();
    404   Table ons = tin(tin.col("SRCTYPE") == Int(0));
    405   Table offs = tin(tin.col("SRCTYPE") == Int(1));
     823  Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
     824  Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
    406825  if ( offs.nrow() == 0 )
    407826    throw(AipsError("No 'off' scans present."));
     
    6411060         //Debug
    6421061         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
     1062         //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ;
     1063         //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
    6431064         meanoff = sum(spoff)/noff;
    6441065         meandiff = sum(spdiff)/ndiff;
     
    7601181      //Debug
    7611182      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
     1183      //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
     1184      //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
    7621185      // fill the result, replay signal tsys by reference tsys
    7631186      outintCol.put(i, resint);
     
    7821205  setInsitu(false);
    7831206  STSelector sel;
    784   std::vector<int> scan1, scan2, beams;
     1207  std::vector<int> scan1, scan2, beams, types;
    7851208  std::vector< vector<int> > scanpair;
    786   std::vector<string> calstate;
     1209  //std::vector<string> calstate;
     1210  std::vector<int> calstate;
    7871211  String msg;
    7881212
     
    8311255  scanpair.push_back(scan1);
    8321256  scanpair.push_back(scan2);
    833   calstate.push_back("*calon");
    834   calstate.push_back("*[^calon]");
     1257  //calstate.push_back("*calon");
     1258  //calstate.push_back("*[^calon]");
     1259  calstate.push_back(SrcType::NODCAL);
     1260  calstate.push_back(SrcType::NOD);
    8351261  CountedPtr< Scantable > ws = getScantable(s, false);
    8361262  uInt l=0;
     
    8411267          sel.reset();
    8421268          sel.setScans(scanpair[i]);
    843           sel.setName(calstate[k]);
     1269          //sel.setName(calstate[k]);
     1270          types.clear();
     1271          types.push_back(calstate[k]);
     1272          sel.setTypes(types);
    8441273          beams.clear();
    8451274          beams.push_back(j);
     
    9261355      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
    9271356      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
     1357      // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
     1358      // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
    9281359      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
    9291360      Array<Float> avtsys =  tsys1;
     
    9521383  CountedPtr< Scantable > ws = getScantable(s, false);
    9531384  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
    954   CountedPtr< Scantable > calsig, calref, out;
     1385  CountedPtr< Scantable > calsig, calref, out, out1, out2;
     1386  Bool nofold=False;
     1387  vector<int> types ;
    9551388
    9561389  //split the data
    957   sel.setName("*_fs");
     1390  //sel.setName("*_fs");
     1391  types.push_back( SrcType::FSON ) ;
     1392  sel.setTypes( types ) ;
    9581393  ws->setSelection(sel);
    9591394  sig = getScantable(ws,false);
    9601395  sel.reset();
    961   sel.setName("*_fs_calon");
     1396  types.clear() ;
     1397  //sel.setName("*_fs_calon");
     1398  types.push_back( SrcType::FONCAL ) ;
     1399  sel.setTypes( types ) ;
    9621400  ws->setSelection(sel);
    9631401  sigwcal = getScantable(ws,false);
    9641402  sel.reset();
    965   sel.setName("*_fsr");
     1403  types.clear() ;
     1404  //sel.setName("*_fsr");
     1405  types.push_back( SrcType::FSOFF ) ;
     1406  sel.setTypes( types ) ;
    9661407  ws->setSelection(sel);
    9671408  ref = getScantable(ws,false);
    9681409  sel.reset();
    969   sel.setName("*_fsr_calon");
     1410  types.clear() ;
     1411  //sel.setName("*_fsr_calon");
     1412  types.push_back( SrcType::FOFFCAL ) ;
     1413  sel.setTypes( types ) ;
    9701414  ws->setSelection(sel);
    9711415  refwcal = getScantable(ws,false);
     1416  sel.reset() ;
     1417  types.clear() ;
    9721418
    9731419  calsig = dototalpower(sigwcal, sig, tcal=tcal);
    9741420  calref = dototalpower(refwcal, ref, tcal=tcal);
    9751421
    976   out=dosigref(calsig,calref,smoothref,tsysv,tau);
    977 
     1422  out1=dosigref(calsig,calref,smoothref,tsysv,tau);
     1423  out2=dosigref(calref,calsig,smoothref,tsysv,tau);
     1424
     1425  Table& tabout1=out1->table();
     1426  Table& tabout2=out2->table();
     1427  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
     1428  ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
     1429  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
     1430  Vector<Float> spec; specCol.get(0, spec);
     1431  uInt nchan = spec.nelements();
     1432  uInt freqid1; freqidCol1.get(0,freqid1);
     1433  uInt freqid2; freqidCol2.get(0,freqid2);
     1434  Double rp1, rp2, rv1, rv2, inc1, inc2;
     1435  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
     1436  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
     1437  //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
     1438  //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
     1439  //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
     1440  if (rp1==rp2) {
     1441    Double foffset = rv1 - rv2;
     1442    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
     1443    if (choffset >= nchan) {
     1444      //cerr<<"out-band frequency switching, no folding"<<endl;
     1445      LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
     1446      os<<"out-band frequency switching, no folding"<<LogIO::POST;
     1447      nofold = True;
     1448    }
     1449  }
     1450
     1451  if (nofold) {
     1452    std::vector< CountedPtr< Scantable > > tabs;
     1453    tabs.push_back(out1);
     1454    tabs.push_back(out2);
     1455    out = merge(tabs);
     1456  }
     1457  else {
     1458    //out = out1;
     1459    Double choffset = ( rv1 - rv2 ) / inc2 ;
     1460    out = dofold( out1, out2, choffset ) ;
     1461  }
     1462   
    9781463  return out;
     1464}
     1465
     1466CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
     1467                                      const CountedPtr<Scantable> &ref,
     1468                                      Double choffset,
     1469                                      Double choffset2 )
     1470{
     1471  LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
     1472  os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
     1473
     1474  // output scantable
     1475  CountedPtr<Scantable> out = getScantable( sig, false ) ;
     1476
     1477  // separate choffset to integer part and decimal part
     1478  Int ioffset = (Int)choffset ;
     1479  Double doffset = choffset - ioffset ;
     1480  Int ioffset2 = (Int)choffset2 ;
     1481  Double doffset2 = choffset2 - ioffset2 ;
     1482  os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
     1483  os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ; 
     1484
     1485  // get column
     1486  ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
     1487  ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
     1488  ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
     1489  ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
     1490  ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
     1491  ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
     1492  ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
     1493  ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
     1494  ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
     1495  ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
     1496
     1497  // check
     1498  if ( ioffset == 0 ) {
     1499    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
     1500    os << "channel offset is zero, no folding" << LogIO::POST ;
     1501    return out ;
     1502  }
     1503  int nchan = ref->nchan() ;
     1504  if ( abs(ioffset) >= nchan ) {
     1505    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
     1506    os << "out-band frequency switching, no folding" << LogIO::POST ;
     1507    return out ;
     1508  }
     1509
     1510  // attach column for output scantable
     1511  ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
     1512  ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
     1513  ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
     1514  ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
     1515  ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
     1516  ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
     1517
     1518  // for each row
     1519  // assume that the data order are same between sig and ref
     1520  RowAccumulator acc( asap::W_TINTSYS ) ;
     1521  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
     1522    // get values
     1523    Vector<Float> spsig ;
     1524    specCol1.get( i, spsig ) ;
     1525    Vector<Float> spref ;
     1526    specCol2.get( i, spref ) ;
     1527    Vector<Float> tsyssig ;
     1528    tsysCol1.get( i, tsyssig ) ;
     1529    Vector<Float> tsysref ;
     1530    tsysCol2.get( i, tsysref ) ;
     1531    Vector<uChar> flagsig ;
     1532    flagCol1.get( i, flagsig ) ;
     1533    Vector<uChar> flagref ;
     1534    flagCol2.get( i, flagref ) ;
     1535    Double timesig ;
     1536    mjdCol1.get( i, timesig ) ;
     1537    Double timeref ;
     1538    mjdCol2.get( i, timeref ) ;
     1539    Double intsig ;
     1540    intervalCol1.get( i, intsig ) ;
     1541    Double intref ;
     1542    intervalCol2.get( i, intref ) ;
     1543
     1544    // shift reference spectra
     1545    int refchan = spref.nelements() ;
     1546    Vector<Float> sspref( spref.nelements() ) ;
     1547    Vector<Float> stsysref( tsysref.nelements() ) ;
     1548    Vector<uChar> sflagref( flagref.nelements() ) ;
     1549    if ( ioffset > 0 ) {
     1550      // SPECTRA and FLAGTRA
     1551      for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
     1552        sspref[j] = spref[j+ioffset] ;
     1553        sflagref[j] = flagref[j+ioffset] ;
     1554      }
     1555      for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
     1556        sspref[j] = spref[j-refchan+ioffset] ;
     1557        sflagref[j] = flagref[j-refchan+ioffset] ;
     1558      }
     1559      spref = sspref.copy() ;
     1560      flagref = sflagref.copy() ;
     1561      for ( int j = 0 ; j < refchan - 1 ; j++ ) {
     1562        sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
     1563        sflagref[j] = flagref[j+1] + flagref[j] ;
     1564      }
     1565      sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
     1566      sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
     1567
     1568      // TSYS
     1569      if ( spref.nelements() == tsysref.nelements() ) {
     1570        for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
     1571          stsysref[j] = tsysref[j+ioffset] ;
     1572        }
     1573        for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
     1574          stsysref[j] = tsysref[j-refchan+ioffset] ;
     1575        }
     1576        tsysref = stsysref.copy() ;
     1577        for ( int j = 0 ; j < refchan - 1 ; j++ ) {
     1578          stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
     1579        }
     1580        stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
     1581      }
     1582    }
     1583    else {
     1584      // SPECTRA and FLAGTRA
     1585      for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
     1586        sspref[j] = spref[refchan+ioffset+j] ;
     1587        sflagref[j] = flagref[refchan+ioffset+j] ;
     1588      }
     1589      for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
     1590        sspref[j] = spref[j+ioffset] ;
     1591        sflagref[j] = flagref[j+ioffset] ;
     1592      }
     1593      spref = sspref.copy() ;
     1594      flagref = sflagref.copy() ;
     1595      sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
     1596      sflagref[0] = flagref[0] + flagref[refchan-1] ;
     1597      for ( int j = 1 ; j < refchan ; j++ ) {
     1598        sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
     1599        sflagref[j] = flagref[j-1] + flagref[j] ;
     1600      }
     1601      // TSYS
     1602      if ( spref.nelements() == tsysref.nelements() ) {
     1603        for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
     1604          stsysref[j] = tsysref[refchan+ioffset+j] ;
     1605        }
     1606        for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
     1607          stsysref[j] = tsysref[j+ioffset] ;
     1608        }
     1609        tsysref = stsysref.copy() ;
     1610        stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
     1611        for ( int j = 1 ; j < refchan ; j++ ) {
     1612          stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
     1613        }
     1614      }
     1615    }
     1616
     1617    // shift signal spectra if necessary (only for APEX?)
     1618    if ( choffset2 != 0.0 ) {
     1619      int sigchan = spsig.nelements() ;
     1620      Vector<Float> sspsig( spsig.nelements() ) ;
     1621      Vector<Float> stsyssig( tsyssig.nelements() ) ;
     1622      Vector<uChar> sflagsig( flagsig.nelements() ) ;
     1623      if ( ioffset2 > 0 ) {
     1624        // SPECTRA and FLAGTRA
     1625        for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
     1626          sspsig[j] = spsig[j+ioffset2] ;
     1627          sflagsig[j] = flagsig[j+ioffset2] ;
     1628        }
     1629        for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
     1630          sspsig[j] = spsig[j-sigchan+ioffset2] ;
     1631          sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
     1632        }
     1633        spsig = sspsig.copy() ;
     1634        flagsig = sflagsig.copy() ;
     1635        for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
     1636          sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
     1637          sflagsig[j] = flagsig[j+1] || flagsig[j] ;
     1638        }
     1639        sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
     1640        sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ;
     1641        // TSTS
     1642        if ( spsig.nelements() == tsyssig.nelements() ) {
     1643          for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
     1644            stsyssig[j] = tsyssig[j+ioffset2] ;
     1645          }
     1646          for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
     1647            stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
     1648          }
     1649          tsyssig = stsyssig.copy() ;
     1650          for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
     1651            stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
     1652          }
     1653          stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
     1654        }
     1655      }
     1656      else {
     1657        // SPECTRA and FLAGTRA
     1658        for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
     1659          sspsig[j] = spsig[sigchan+ioffset2+j] ;
     1660          sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
     1661        }
     1662        for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
     1663          sspsig[j] = spsig[j+ioffset2] ;
     1664          sflagsig[j] = flagsig[j+ioffset2] ;
     1665        }
     1666        spsig = sspsig.copy() ;
     1667        flagsig = sflagsig.copy() ;
     1668        sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
     1669        sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
     1670        for ( int j = 1 ; j < sigchan ; j++ ) {
     1671          sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
     1672          sflagsig[j] = flagsig[j-1] + flagsig[j] ;
     1673        }
     1674        // TSYS
     1675        if ( spsig.nelements() == tsyssig.nelements() ) {
     1676          for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
     1677            stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
     1678          }
     1679          for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
     1680            stsyssig[j] = tsyssig[j+ioffset2] ;
     1681          }
     1682          tsyssig = stsyssig.copy() ;
     1683          stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
     1684          for ( int j = 1 ; j < sigchan ; j++ ) {
     1685            stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
     1686          }
     1687        }
     1688      }
     1689    }
     1690
     1691    // folding
     1692    acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
     1693    acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
     1694   
     1695    // put result
     1696    specColOut.put( i, acc.getSpectrum() ) ;
     1697    const Vector<Bool> &msk = acc.getMask() ;
     1698    Vector<uChar> flg( msk.shape() ) ;
     1699    convertArray( flg, !msk ) ;
     1700    flagColOut.put( i, flg ) ;
     1701    tsysColOut.put( i, acc.getTsys() ) ;
     1702    intervalColOut.put( i, acc.getInterval() ) ;
     1703    mjdColOut.put( i, acc.getTime() ) ;
     1704    // change FREQ_ID to unshifted IF setting (only for APEX?)
     1705    if ( choffset2 != 0.0 ) {
     1706      uInt freqid = fidColOut( 0 ) ; // assume single-IF data
     1707      double refpix, refval, increment ;
     1708      out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
     1709      refval -= choffset * increment ;
     1710      uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
     1711      Vector<uInt> freqids = fidColOut.getColumn() ;
     1712      for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
     1713        if ( freqids[j] == freqid )
     1714          freqids[j] = newfreqid ;
     1715      }
     1716      fidColOut.putColumn( freqids ) ;
     1717    }
     1718
     1719    acc.reset() ;
     1720  }
     1721
     1722  return out ;
    9791723}
    9801724
     
    10511795    }
    10521796    out.push_back(outstat);
     1797  }
     1798  return out;
     1799}
     1800
     1801std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
     1802                                        const std::vector< bool > & mask,
     1803                                        const std::string& which )
     1804{
     1805
     1806  Vector<Bool> m(mask);
     1807  const Table& tab = in->table();
     1808  ROArrayColumn<Float> specCol(tab, "SPECTRA");
     1809  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     1810  std::vector<int> out;
     1811  for (uInt i=0; i < tab.nrow(); ++i ) {
     1812    Vector<Float> spec; specCol.get(i, spec);
     1813    Vector<uChar> flag; flagCol.get(i, flag);
     1814    MaskedArray<Float> ma  = maskedArray(spec, flag);
     1815    if (ma.ndim() != 1) {
     1816      throw (ArrayError(
     1817          "std::vector<int> STMath::minMaxChan("
     1818          "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
     1819          " std::string &which)"
     1820          " - MaskedArray is not 1D"));
     1821    }
     1822    IPosition outpos(1,0);
     1823    if ( spec.nelements() == m.nelements() ) {
     1824      outpos = mathutil::minMaxPos(which, ma(m));
     1825    } else {
     1826      outpos = mathutil::minMaxPos(which, ma);
     1827    }
     1828    out.push_back(outpos[0]);
    10531829  }
    10541830  return out;
     
    15712347    if ( ! (*it)->conformant(*out) ) {
    15722348      // non conformant.
    1573       pushLog(String("Warning: Can't merge scantables as header info differs."));
     2349      //pushLog(String("Warning: Can't merge scantables as header info differs."));
     2350      LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
     2351      os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
    15742352    }
    15752353    out->appendToHistoryTable((*it)->history());
     
    15932371          id = out->frequencies().addEntry(rp, rv, inc);
    15942372          freqidcol.put(k,id);
    1595           String name,fname;Double rf;
     2373          //String name,fname;Double rf;
     2374          Vector<String> name,fname;Vector<Double> rf;
    15962375          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
    15972376          id = out->molecules().addEntry(rf, name, fname);
     
    19972776      int fstart = -1;
    19982777      int fend = -1;
    1999       for (int k=0; k < flag.nelements(); ++k ) {
     2778      for (unsigned int k=0; k < flag.nelements(); ++k ) {
    20002779        if (flag[k] > 0) {
    20012780          fstart = k;
     
    20512830  return out;
    20522831}
     2832
     2833// Averaging spectra with different channel/resolution
     2834CountedPtr<Scantable>
     2835STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
     2836                     const bool& compel,
     2837                     const std::vector<bool>& mask,
     2838                     const std::string& weight,
     2839                     const std::string& avmode )
     2840  throw ( casa::AipsError )
     2841{
     2842  LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
     2843  if ( avmode == "SCAN" && in.size() != 1 )
     2844    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
     2845                    "Use merge first."));
     2846 
     2847  // check if OTF observation
     2848  String obstype = in[0]->getHeader().obstype ;
     2849  Double tol = 0.0 ;
     2850  if ( obstype.find( "OTF" ) != String::npos ) {
     2851    tol = TOL_OTF ;
     2852  }
     2853  else {
     2854    tol = TOL_POINT ;
     2855  }
     2856
     2857  CountedPtr<Scantable> out ;     // processed result
     2858  if ( compel ) {
     2859    std::vector< CountedPtr<Scantable> > newin ; // input for average process
     2860    uInt insize = in.size() ;    // number of input scantables
     2861
     2862    // TEST: do normal average in each table before IF grouping
     2863    os << "Do preliminary averaging" << LogIO::POST ;
     2864    vector< CountedPtr<Scantable> > tmpin( insize ) ;
     2865    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2866      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
     2867      tmpin[itable] = average( v, mask, weight, avmode ) ;
     2868    }
     2869
     2870    // warning
     2871    os << "Average spectra with different spectral resolution" << LogIO::POST ;
     2872
     2873    // temporarily set coordinfo
     2874    vector<string> oldinfo( insize ) ;
     2875    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2876      vector<string> coordinfo = in[itable]->getCoordInfo() ;
     2877      oldinfo[itable] = coordinfo[0] ;
     2878      coordinfo[0] = "Hz" ;
     2879      tmpin[itable]->setCoordInfo( coordinfo ) ;
     2880    }
     2881
     2882    // columns
     2883    ScalarColumn<uInt> freqIDCol ;
     2884    ScalarColumn<uInt> ifnoCol ;
     2885    ScalarColumn<uInt> scannoCol ;
     2886
     2887
     2888    // check IF frequency coverage
     2889    // freqid: list of FREQ_ID, which is used, in each table 
     2890    // iffreq: list of minimum and maximum frequency for each FREQ_ID in
     2891    //         each table
     2892    // freqid[insize][numIF]
     2893    // freqid: [[id00, id01, ...],
     2894    //          [id10, id11, ...],
     2895    //          ...
     2896    //          [idn0, idn1, ...]]
     2897    // iffreq[insize][numIF*2]
     2898    // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
     2899    //          [min_id10, max_id10, min_id11, max_id11, ...],
     2900    //          ...
     2901    //          [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
     2902    //os << "Check IF settings in each table" << LogIO::POST ;
     2903    vector< vector<uInt> > freqid( insize );
     2904    vector< vector<double> > iffreq( insize ) ;
     2905    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2906      uInt rows = tmpin[itable]->nrow() ;
     2907      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
     2908      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
     2909        if ( freqid[itable].size() == freqnrows ) {
     2910          break ;
     2911        }
     2912        else {
     2913          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
     2914          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
     2915          uInt id = freqIDCol( irow ) ;
     2916          if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
     2917            //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
     2918            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
     2919            freqid[itable].push_back( id ) ;
     2920            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
     2921            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
     2922          }
     2923        }
     2924      }
     2925    }
     2926
     2927    // debug
     2928    //os << "IF settings summary:" << endl ;
     2929    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
     2930    //os << "   Table" << i << endl ;
     2931    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
     2932    //os << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
     2933    //}
     2934    //}
     2935    //os << endl ;
     2936    //os.post() ;
     2937
     2938    // IF grouping based on their frequency coverage
     2939    // ifgrp: list of table index and FREQ_ID for all members in each IF group
     2940    // ifgfreq: list of minimum and maximum frequency in each IF group
     2941    // ifgrp[numgrp][nummember*2]
     2942    // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
     2943    //         [table10, freqrow10, table11, freqrow11, ...],
     2944    //         ...
     2945    //         [tablen0, freqrown0, tablen1, freqrown1, ...]]
     2946    // ifgfreq[numgrp*2]
     2947    // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
     2948    //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
     2949    vector< vector<uInt> > ifgrp ;
     2950    vector<double> ifgfreq ;
     2951
     2952    // parameter for IF grouping
     2953    // groupmode = OR    retrieve all region
     2954    //             AND   only retrieve overlaped region
     2955    //string groupmode = "AND" ;
     2956    string groupmode = "OR" ;
     2957    uInt sizecr = 0 ;
     2958    if ( groupmode == "AND" )
     2959      sizecr = 2 ;
     2960    else if ( groupmode == "OR" )
     2961      sizecr = 0 ;
     2962
     2963    vector<double> sortedfreq ;
     2964    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
     2965      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
     2966        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
     2967          sortedfreq.push_back( iffreq[i][j] ) ;
     2968      }
     2969    }
     2970    sort( sortedfreq.begin(), sortedfreq.end() ) ;
     2971    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
     2972      ifgfreq.push_back( *i ) ;
     2973      ifgfreq.push_back( *(i+1) ) ;
     2974    }
     2975    ifgrp.resize( ifgfreq.size()/2 ) ;
     2976    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     2977      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
     2978        double range0 = iffreq[itable][2*iif] ;
     2979        double range1 = iffreq[itable][2*iif+1] ;
     2980        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
     2981          double fmin = max( range0, ifgfreq[2*j] ) ;
     2982          double fmax = min( range1, ifgfreq[2*j+1] ) ;
     2983          if ( fmin < fmax ) {
     2984            ifgrp[j].push_back( itable ) ;
     2985            ifgrp[j].push_back( freqid[itable][iif] ) ;
     2986          }
     2987        }
     2988      }
     2989    }
     2990    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
     2991    vector<double>::iterator giter = ifgfreq.begin() ;
     2992    while( fiter != ifgrp.end() ) {
     2993      if ( fiter->size() <= sizecr ) {
     2994        fiter = ifgrp.erase( fiter ) ;
     2995        giter = ifgfreq.erase( giter ) ;
     2996        giter = ifgfreq.erase( giter ) ;
     2997      }
     2998      else {
     2999        fiter++ ;
     3000        advance( giter, 2 ) ;
     3001      }
     3002    }
     3003
     3004    // Grouping continuous IF groups (without frequency gap)
     3005    // freqgrp: list of IF group indexes in each frequency group
     3006    // freqrange: list of minimum and maximum frequency in each frequency group
     3007    // freqgrp[numgrp][nummember]
     3008    // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
     3009    //           [ifgrp10, ifgrp11, ifgrp12, ...],
     3010    //           ...
     3011    //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
     3012    // freqrange[numgrp*2]
     3013    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
     3014    vector< vector<uInt> > freqgrp ;
     3015    double freqrange = 0.0 ;
     3016    uInt grpnum = 0 ;
     3017    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     3018      // Assumed that ifgfreq was sorted
     3019      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
     3020        freqgrp[grpnum-1].push_back( i ) ;
     3021      }
     3022      else {
     3023        vector<uInt> grp0( 1, i ) ;
     3024        freqgrp.push_back( grp0 ) ;
     3025        grpnum++ ;
     3026      }
     3027      freqrange = ifgfreq[2*i+1] ;
     3028    }
     3029       
     3030
     3031    // print IF groups
     3032    ostringstream oss ;
     3033    oss << "IF Group summary: " << endl ;
     3034    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
     3035    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     3036      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
     3037      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
     3038        oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
     3039      }
     3040      oss << endl ;
     3041    }
     3042    oss << endl ;
     3043    os << oss.str() << LogIO::POST ;
     3044   
     3045    // print frequency group
     3046    oss.str("") ;
     3047    oss << "Frequency Group summary: " << endl ;
     3048    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
     3049    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
     3050      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
     3051      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
     3052        oss << freqgrp[i][j] << " " ;
     3053      }
     3054      oss << endl ;
     3055    }
     3056    oss << endl ;
     3057    os << oss.str() << LogIO::POST ;
     3058
     3059    // membership check
     3060    // groups: list of IF group indexes whose frequency range overlaps with
     3061    //         that of each table and IF
     3062    // groups[numtable][numIF][nummembership]
     3063    // groups: [[[grp, grp,...], [grp, grp,...],...],
     3064    //          [[grp, grp,...], [grp, grp,...],...],
     3065    //          ...
     3066    //          [[grp, grp,...], [grp, grp,...],...]]
     3067    vector< vector< vector<uInt> > > groups( insize ) ;
     3068    for ( uInt i = 0 ; i < insize ; i++ ) {
     3069      groups[i].resize( freqid[i].size() ) ;
     3070    }
     3071    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
     3072      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     3073        uInt tableid = ifgrp[igrp][2*imem] ;
     3074        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
     3075        if ( iter != freqid[tableid].end() ) {
     3076          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
     3077          groups[tableid][rowid].push_back( igrp ) ;
     3078        }
     3079      }
     3080    }
     3081
     3082    // print membership
     3083    //oss.str("") ;
     3084    //for ( uInt i = 0 ; i < insize ; i++ ) {
     3085    //oss << "Table " << i << endl ;
     3086    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
     3087    //oss << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
     3088    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
     3089    //oss << setw( 2 ) << groups[i][j][k] << " " ;
     3090    //}
     3091    //oss << endl ;
     3092    //}
     3093    //}
     3094    //os << oss.str() << LogIO::POST ;
     3095
     3096    // set back coordinfo
     3097    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3098      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
     3099      coordinfo[0] = oldinfo[itable] ;
     3100      tmpin[itable]->setCoordInfo( coordinfo ) ;
     3101    }
     3102
     3103    // Create additional table if needed
     3104    bool oldInsitu = insitu_ ;
     3105    setInsitu( false ) ;
     3106    vector< vector<uInt> > addrow( insize ) ;
     3107    vector<uInt> addtable( insize, 0 ) ;
     3108    vector<uInt> newtableids( insize ) ;
     3109    vector<uInt> newifids( insize, 0 ) ;
     3110    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3111      //os << "Table " << itable << ": " ;
     3112      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
     3113        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
     3114        //os << addrow[itable][ifrow] << " " ;
     3115      }
     3116      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
     3117      //os << "(" << addtable[itable] << ")" << LogIO::POST ;
     3118    }
     3119    newin.resize( insize ) ;
     3120    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
     3121    for ( uInt i = 0 ; i < insize ; i++ ) {
     3122      newtableids[i] = i ;
     3123    }
     3124    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3125      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
     3126        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
     3127        vector<int> freqidlist ;
     3128        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
     3129          if ( groups[itable][i].size() > iadd + 1 ) {
     3130            freqidlist.push_back( freqid[itable][i] ) ;
     3131          }
     3132        }
     3133        stringstream taqlstream ;
     3134        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
     3135        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
     3136          taqlstream << i ;
     3137          if ( i < freqidlist.size() - 1 )
     3138            taqlstream << "," ;
     3139          else
     3140            taqlstream << "]" ;
     3141        }
     3142        string taql = taqlstream.str() ;
     3143        //os << "taql = " << taql << LogIO::POST ;
     3144        STSelector selector = STSelector() ;
     3145        selector.setTaQL( taql ) ;
     3146        add->setSelection( selector ) ;
     3147        newin.push_back( add ) ;
     3148        newtableids.push_back( itable ) ;
     3149        newifids.push_back( iadd + 1 ) ;
     3150      }
     3151    }
     3152
     3153    // udpate ifgrp
     3154    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3155      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
     3156        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
     3157          if ( groups[itable][ifrow].size() > iadd + 1 ) {
     3158            uInt igrp = groups[itable][ifrow][iadd+1] ;
     3159            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     3160              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
     3161                ifgrp[igrp][2*imem] = insize + iadd ;
     3162              }
     3163            }
     3164          }
     3165        }
     3166      }
     3167    }
     3168
     3169    // print IF groups again for debug
     3170    //oss.str( "" ) ;
     3171    //oss << "IF Group summary: " << endl ;
     3172    //oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
     3173    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
     3174    //oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
     3175    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
     3176    //oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
     3177    //}
     3178    //oss << endl ;
     3179    //}
     3180    //oss << endl ;
     3181    //os << oss.str() << LogIO::POST ;
     3182
     3183    // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
     3184    os << "All scan number is set to 0" << LogIO::POST ;
     3185    //os << "All IF number is set to IF group index" << LogIO::POST ;
     3186    insize = newin.size() ;
     3187    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3188      uInt rows = newin[itable]->nrow() ;
     3189      Table &tmpt = newin[itable]->table() ;
     3190      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
     3191      scannoCol.attach( tmpt, "SCANNO" ) ;
     3192      ifnoCol.attach( tmpt, "IFNO" ) ;
     3193      for ( uInt irow=0 ; irow < rows ; irow++ ) {
     3194        scannoCol.put( irow, 0 ) ;
     3195        uInt freqID = freqIDCol( irow ) ;
     3196        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
     3197        if ( iter != freqid[newtableids[itable]].end() ) {
     3198          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
     3199          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
     3200        }
     3201        else {
     3202          throw(AipsError("IF grouping was wrong in additional tables.")) ;
     3203        }
     3204      }
     3205    }
     3206    oldinfo.resize( insize ) ;
     3207    setInsitu( oldInsitu ) ;
     3208
     3209    // temporarily set coordinfo
     3210    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3211      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
     3212      oldinfo[itable] = coordinfo[0] ;
     3213      coordinfo[0] = "Hz" ;
     3214      newin[itable]->setCoordInfo( coordinfo ) ;
     3215    }
     3216
     3217    // save column values in the vector
     3218    vector< vector<uInt> > freqTableIdVec( insize ) ;
     3219    vector< vector<uInt> > freqIdVec( insize ) ;
     3220    vector< vector<uInt> > ifNoVec( insize ) ;
     3221    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3222      ScalarColumn<uInt> freqIDs ;
     3223      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
     3224      ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
     3225      freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
     3226      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
     3227        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
     3228      }
     3229      for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
     3230        freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
     3231        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
     3232      }
     3233    }
     3234
     3235    // reset spectra and flagtra: pick up common part of frequency coverage
     3236    //os << "Pick common frequency range and align resolution" << LogIO::POST ;
     3237    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3238      uInt rows = newin[itable]->nrow() ;
     3239      int nminchan = -1 ;
     3240      int nmaxchan = -1 ;
     3241      vector<uInt> freqIdUpdate ;
     3242      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
     3243        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
     3244        double minfreq = ifgfreq[2*ifno] ;
     3245        double maxfreq = ifgfreq[2*ifno+1] ;
     3246        //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ;
     3247        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
     3248        int nchan = abcissa.size() ;
     3249        double resol = abcissa[1] - abcissa[0] ;
     3250        //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
     3251        if ( minfreq <= abcissa[0] )
     3252          nminchan = 0 ;
     3253        else {
     3254          //double cfreq = ( minfreq - abcissa[0] ) / resol ;
     3255          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
     3256          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
     3257        }
     3258        if ( maxfreq >= abcissa[abcissa.size()-1] )
     3259          nmaxchan = abcissa.size() - 1 ;
     3260        else {
     3261          //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
     3262          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
     3263          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
     3264        }
     3265        //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
     3266        if ( nmaxchan > nminchan ) {
     3267          newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
     3268          int newchan = nmaxchan - nminchan + 1 ;
     3269          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
     3270            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
     3271        }
     3272        else {
     3273          throw(AipsError("Failed to pick up common part of frequency range.")) ;
     3274        }
     3275      }
     3276      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
     3277        uInt freqId = freqIdUpdate[i] ;
     3278        Double refpix ;
     3279        Double refval ;
     3280        Double increment ;
     3281       
     3282        // update row
     3283        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
     3284        refval = refval - ( refpix - nminchan ) * increment ;
     3285        refpix = 0 ;
     3286        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
     3287      }   
     3288    }
     3289
     3290   
     3291    // reset spectra and flagtra: align spectral resolution
     3292    //os << "Align spectral resolution" << LogIO::POST ;
     3293    // gmaxdnu: the coarsest frequency resolution in the frequency group
     3294    // gmemid: member index that have a resolution equal to gmaxdnu
     3295    // gmaxdnu[numfreqgrp]
     3296    // gmaxdnu: [dnu0, dnu1, ...]
     3297    // gmemid[numfreqgrp]
     3298    // gmemid: [id0, id1, ...]
     3299    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
     3300    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
     3301    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
     3302      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution
     3303      int minchan = INT_MAX ;     // minimum channel number
     3304      Double refpixref = -1 ;     // reference of 'reference pixel'
     3305      Double refvalref = -1 ;     // reference of 'reference frequency'
     3306      Double refinc = -1 ;        // reference frequency resolution
     3307      uInt refreqid ;
     3308      uInt reftable = INT_MAX;
     3309      // process only if group member > 1
     3310      if ( ifgrp[igrp].size() > 2 ) {
     3311        // find minchan and maxdnu in each group
     3312        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     3313          uInt tableid = ifgrp[igrp][2*imem] ;
     3314          uInt rowid = ifgrp[igrp][2*imem+1] ;
     3315          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
     3316          if ( iter != freqIdVec[tableid].end() ) {
     3317            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
     3318            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
     3319            int nchan = abcissa.size() ;
     3320            double dnu = abcissa[1] - abcissa[0] ;
     3321            //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ;
     3322            if ( nchan < minchan ) {
     3323              minchan = nchan ;
     3324              maxdnu = dnu ;
     3325              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
     3326              refreqid = rowid ;
     3327              reftable = tableid ;
     3328            }
     3329          }
     3330        }
     3331        // regrid spectra in each group
     3332        os << "GROUP " << igrp << endl ;
     3333        os << "   Channel number is adjusted to " << minchan << endl ;
     3334        os << "   Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ;
     3335        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
     3336          uInt tableid = ifgrp[igrp][2*imem] ;
     3337          uInt rowid = ifgrp[igrp][2*imem+1] ;
     3338          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
     3339          //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ;
     3340          //os << "   regridChannel applied to " ;
     3341          if ( tableid != reftable )
     3342            refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
     3343          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
     3344            uInt tfreqid = freqIdVec[tableid][irow] ;
     3345            if ( tfreqid == rowid ) {     
     3346              //os << irow << " " ;
     3347              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
     3348              freqIDCol.put( irow, refreqid ) ;
     3349              freqIdVec[tableid][irow] = refreqid ;
     3350            }
     3351          }
     3352          //os << LogIO::POST ;
     3353        }
     3354      }
     3355      else {
     3356        uInt tableid = ifgrp[igrp][0] ;
     3357        uInt rowid = ifgrp[igrp][1] ;
     3358        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
     3359        if ( iter != freqIdVec[tableid].end() ) {
     3360          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
     3361          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
     3362          minchan = abcissa.size() ;
     3363          maxdnu = abcissa[1] - abcissa[0] ;
     3364        }
     3365      }
     3366      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
     3367        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
     3368          if ( maxdnu > gmaxdnu[i] ) {
     3369            gmaxdnu[i] = maxdnu ;
     3370            gmemid[i] = igrp ;
     3371          }
     3372          break ;
     3373        }
     3374      }
     3375    }
     3376
     3377    // set back coordinfo
     3378    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3379      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
     3380      coordinfo[0] = oldinfo[itable] ;
     3381      newin[itable]->setCoordInfo( coordinfo ) ;
     3382    }     
     3383
     3384    // accumulate all rows into the first table
     3385    // NOTE: assumed in.size() = 1
     3386    vector< CountedPtr<Scantable> > tmp( 1 ) ;
     3387    if ( newin.size() == 1 )
     3388      tmp[0] = newin[0] ;
     3389    else
     3390      tmp[0] = merge( newin ) ;
     3391
     3392    //return tmp[0] ;
     3393
     3394    // average
     3395    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
     3396
     3397    //return tmpout ;
     3398
     3399    // combine frequency group
     3400    os << "Combine spectra based on frequency grouping" << LogIO::POST ;
     3401    os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
     3402    vector<string> coordinfo = tmpout->getCoordInfo() ;
     3403    oldinfo[0] = coordinfo[0] ;
     3404    coordinfo[0] = "Hz" ;
     3405    tmpout->setCoordInfo( coordinfo ) ;
     3406    // create proformas of output table
     3407    stringstream taqlstream ;
     3408    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
     3409    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
     3410      taqlstream << gmemid[i] ;
     3411      if ( i < gmemid.size() - 1 )
     3412        taqlstream << "," ;
     3413      else
     3414        taqlstream << "]" ;
     3415    }
     3416    string taql = taqlstream.str() ;
     3417    //os << "taql = " << taql << LogIO::POST ;
     3418    STSelector selector = STSelector() ;
     3419    selector.setTaQL( taql ) ;
     3420    oldInsitu = insitu_ ;
     3421    setInsitu( false ) ;
     3422    out = getScantable( tmpout, false ) ;
     3423    setInsitu( oldInsitu ) ;
     3424    out->setSelection( selector ) ;
     3425    // regrid rows
     3426    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
     3427    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
     3428      uInt ifno = ifnoCol( irow ) ;
     3429      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3430        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
     3431          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
     3432          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
     3433          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
     3434          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
     3435          break ;
     3436        }
     3437      }
     3438    }
     3439    // combine spectra
     3440    ArrayColumn<Float> specColOut ;
     3441    specColOut.attach( out->table(), "SPECTRA" ) ;
     3442    ArrayColumn<uChar> flagColOut ;
     3443    flagColOut.attach( out->table(), "FLAGTRA" ) ;
     3444    ScalarColumn<uInt> ifnoColOut ;
     3445    ifnoColOut.attach( out->table(), "IFNO" ) ;
     3446    ScalarColumn<uInt> polnoColOut ;
     3447    polnoColOut.attach( out->table(), "POLNO" ) ;
     3448    ScalarColumn<uInt> freqidColOut ;
     3449    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
     3450    MDirection::ScalarColumn dirColOut ;
     3451    dirColOut.attach( out->table(), "DIRECTION" ) ;
     3452    Table &tab = tmpout->table() ;
     3453    Block<String> cols(1);
     3454    cols[0] = String("POLNO") ;
     3455    TableIterator iter( tab, cols ) ;
     3456    bool done = false ;
     3457    vector< vector<uInt> > sizes( freqgrp.size() ) ;
     3458    while( !iter.pastEnd() ) {
     3459      vector< vector<Float> > specout( freqgrp.size() ) ;
     3460      vector< vector<uChar> > flagout( freqgrp.size() ) ;
     3461      ArrayColumn<Float> specCols ;
     3462      specCols.attach( iter.table(), "SPECTRA" ) ;
     3463      ArrayColumn<uChar> flagCols ;
     3464      flagCols.attach( iter.table(), "FLAGTRA" ) ;
     3465      ifnoCol.attach( iter.table(), "IFNO" ) ;
     3466      ScalarColumn<uInt> polnos ;
     3467      polnos.attach( iter.table(), "POLNO" ) ;
     3468      MDirection::ScalarColumn dircol ;
     3469      dircol.attach( iter.table(), "DIRECTION" ) ;
     3470      uInt polno = polnos( 0 ) ;
     3471      //os << "POLNO iteration: " << polno << LogIO::POST ;
     3472//       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3473//      sizes[igrp].resize( freqgrp[igrp].size() ) ;
     3474//      for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     3475//        for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
     3476//          uInt ifno = ifnoCol( irow ) ;
     3477//          if ( ifno == freqgrp[igrp][imem] ) {
     3478//            Vector<Float> spec = specCols( irow ) ;
     3479//            Vector<uChar> flag = flagCols( irow ) ;
     3480//            vector<Float> svec ;
     3481//            spec.tovector( svec ) ;
     3482//            vector<uChar> fvec ;
     3483//            flag.tovector( fvec ) ;
     3484//            //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
     3485//            specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
     3486//            flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
     3487//            //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
     3488//            sizes[igrp][imem] = spec.nelements() ;
     3489//          }
     3490//        }
     3491//      }
     3492//      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     3493//        uInt ifout = ifnoColOut( irow ) ;
     3494//        uInt polout = polnoColOut( irow ) ;
     3495//        if ( ifout == gmemid[igrp] && polout == polno ) {
     3496//          // set SPECTRA and FRAGTRA
     3497//          Vector<Float> newspec( specout[igrp] ) ;
     3498//          Vector<uChar> newflag( flagout[igrp] ) ;
     3499//          specColOut.put( irow, newspec ) ;
     3500//          flagColOut.put( irow, newflag ) ;
     3501//          // IFNO renumbering
     3502//          ifnoColOut.put( irow, igrp ) ;
     3503//        }
     3504//      }
     3505//       }
     3506      // get a list of number of channels for each frequency group member
     3507      if ( !done ) {
     3508        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3509          sizes[igrp].resize( freqgrp[igrp].size() ) ;
     3510          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     3511            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
     3512              uInt ifno = ifnoCol( irow ) ;
     3513              if ( ifno == freqgrp[igrp][imem] ) {
     3514                Vector<Float> spec = specCols( irow ) ;
     3515                sizes[igrp][imem] = spec.nelements() ;
     3516                break ;
     3517              }               
     3518            }
     3519          }
     3520        }
     3521        done = true ;
     3522      }
     3523      // combine spectra
     3524      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     3525        uInt polout = polnoColOut( irow ) ;
     3526        if ( polout == polno ) {
     3527          uInt ifout = ifnoColOut( irow ) ;
     3528          Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
     3529          uInt igrp ;
     3530          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
     3531            if ( ifout == gmemid[jgrp] ) {
     3532              igrp = jgrp ;
     3533              break ;
     3534            }
     3535          }
     3536          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     3537            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
     3538              uInt ifno = ifnoCol( jrow ) ;
     3539              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     3540              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
     3541              Double dx = tdir[0] - direction[0] ;
     3542              Double dy = tdir[1] - direction[1] ;
     3543              Double dd = sqrt( dx * dx + dy * dy ) ;
     3544              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
     3545              if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
     3546                Vector<Float> spec = specCols( jrow ) ;
     3547                Vector<uChar> flag = flagCols( jrow ) ;
     3548                vector<Float> svec ;
     3549                spec.tovector( svec ) ;
     3550                vector<uChar> fvec ;
     3551                flag.tovector( fvec ) ;
     3552                //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
     3553                specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
     3554                flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
     3555                //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
     3556              }
     3557            }
     3558          }
     3559          // set SPECTRA and FRAGTRA
     3560          Vector<Float> newspec( specout[igrp] ) ;
     3561          Vector<uChar> newflag( flagout[igrp] ) ;
     3562          specColOut.put( irow, newspec ) ;
     3563          flagColOut.put( irow, newflag ) ;
     3564          // IFNO renumbering
     3565          ifnoColOut.put( irow, igrp ) ;
     3566        }
     3567      }
     3568      iter++ ;
     3569    }
     3570    // update FREQUENCIES subtable
     3571    vector<bool> updated( freqgrp.size(), false ) ;
     3572    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3573      uInt index = 0 ;
     3574      uInt pixShift = 0 ;
     3575      while ( freqgrp[igrp][index] != gmemid[igrp] ) {
     3576        pixShift += sizes[igrp][index++] ;
     3577      }
     3578      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     3579        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
     3580          uInt freqidOut = freqidColOut( irow ) ;
     3581          //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
     3582          double refpix ;
     3583          double refval ;
     3584          double increm ;
     3585          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
     3586          refpix += pixShift ;
     3587          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
     3588          updated[igrp] = true ;
     3589        }
     3590      }
     3591    }
     3592
     3593    //out = tmpout ;
     3594
     3595    coordinfo = tmpout->getCoordInfo() ;
     3596    coordinfo[0] = oldinfo[0] ;
     3597    tmpout->setCoordInfo( coordinfo ) ;
     3598  }
     3599  else {
     3600    // simple average
     3601    out =  average( in, mask, weight, avmode ) ;
     3602  }
     3603 
     3604  return out ;
     3605}
     3606
     3607CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
     3608                                     const String calmode,
     3609                                     const String antname )
     3610{
     3611  // frequency switch
     3612  if ( calmode == "fs" ) {
     3613    return cwcalfs( s, antname ) ;
     3614  }
     3615  else {
     3616    vector<bool> masks = s->getMask( 0 ) ;
     3617    vector<int> types ;
     3618
     3619    // sky scan
     3620    STSelector sel = STSelector() ;
     3621    types.push_back( SrcType::SKY ) ;
     3622    sel.setTypes( types ) ;
     3623    s->setSelection( sel ) ;
     3624    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3625    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
     3626    s->unsetSelection() ;
     3627    sel.reset() ;
     3628    types.clear() ;
     3629
     3630    // hot scan
     3631    types.push_back( SrcType::HOT ) ;
     3632    sel.setTypes( types ) ;
     3633    s->setSelection( sel ) ;
     3634    tmp.clear() ;
     3635    tmp.push_back( getScantable( s, false ) ) ;
     3636    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
     3637    s->unsetSelection() ;
     3638    sel.reset() ;
     3639    types.clear() ;
     3640   
     3641    // cold scan
     3642    CountedPtr<Scantable> acold ;
     3643//     types.push_back( SrcType::COLD ) ;
     3644//     sel.setTypes( types ) ;
     3645//     s->setSelection( sel ) ;
     3646//     tmp.clear() ;
     3647//     tmp.push_back( getScantable( s, false ) ) ;
     3648//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
     3649//     s->unsetSelection() ;
     3650//     sel.reset() ;
     3651//     types.clear() ;
     3652
     3653    // off scan
     3654    types.push_back( SrcType::PSOFF ) ;
     3655    sel.setTypes( types ) ;
     3656    s->setSelection( sel ) ;
     3657    tmp.clear() ;
     3658    tmp.push_back( getScantable( s, false ) ) ;
     3659    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
     3660    s->unsetSelection() ;
     3661    sel.reset() ;
     3662    types.clear() ;
     3663   
     3664    // on scan
     3665    bool insitu = insitu_ ;
     3666    insitu_ = false ;
     3667    CountedPtr<Scantable> out = getScantable( s, true ) ;
     3668    insitu_ = insitu ;
     3669    types.push_back( SrcType::PSON ) ;
     3670    sel.setTypes( types ) ;
     3671    s->setSelection( sel ) ;
     3672    TableCopy::copyRows( out->table(), s->table() ) ;
     3673    s->unsetSelection() ;
     3674    sel.reset() ;
     3675    types.clear() ;
     3676   
     3677    // process each on scan
     3678    ArrayColumn<Float> tsysCol ;
     3679    tsysCol.attach( out->table(), "TSYS" ) ;
     3680    for ( int i = 0 ; i < out->nrow() ; i++ ) {
     3681      vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
     3682      out->setSpectrum( sp, i ) ;
     3683      string reftime = out->getTime( i ) ;
     3684      vector<int> ii( 1, out->getIF( i ) ) ;
     3685      vector<int> ib( 1, out->getBeam( i ) ) ;
     3686      vector<int> ip( 1, out->getPol( i ) ) ;
     3687      sel.setIFs( ii ) ;
     3688      sel.setBeams( ib ) ;
     3689      sel.setPolarizations( ip ) ;
     3690      asky->setSelection( sel ) ;   
     3691      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
     3692      const Vector<Float> Vtsys( sptsys ) ;
     3693      tsysCol.put( i, Vtsys ) ;
     3694      asky->unsetSelection() ;
     3695      sel.reset() ;
     3696    }
     3697
     3698    // flux unit
     3699    out->setFluxUnit( "K" ) ;
     3700
     3701    return out ;
     3702  }
     3703}
     3704 
     3705CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
     3706                                       const String calmode )
     3707{
     3708  // frequency switch
     3709  if ( calmode == "fs" ) {
     3710    return almacalfs( s ) ;
     3711  }
     3712  else {
     3713    vector<bool> masks = s->getMask( 0 ) ;
     3714   
     3715    // off scan
     3716    STSelector sel = STSelector() ;
     3717    vector<int> types ;
     3718    types.push_back( SrcType::PSOFF ) ;
     3719    sel.setTypes( types ) ;
     3720    s->setSelection( sel ) ;
     3721    // TODO 2010/01/08 TN
     3722    // Grouping by time should be needed before averaging.
     3723    // Each group must have own unique SCANNO (should be renumbered).
     3724    // See PIPELINE/SDCalibration.py
     3725    CountedPtr<Scantable> soff = getScantable( s, false ) ;
     3726    Table ttab = soff->table() ;
     3727    ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
     3728    uInt nrow = timeCol.nrow() ;
     3729    Vector<Double> timeSep( nrow - 1 ) ;
     3730    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
     3731      timeSep[i] = timeCol(i+1) - timeCol(i) ;
     3732    }
     3733    ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
     3734    Vector<Double> interval = intervalCol.getColumn() ;
     3735    interval /= 86400.0 ;
     3736    ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
     3737    vector<uInt> glist ;
     3738    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
     3739      double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
     3740      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
     3741      if ( gap > 1.1 ) {
     3742        glist.push_back( i ) ;
     3743      }
     3744    }
     3745    Vector<uInt> gaplist( glist ) ;
     3746    //cout << "gaplist = " << gaplist << endl ;
     3747    uInt newid = 0 ;
     3748    for ( uInt i = 0 ; i < nrow ; i++ ) {
     3749      scanCol.put( i, newid ) ;
     3750      if ( i == gaplist[newid] ) {
     3751        newid++ ;
     3752      }
     3753    }
     3754    //cout << "new scancol = " << scanCol.getColumn() << endl ;
     3755    vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
     3756    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
     3757    //cout << "aoff.nrow = " << aoff->nrow() << endl ;
     3758    s->unsetSelection() ;
     3759    sel.reset() ;
     3760    types.clear() ;
     3761   
     3762    // on scan
     3763    bool insitu = insitu_ ;
     3764    insitu_ = false ;
     3765    CountedPtr<Scantable> out = getScantable( s, true ) ;
     3766    insitu_ = insitu ;
     3767    types.push_back( SrcType::PSON ) ;
     3768    sel.setTypes( types ) ;
     3769    s->setSelection( sel ) ;
     3770    TableCopy::copyRows( out->table(), s->table() ) ;
     3771    s->unsetSelection() ;
     3772    sel.reset() ;
     3773    types.clear() ;
     3774   
     3775    // process each on scan
     3776    ArrayColumn<Float> tsysCol ;
     3777    tsysCol.attach( out->table(), "TSYS" ) ;
     3778    for ( int i = 0 ; i < out->nrow() ; i++ ) {
     3779      vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
     3780      out->setSpectrum( sp, i ) ;
     3781    }
     3782
     3783    // flux unit
     3784    out->setFluxUnit( "K" ) ;
     3785
     3786    return out ;
     3787  }
     3788}
     3789
     3790CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
     3791                                       const String antname )
     3792{
     3793  vector<int> types ;
     3794
     3795  // APEX calibration mode
     3796  int apexcalmode = 1 ;
     3797 
     3798  if ( antname.find( "APEX" ) != string::npos ) {
     3799    // check if off scan exists or not
     3800    STSelector sel = STSelector() ;
     3801    //sel.setName( offstr1 ) ;
     3802    types.push_back( SrcType::FLOOFF ) ;
     3803    sel.setTypes( types ) ;
     3804    try {
     3805      s->setSelection( sel ) ;
     3806    }
     3807    catch ( AipsError &e ) {
     3808      apexcalmode = 0 ;
     3809    }
     3810    sel.reset() ;
     3811  }
     3812  s->unsetSelection() ;
     3813  types.clear() ;
     3814
     3815  vector<bool> masks = s->getMask( 0 ) ;
     3816  CountedPtr<Scantable> ssig, sref ;
     3817  CountedPtr<Scantable> out ;
     3818
     3819  if ( antname.find( "APEX" ) != string::npos ) {
     3820    // APEX calibration
     3821    // sky scan
     3822    STSelector sel = STSelector() ;
     3823    types.push_back( SrcType::FLOSKY ) ;
     3824    sel.setTypes( types ) ;
     3825    s->setSelection( sel ) ;
     3826    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     3827    CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
     3828    s->unsetSelection() ;
     3829    sel.reset() ;
     3830    types.clear() ;
     3831    types.push_back( SrcType::FHISKY ) ;
     3832    sel.setTypes( types ) ;
     3833    s->setSelection( sel ) ;
     3834    tmp.clear() ;
     3835    tmp.push_back( getScantable( s, false ) ) ;
     3836    CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3837    s->unsetSelection() ;
     3838    sel.reset() ;
     3839    types.clear() ;
     3840   
     3841    // hot scan
     3842    types.push_back( SrcType::FLOHOT ) ;
     3843    sel.setTypes( types ) ;
     3844    s->setSelection( sel ) ;
     3845    tmp.clear() ;
     3846    tmp.push_back( getScantable( s, false ) ) ;
     3847    CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3848    s->unsetSelection() ;
     3849    sel.reset() ;
     3850    types.clear() ;
     3851    types.push_back( SrcType::FHIHOT ) ;
     3852    sel.setTypes( types ) ;
     3853    s->setSelection( sel ) ;
     3854    tmp.clear() ;
     3855    tmp.push_back( getScantable( s, false ) ) ;
     3856    CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
     3857    s->unsetSelection() ;
     3858    sel.reset() ;
     3859    types.clear() ;
     3860   
     3861    // cold scan
     3862    CountedPtr<Scantable> acoldlo, acoldhi ;
     3863//     types.push_back( SrcType::FLOCOLD ) ;
     3864//     sel.setTypes( types ) ;
     3865//     s->setSelection( sel ) ;
     3866//     tmp.clear() ;
     3867//     tmp.push_back( getScantable( s, false ) ) ;
     3868//     CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
     3869//     s->unsetSelection() ;
     3870//     sel.reset() ;
     3871//     types.clear() ;
     3872//     types.push_back( SrcType::FHICOLD ) ;
     3873//     sel.setTypes( types ) ;
     3874//     s->setSelection( sel ) ;
     3875//     tmp.clear() ;
     3876//     tmp.push_back( getScantable( s, false ) ) ;
     3877//     CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3878//     s->unsetSelection() ;
     3879//     sel.reset() ;
     3880//     types.clear() ;
     3881
     3882    // ref scan
     3883    bool insitu = insitu_ ;
     3884    insitu_ = false ;
     3885    sref = getScantable( s, true ) ;
     3886    insitu_ = insitu ;
     3887    types.push_back( SrcType::FSLO ) ;
     3888    sel.setTypes( types ) ;
     3889    s->setSelection( sel ) ;
     3890    TableCopy::copyRows( sref->table(), s->table() ) ;
     3891    s->unsetSelection() ;
     3892    sel.reset() ;
     3893    types.clear() ;
     3894   
     3895    // sig scan
     3896    insitu_ = false ;
     3897    ssig = getScantable( s, true ) ;
     3898    insitu_ = insitu ;
     3899    types.push_back( SrcType::FSHI ) ;
     3900    sel.setTypes( types ) ;
     3901    s->setSelection( sel ) ;
     3902    TableCopy::copyRows( ssig->table(), s->table() ) ;
     3903    s->unsetSelection() ;
     3904    sel.reset() ; 
     3905    types.clear() ;
     3906         
     3907    if ( apexcalmode == 0 ) {
     3908      // APEX fs data without off scan
     3909      // process each sig and ref scan
     3910      ArrayColumn<Float> tsysCollo ;
     3911      tsysCollo.attach( ssig->table(), "TSYS" ) ;
     3912      ArrayColumn<Float> tsysColhi ;
     3913      tsysColhi.attach( sref->table(), "TSYS" ) ;
     3914      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     3915        vector< CountedPtr<Scantable> > sky( 2 ) ;
     3916        sky[0] = askylo ;
     3917        sky[1] = askyhi ;
     3918        vector< CountedPtr<Scantable> > hot( 2 ) ;
     3919        hot[0] = ahotlo ;
     3920        hot[1] = ahothi ;
     3921        vector< CountedPtr<Scantable> > cold( 2 ) ;
     3922        //cold[0] = acoldlo ;
     3923        //cold[1] = acoldhi ;
     3924        vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
     3925        ssig->setSpectrum( sp, i ) ;
     3926        string reftime = ssig->getTime( i ) ;
     3927        vector<int> ii( 1, ssig->getIF( i ) ) ;
     3928        vector<int> ib( 1, ssig->getBeam( i ) ) ;
     3929        vector<int> ip( 1, ssig->getPol( i ) ) ;
     3930        sel.setIFs( ii ) ;
     3931        sel.setBeams( ib ) ;
     3932        sel.setPolarizations( ip ) ;
     3933        askylo->setSelection( sel ) ;
     3934        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
     3935        const Vector<Float> Vtsyslo( sptsys ) ;
     3936        tsysCollo.put( i, Vtsyslo ) ;
     3937        askylo->unsetSelection() ;
     3938        sel.reset() ;
     3939        sky[0] = askyhi ;
     3940        sky[1] = askylo ;
     3941        hot[0] = ahothi ;
     3942        hot[1] = ahotlo ;
     3943        cold[0] = acoldhi ;
     3944        cold[1] = acoldlo ;
     3945        sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
     3946        sref->setSpectrum( sp, i ) ;
     3947        reftime = sref->getTime( i ) ;
     3948        ii[0] = sref->getIF( i )  ;
     3949        ib[0] = sref->getBeam( i ) ;
     3950        ip[0] = sref->getPol( i ) ;
     3951        sel.setIFs( ii ) ;
     3952        sel.setBeams( ib ) ;
     3953        sel.setPolarizations( ip ) ;
     3954        askyhi->setSelection( sel ) ;   
     3955        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
     3956        const Vector<Float> Vtsyshi( sptsys ) ;
     3957        tsysColhi.put( i, Vtsyshi ) ;
     3958        askyhi->unsetSelection() ;
     3959        sel.reset() ;
     3960      }
     3961    }
     3962    else if ( apexcalmode == 1 ) {
     3963      // APEX fs data with off scan
     3964      // off scan
     3965      types.push_back( SrcType::FLOOFF ) ;
     3966      sel.setTypes( types ) ;
     3967      s->setSelection( sel ) ;
     3968      tmp.clear() ;
     3969      tmp.push_back( getScantable( s, false ) ) ;
     3970      CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
     3971      s->unsetSelection() ;
     3972      sel.reset() ;
     3973      types.clear() ;
     3974      types.push_back( SrcType::FHIOFF ) ;
     3975      sel.setTypes( types ) ;
     3976      s->setSelection( sel ) ;
     3977      tmp.clear() ;
     3978      tmp.push_back( getScantable( s, false ) ) ;
     3979      CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
     3980      s->unsetSelection() ;
     3981      sel.reset() ;
     3982      types.clear() ;
     3983     
     3984      // process each sig and ref scan
     3985      ArrayColumn<Float> tsysCollo ;
     3986      tsysCollo.attach( ssig->table(), "TSYS" ) ;
     3987      ArrayColumn<Float> tsysColhi ;
     3988      tsysColhi.attach( sref->table(), "TSYS" ) ;
     3989      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     3990        vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
     3991        ssig->setSpectrum( sp, i ) ;
     3992        sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
     3993        string reftime = ssig->getTime( i ) ;
     3994        vector<int> ii( 1, ssig->getIF( i ) ) ;
     3995        vector<int> ib( 1, ssig->getBeam( i ) ) ;
     3996        vector<int> ip( 1, ssig->getPol( i ) ) ;
     3997        sel.setIFs( ii ) ;
     3998        sel.setBeams( ib ) ;
     3999        sel.setPolarizations( ip ) ;
     4000        askylo->setSelection( sel ) ;
     4001        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
     4002        const Vector<Float> Vtsyslo( sptsys ) ;
     4003        tsysCollo.put( i, Vtsyslo ) ;
     4004        askylo->unsetSelection() ;
     4005        sel.reset() ;
     4006        sref->setSpectrum( sp, i ) ;
     4007        reftime = sref->getTime( i ) ;
     4008        ii[0] = sref->getIF( i )  ;
     4009        ib[0] = sref->getBeam( i ) ;
     4010        ip[0] = sref->getPol( i ) ;
     4011        sel.setIFs( ii ) ;
     4012        sel.setBeams( ib ) ;
     4013        sel.setPolarizations( ip ) ;
     4014        askyhi->setSelection( sel ) ;   
     4015        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
     4016        const Vector<Float> Vtsyshi( sptsys ) ;
     4017        tsysColhi.put( i, Vtsyshi ) ;
     4018        askyhi->unsetSelection() ;
     4019        sel.reset() ;
     4020      }
     4021    }
     4022  }
     4023  else {
     4024    // non-APEX fs data
     4025    // sky scan
     4026    STSelector sel = STSelector() ;
     4027    types.push_back( SrcType::SKY ) ;
     4028    sel.setTypes( types ) ;
     4029    s->setSelection( sel ) ;
     4030    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
     4031    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
     4032    s->unsetSelection() ;
     4033    sel.reset() ;
     4034    types.clear() ;
     4035   
     4036    // hot scan
     4037    types.push_back( SrcType::HOT ) ;
     4038    sel.setTypes( types ) ;
     4039    s->setSelection( sel ) ;
     4040    tmp.clear() ;
     4041    tmp.push_back( getScantable( s, false ) ) ;
     4042    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
     4043    s->unsetSelection() ;
     4044    sel.reset() ;
     4045    types.clear() ;
     4046
     4047    // cold scan
     4048    CountedPtr<Scantable> acold ;
     4049//     types.push_back( SrcType::COLD ) ;
     4050//     sel.setTypes( types ) ;
     4051//     s->setSelection( sel ) ;
     4052//     tmp.clear() ;
     4053//     tmp.push_back( getScantable( s, false ) ) ;
     4054//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
     4055//     s->unsetSelection() ;
     4056//     sel.reset() ;
     4057//     types.clear() ;
     4058   
     4059    // ref scan
     4060    bool insitu = insitu_ ;
     4061    insitu_ = false ;
     4062    sref = getScantable( s, true ) ;
     4063    insitu_ = insitu ;
     4064    types.push_back( SrcType::FSOFF ) ;
     4065    sel.setTypes( types ) ;
     4066    s->setSelection( sel ) ;
     4067    TableCopy::copyRows( sref->table(), s->table() ) ;
     4068    s->unsetSelection() ;
     4069    sel.reset() ;
     4070    types.clear() ;
     4071   
     4072    // sig scan
     4073    insitu_ = false ;
     4074    ssig = getScantable( s, true ) ;
     4075    insitu_ = insitu ;
     4076    types.push_back( SrcType::FSON ) ;
     4077    sel.setTypes( types ) ;
     4078    s->setSelection( sel ) ;
     4079    TableCopy::copyRows( ssig->table(), s->table() ) ;
     4080    s->unsetSelection() ;
     4081    sel.reset() ;
     4082    types.clear() ;
     4083
     4084    // process each sig and ref scan
     4085    ArrayColumn<Float> tsysColsig ;
     4086    tsysColsig.attach( ssig->table(), "TSYS" ) ;
     4087    ArrayColumn<Float> tsysColref ;
     4088    tsysColref.attach( ssig->table(), "TSYS" ) ;
     4089    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
     4090      vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
     4091      ssig->setSpectrum( sp, i ) ;
     4092      string reftime = ssig->getTime( i ) ;
     4093      vector<int> ii( 1, ssig->getIF( i ) ) ;
     4094      vector<int> ib( 1, ssig->getBeam( i ) ) ;
     4095      vector<int> ip( 1, ssig->getPol( i ) ) ;
     4096      sel.setIFs( ii ) ;
     4097      sel.setBeams( ib ) ;
     4098      sel.setPolarizations( ip ) ;
     4099      asky->setSelection( sel ) ;
     4100      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
     4101      const Vector<Float> Vtsys( sptsys ) ;
     4102      tsysColsig.put( i, Vtsys ) ;
     4103      asky->unsetSelection() ;
     4104      sel.reset() ;
     4105      sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
     4106      sref->setSpectrum( sp, i ) ;
     4107      tsysColref.put( i, Vtsys ) ;
     4108    }
     4109  }
     4110
     4111  // do folding if necessary
     4112  Table sigtab = ssig->table() ;
     4113  Table reftab = sref->table() ;
     4114  ScalarColumn<uInt> sigifnoCol ;
     4115  ScalarColumn<uInt> refifnoCol ;
     4116  ScalarColumn<uInt> sigfidCol ;
     4117  ScalarColumn<uInt> reffidCol ;
     4118  Int nchan = (Int)ssig->nchan() ;
     4119  sigifnoCol.attach( sigtab, "IFNO" ) ;
     4120  refifnoCol.attach( reftab, "IFNO" ) ;
     4121  sigfidCol.attach( sigtab, "FREQ_ID" ) ;
     4122  reffidCol.attach( reftab, "FREQ_ID" ) ;
     4123  Vector<uInt> sfids( sigfidCol.getColumn() ) ;
     4124  Vector<uInt> rfids( reffidCol.getColumn() ) ;
     4125  vector<uInt> sfids_unique ;
     4126  vector<uInt> rfids_unique ;
     4127  vector<uInt> sifno_unique ;
     4128  vector<uInt> rifno_unique ;
     4129  for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
     4130    if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
     4131      sfids_unique.push_back( sfids[i] ) ;
     4132      sifno_unique.push_back( ssig->getIF( i ) ) ;
     4133    }
     4134    if ( count( rfids_unique.begin(), rfids_unique.end(),  rfids[i] ) == 0 ) {
     4135      rfids_unique.push_back( rfids[i] ) ;
     4136      rifno_unique.push_back( sref->getIF( i ) ) ;
     4137    }
     4138  }
     4139  double refpix_sig, refval_sig, increment_sig ;
     4140  double refpix_ref, refval_ref, increment_ref ;
     4141  vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
     4142  for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
     4143    ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
     4144    sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
     4145    if ( refpix_sig == refpix_ref ) {
     4146      double foffset = refval_ref - refval_sig ;
     4147      int choffset = static_cast<int>(foffset/increment_sig) ;
     4148      double doffset = foffset / increment_sig ;
     4149      if ( abs(choffset) >= nchan ) {
     4150        LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
     4151        os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
     4152        os << "Just return signal data" << LogIO::POST ;
     4153        //std::vector< CountedPtr<Scantable> > tabs ;
     4154        //tabs.push_back( ssig ) ;
     4155        //tabs.push_back( sref ) ;
     4156        //out = merge( tabs ) ;
     4157        tmp[i] = ssig ;
     4158      }
     4159      else {
     4160        STSelector sel = STSelector() ;
     4161        vector<int> v( 1, sifno_unique[i] ) ;
     4162        sel.setIFs( v ) ;
     4163        ssig->setSelection( sel ) ;
     4164        sel.reset() ;
     4165        v[0] = rifno_unique[i] ;
     4166        sel.setIFs( v ) ;
     4167        sref->setSelection( sel ) ;
     4168        sel.reset() ;
     4169        if ( antname.find( "APEX" ) != string::npos ) {
     4170          tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
     4171          //tmp[i] = dofold( ssig, sref, doffset ) ;
     4172        }
     4173        else {
     4174          tmp[i] = dofold( ssig, sref, doffset ) ;
     4175        }
     4176        ssig->unsetSelection() ;
     4177        sref->unsetSelection() ;
     4178      }
     4179    }
     4180  }
     4181
     4182  if ( tmp.size() > 1 ) {
     4183    out = merge( tmp ) ;
     4184  }
     4185  else {
     4186    out = tmp[0] ;
     4187  }
     4188
     4189  // flux unit
     4190  out->setFluxUnit( "K" ) ;
     4191
     4192  return out ;
     4193}
     4194
     4195CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
     4196{
     4197  CountedPtr<Scantable> out ;
     4198
     4199  return out ;
     4200}
     4201
     4202vector<float> STMath::getSpectrumFromTime( string reftime,
     4203                                           CountedPtr<Scantable>& s,
     4204                                           string mode )
     4205{
     4206  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
     4207  vector<float> sp ;
     4208
     4209  if ( s->nrow() == 0 ) {
     4210    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
     4211    return sp ;
     4212  }
     4213  else if ( s->nrow() == 1 ) {
     4214    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
     4215    return s->getSpectrum( 0 ) ;
     4216  }
     4217  else {
     4218    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     4219    if ( mode == "before" ) {
     4220      int id = -1 ;
     4221      if ( idx[0] != -1 ) {
     4222        id = idx[0] ;
     4223      }
     4224      else if ( idx[1] != -1 ) {
     4225        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     4226        id = idx[1] ;
     4227      }
     4228      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4229      sp = s->getSpectrum( id ) ;
     4230    }
     4231    else if ( mode == "after" ) {
     4232      int id = -1 ;
     4233      if ( idx[1] != -1 ) {
     4234        id = idx[1] ;
     4235      }
     4236      else if ( idx[0] != -1 ) {
     4237        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     4238        id = idx[1] ;
     4239      }
     4240      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4241      sp = s->getSpectrum( id ) ;
     4242    }
     4243    else if ( mode == "nearest" ) {
     4244      int id = -1 ;
     4245      if ( idx[0] == -1 ) {
     4246        id = idx[1] ;
     4247      }
     4248      else if ( idx[1] == -1 ) {
     4249        id = idx[0] ;
     4250      }
     4251      else if ( idx[0] == idx[1] ) {
     4252        id = idx[0] ;
     4253      }
     4254      else {
     4255        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4256        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4257        double tref = getMJD( reftime ) ;
     4258        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4259          id = idx[1] ;
     4260        }
     4261        else {
     4262          id = idx[0] ;
     4263        }
     4264      }
     4265      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4266      sp = s->getSpectrum( id ) ;     
     4267    }
     4268    else if ( mode == "linear" ) {
     4269      if ( idx[0] == -1 ) {
     4270        // use after
     4271        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4272        int id = idx[1] ;
     4273        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4274        sp = s->getSpectrum( id ) ;
     4275      }
     4276      else if ( idx[1] == -1 ) {
     4277        // use before
     4278        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4279        int id = idx[0] ;
     4280        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4281        sp = s->getSpectrum( id ) ;
     4282      }
     4283      else if ( idx[0] == idx[1] ) {
     4284        // use before
     4285        //os << "No need to interporate." << LogIO::POST ;
     4286        int id = idx[0] ;
     4287        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
     4288        sp = s->getSpectrum( id ) ;
     4289      }
     4290      else {
     4291        // do interpolation
     4292        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4293        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4294        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4295        double tref = getMJD( reftime ) ;
     4296        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
     4297        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
     4298        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
     4299          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
     4300          sp.push_back( v ) ;
     4301        }
     4302      }
     4303    }
     4304    else {
     4305      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4306    }
     4307    return sp ;
     4308  }
     4309}
     4310
     4311double STMath::getMJD( string strtime )
     4312{
     4313  if ( strtime.find("/") == string::npos ) {
     4314    // MJD time string
     4315    return atof( strtime.c_str() ) ;
     4316  }
     4317  else {
     4318    // string in YYYY/MM/DD/HH:MM:SS format
     4319    uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
     4320    uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
     4321    uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
     4322    uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
     4323    uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
     4324    uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
     4325    Time t( year, month, day, hour, minute, sec ) ;
     4326    return t.modifiedJulianDay() ;
     4327  }
     4328}
     4329
     4330vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
     4331{
     4332  double reft = getMJD( reftime ) ;
     4333  double dtmin = 1.0e100 ;
     4334  double dtmax = -1.0e100 ;
     4335  vector<double> dt ;
     4336  int just_before = -1 ;
     4337  int just_after = -1 ;
     4338  for ( int i = 0 ; i < s->nrow() ; i++ ) {
     4339    dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
     4340  }
     4341  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
     4342    if ( dt[i] > 0.0 ) {
     4343      // after reftime
     4344      if ( dt[i] < dtmin ) {
     4345        just_after = i ;
     4346        dtmin = dt[i] ;
     4347      }
     4348    }
     4349    else if ( dt[i] < 0.0 ) {
     4350      // before reftime
     4351      if ( dt[i] > dtmax ) {
     4352        just_before = i ;
     4353        dtmax = dt[i] ;
     4354      }
     4355    }
     4356    else {
     4357      // just a reftime
     4358      just_before = i ;
     4359      just_after = i ;
     4360      dtmax = 0 ;
     4361      dtmin = 0 ;
     4362      break ;
     4363    }
     4364  }
     4365
     4366  vector<int> v ;
     4367  v.push_back( just_before ) ;
     4368  v.push_back( just_after ) ;
     4369
     4370  return v ;
     4371}
     4372
     4373vector<float> STMath::getTcalFromTime( string reftime,
     4374                                       CountedPtr<Scantable>& s,
     4375                                       string mode )
     4376{
     4377  LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
     4378  vector<float> tcal ;
     4379  STTcal tcalTable = s->tcal() ;
     4380  String time ;
     4381  Vector<Float> tcalval ;
     4382  if ( s->nrow() == 0 ) {
     4383    os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
     4384    return tcal ;
     4385  }
     4386  else if ( s->nrow() == 1 ) {
     4387    uInt tcalid = s->getTcalId( 0 ) ;
     4388    //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4389    tcalTable.getEntry( time, tcalval, tcalid ) ;
     4390    tcalval.tovector( tcal ) ;
     4391    return tcal ;
     4392  }
     4393  else {
     4394    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     4395    if ( mode == "before" ) {
     4396      int id = -1 ;
     4397      if ( idx[0] != -1 ) {
     4398        id = idx[0] ;
     4399      }
     4400      else if ( idx[1] != -1 ) {
     4401        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     4402        id = idx[1] ;
     4403      }
     4404      uInt tcalid = s->getTcalId( id ) ;
     4405      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4406      tcalTable.getEntry( time, tcalval, tcalid ) ;
     4407      tcalval.tovector( tcal ) ;
     4408    }
     4409    else if ( mode == "after" ) {
     4410      int id = -1 ;
     4411      if ( idx[1] != -1 ) {
     4412        id = idx[1] ;
     4413      }
     4414      else if ( idx[0] != -1 ) {
     4415        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     4416        id = idx[1] ;
     4417      }
     4418      uInt tcalid = s->getTcalId( id ) ;
     4419      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4420      tcalTable.getEntry( time, tcalval, tcalid ) ;
     4421      tcalval.tovector( tcal ) ;
     4422    }
     4423    else if ( mode == "nearest" ) {
     4424      int id = -1 ;
     4425      if ( idx[0] == -1 ) {
     4426        id = idx[1] ;
     4427      }
     4428      else if ( idx[1] == -1 ) {
     4429        id = idx[0] ;
     4430      }
     4431      else if ( idx[0] == idx[1] ) {
     4432        id = idx[0] ;
     4433      }
     4434      else {
     4435        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4436        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4437        double tref = getMJD( reftime ) ;
     4438        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4439          id = idx[1] ;
     4440        }
     4441        else {
     4442          id = idx[0] ;
     4443        }
     4444      }
     4445      uInt tcalid = s->getTcalId( id ) ;
     4446      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4447      tcalTable.getEntry( time, tcalval, tcalid ) ;
     4448      tcalval.tovector( tcal ) ;
     4449    }
     4450    else if ( mode == "linear" ) {
     4451      if ( idx[0] == -1 ) {
     4452        // use after
     4453        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4454        int id = idx[1] ;
     4455        uInt tcalid = s->getTcalId( id ) ;
     4456        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4457        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4458        tcalval.tovector( tcal ) ;
     4459      }
     4460      else if ( idx[1] == -1 ) {
     4461        // use before
     4462        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4463        int id = idx[0] ;
     4464        uInt tcalid = s->getTcalId( id ) ;
     4465        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4466        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4467        tcalval.tovector( tcal ) ;
     4468      }
     4469      else if ( idx[0] == idx[1] ) {
     4470        // use before
     4471        //os << "No need to interporate." << LogIO::POST ;
     4472        int id = idx[0] ;
     4473        uInt tcalid = s->getTcalId( id ) ;
     4474        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
     4475        tcalTable.getEntry( time, tcalval, tcalid ) ;
     4476        tcalval.tovector( tcal ) ;
     4477      }
     4478      else {
     4479        // do interpolation
     4480        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4481        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4482        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4483        double tref = getMJD( reftime ) ;
     4484        vector<float> tcal0 ;
     4485        vector<float> tcal1 ;
     4486        uInt tcalid0 = s->getTcalId( idx[0] ) ;
     4487        uInt tcalid1 = s->getTcalId( idx[1] ) ;
     4488        tcalTable.getEntry( time, tcalval, tcalid0 ) ;
     4489        tcalval.tovector( tcal0 ) ;
     4490        tcalTable.getEntry( time, tcalval, tcalid1 ) ;
     4491        tcalval.tovector( tcal1 ) ;       
     4492        for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
     4493          float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
     4494          tcal.push_back( v ) ;
     4495        }
     4496      }
     4497    }
     4498    else {
     4499      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4500    }
     4501    return tcal ;
     4502  }
     4503}
     4504
     4505vector<float> STMath::getTsysFromTime( string reftime,
     4506                                       CountedPtr<Scantable>& s,
     4507                                       string mode )
     4508{
     4509  LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
     4510  ArrayColumn<Float> tsysCol ;
     4511  tsysCol.attach( s->table(), "TSYS" ) ;
     4512  vector<float> tsys ;
     4513  String time ;
     4514  Vector<Float> tsysval ;
     4515  if ( s->nrow() == 0 ) {
     4516    os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
     4517    return tsys ;
     4518  }
     4519  else if ( s->nrow() == 1 ) {
     4520    //os << "use row " << 0 << LogIO::POST ;
     4521    tsysval = tsysCol( 0 ) ;
     4522    tsysval.tovector( tsys ) ;
     4523    return tsys ;
     4524  }
     4525  else {
     4526    vector<int> idx = getRowIdFromTime( reftime, s ) ;
     4527    if ( mode == "before" ) {
     4528      int id = -1 ;
     4529      if ( idx[0] != -1 ) {
     4530        id = idx[0] ;
     4531      }
     4532      else if ( idx[1] != -1 ) {
     4533        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
     4534        id = idx[1] ;
     4535      }
     4536      //os << "use row " << id << LogIO::POST ;
     4537      tsysval = tsysCol( id ) ;
     4538      tsysval.tovector( tsys ) ;
     4539    }
     4540    else if ( mode == "after" ) {
     4541      int id = -1 ;
     4542      if ( idx[1] != -1 ) {
     4543        id = idx[1] ;
     4544      }
     4545      else if ( idx[0] != -1 ) {
     4546        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
     4547        id = idx[1] ;
     4548      }
     4549      //os << "use row " << id << LogIO::POST ;
     4550      tsysval = tsysCol( id ) ;
     4551      tsysval.tovector( tsys ) ;
     4552    }
     4553    else if ( mode == "nearest" ) {
     4554      int id = -1 ;
     4555      if ( idx[0] == -1 ) {
     4556        id = idx[1] ;
     4557      }
     4558      else if ( idx[1] == -1 ) {
     4559        id = idx[0] ;
     4560      }
     4561      else if ( idx[0] == idx[1] ) {
     4562        id = idx[0] ;
     4563      }
     4564      else {
     4565        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4566        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4567        double tref = getMJD( reftime ) ;
     4568        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
     4569          id = idx[1] ;
     4570        }
     4571        else {
     4572          id = idx[0] ;
     4573        }
     4574      }
     4575      //os << "use row " << id << LogIO::POST ;
     4576      tsysval = tsysCol( id ) ;
     4577      tsysval.tovector( tsys ) ;
     4578    }
     4579    else if ( mode == "linear" ) {
     4580      if ( idx[0] == -1 ) {
     4581        // use after
     4582        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
     4583        int id = idx[1] ;
     4584        //os << "use row " << id << LogIO::POST ;
     4585        tsysval = tsysCol( id ) ;
     4586        tsysval.tovector( tsys ) ;
     4587      }
     4588      else if ( idx[1] == -1 ) {
     4589        // use before
     4590        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
     4591        int id = idx[0] ;
     4592        //os << "use row " << id << LogIO::POST ;
     4593        tsysval = tsysCol( id ) ;
     4594        tsysval.tovector( tsys ) ;
     4595      }
     4596      else if ( idx[0] == idx[1] ) {
     4597        // use before
     4598        //os << "No need to interporate." << LogIO::POST ;
     4599        int id = idx[0] ;
     4600        //os << "use row " << id << LogIO::POST ;
     4601        tsysval = tsysCol( id ) ;
     4602        tsysval.tovector( tsys ) ;
     4603      }
     4604      else {
     4605        // do interpolation
     4606        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
     4607        double t0 = getMJD( s->getTime( idx[0] ) ) ;
     4608        double t1 = getMJD( s->getTime( idx[1] ) ) ;
     4609        double tref = getMJD( reftime ) ;
     4610        vector<float> tsys0 ;
     4611        vector<float> tsys1 ;
     4612        tsysval = tsysCol( idx[0] ) ;
     4613        tsysval.tovector( tsys0 ) ;
     4614        tsysval = tsysCol( idx[1] ) ;
     4615        tsysval.tovector( tsys1 ) ;       
     4616        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
     4617          float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
     4618          tsys.push_back( v ) ;
     4619        }
     4620      }
     4621    }
     4622    else {
     4623      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
     4624    }
     4625    return tsys ;
     4626  }
     4627}
     4628
     4629vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
     4630                                            CountedPtr<Scantable>& off,
     4631                                            CountedPtr<Scantable>& sky,
     4632                                            CountedPtr<Scantable>& hot,
     4633                                            CountedPtr<Scantable>& cold,
     4634                                            int index,
     4635                                            string antname )
     4636{
     4637  string reftime = on->getTime( index ) ;
     4638  vector<int> ii( 1, on->getIF( index ) ) ;
     4639  vector<int> ib( 1, on->getBeam( index ) ) ;
     4640  vector<int> ip( 1, on->getPol( index ) ) ;
     4641  vector<int> ic( 1, on->getScan( index ) ) ;
     4642  STSelector sel = STSelector() ;
     4643  sel.setIFs( ii ) ;
     4644  sel.setBeams( ib ) ;
     4645  sel.setPolarizations( ip ) ;
     4646  sky->setSelection( sel ) ;
     4647  hot->setSelection( sel ) ;
     4648  //cold->setSelection( sel ) ;
     4649  off->setSelection( sel ) ;
     4650  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
     4651  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
     4652  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
     4653  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
     4654  vector<float> spec = on->getSpectrum( index ) ;
     4655  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4656  vector<float> sp( tcal.size() ) ;
     4657  if ( antname.find( "APEX" ) != string::npos ) {
     4658    // using gain array
     4659    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4660      float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
     4661        * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
     4662      sp[j] = v ;
     4663    }
     4664  }
     4665  else {
     4666    // Chopper-Wheel calibration (Ulich & Haas 1976)
     4667    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4668      float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
     4669      sp[j] = v ;
     4670    }
     4671  }
     4672  sel.reset() ;
     4673  sky->unsetSelection() ;
     4674  hot->unsetSelection() ;
     4675  //cold->unsetSelection() ;
     4676  off->unsetSelection() ;
     4677
     4678  return sp ;
     4679}
     4680
     4681vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
     4682                                            CountedPtr<Scantable>& off,
     4683                                            int index )
     4684{
     4685  string reftime = on->getTime( index ) ;
     4686  vector<int> ii( 1, on->getIF( index ) ) ;
     4687  vector<int> ib( 1, on->getBeam( index ) ) ;
     4688  vector<int> ip( 1, on->getPol( index ) ) ;
     4689  vector<int> ic( 1, on->getScan( index ) ) ;
     4690  STSelector sel = STSelector() ;
     4691  sel.setIFs( ii ) ;
     4692  sel.setBeams( ib ) ;
     4693  sel.setPolarizations( ip ) ;
     4694  off->setSelection( sel ) ;
     4695  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
     4696  vector<float> spec = on->getSpectrum( index ) ;
     4697  //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4698  //vector<float> tsys = on->getTsysVec( index ) ;
     4699  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
     4700  Vector<Float> tsys = tsysCol( index ) ;
     4701  vector<float> sp( spec.size() ) ;
     4702  // ALMA Calibration
     4703  //
     4704  // Ta* = Tsys * ( ON - OFF ) / OFF
     4705  //
     4706  // 2010/01/07 Takeshi Nakazato
     4707  unsigned int tsyssize = tsys.nelements() ;
     4708  unsigned int spsize = sp.size() ;
     4709  for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
     4710    float tscale = 0.0 ;
     4711    if ( tsyssize == spsize )
     4712      tscale = tsys[j] ;
     4713    else
     4714      tscale = tsys[0] ;
     4715    float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
     4716    sp[j] = v ;
     4717  }
     4718  sel.reset() ;
     4719  off->unsetSelection() ;
     4720
     4721  return sp ;
     4722}
     4723
     4724vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
     4725                                              CountedPtr<Scantable>& ref,
     4726                                              CountedPtr<Scantable>& sky,
     4727                                              CountedPtr<Scantable>& hot,
     4728                                              CountedPtr<Scantable>& cold,
     4729                                              int index )
     4730{
     4731  string reftime = sig->getTime( index ) ;
     4732  vector<int> ii( 1, sig->getIF( index ) ) ;
     4733  vector<int> ib( 1, sig->getBeam( index ) ) ;
     4734  vector<int> ip( 1, sig->getPol( index ) ) ;
     4735  vector<int> ic( 1, sig->getScan( index ) ) ;
     4736  STSelector sel = STSelector() ;
     4737  sel.setIFs( ii ) ;
     4738  sel.setBeams( ib ) ;
     4739  sel.setPolarizations( ip ) ;
     4740  sky->setSelection( sel ) ;
     4741  hot->setSelection( sel ) ;
     4742  //cold->setSelection( sel ) ;
     4743  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
     4744  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
     4745  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
     4746  vector<float> spref = ref->getSpectrum( index ) ;
     4747  vector<float> spsig = sig->getSpectrum( index ) ;
     4748  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
     4749  vector<float> sp( tcal.size() ) ;
     4750  for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
     4751    float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
     4752    sp[j] = v ;
     4753  }
     4754  sel.reset() ;
     4755  sky->unsetSelection() ;
     4756  hot->unsetSelection() ;
     4757  //cold->unsetSelection() ;
     4758
     4759  return sp ;
     4760}
     4761
     4762vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
     4763                                              CountedPtr<Scantable>& ref,
     4764                                              vector< CountedPtr<Scantable> >& sky,
     4765                                              vector< CountedPtr<Scantable> >& hot,
     4766                                              vector< CountedPtr<Scantable> >& cold,
     4767                                              int index )
     4768{
     4769  string reftime = sig->getTime( index ) ;
     4770  vector<int> ii( 1, sig->getIF( index ) ) ;
     4771  vector<int> ib( 1, sig->getBeam( index ) ) ;
     4772  vector<int> ip( 1, sig->getPol( index ) ) ;
     4773  vector<int> ic( 1, sig->getScan( index ) ) ;
     4774  STSelector sel = STSelector() ;
     4775  sel.setIFs( ii ) ;
     4776  sel.setBeams( ib ) ;
     4777  sel.setPolarizations( ip ) ;
     4778  sky[0]->setSelection( sel ) ;
     4779  hot[0]->setSelection( sel ) ;
     4780  //cold[0]->setSelection( sel ) ;
     4781  vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
     4782  vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
     4783  //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
     4784  vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
     4785  sel.reset() ;
     4786  ii[0] = ref->getIF( index ) ;
     4787  sel.setIFs( ii ) ;
     4788  sel.setBeams( ib ) ;
     4789  sel.setPolarizations( ip ) ;
     4790  sky[1]->setSelection( sel ) ;
     4791  hot[1]->setSelection( sel ) ;
     4792  //cold[1]->setSelection( sel ) ;
     4793  vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
     4794  vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
     4795  //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
     4796  vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 
     4797  vector<float> spref = ref->getSpectrum( index ) ;
     4798  vector<float> spsig = sig->getSpectrum( index ) ;
     4799  vector<float> sp( tcals.size() ) ;
     4800  for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
     4801    float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
     4802    sp[j] = v ;
     4803  }
     4804  sel.reset() ;
     4805  sky[0]->unsetSelection() ;
     4806  hot[0]->unsetSelection() ;
     4807  //cold[0]->unsetSelection() ;
     4808  sky[1]->unsetSelection() ;
     4809  hot[1]->unsetSelection() ;
     4810  //cold[1]->unsetSelection() ;
     4811
     4812  return sp ;
     4813}
  • branches/mergetest/src/STMath.h

    r1689 r1789  
    119119                  const std::string& mode, bool tsys=false );
    120120
     121  // array operation
     122  casa::CountedPtr<Scantable>
     123    arrayOperate( const casa::CountedPtr<Scantable>& in,
     124                  const std::vector<float> val,
     125                  const std::string& mode,
     126                  const std::string& opmode="channel", 
     127                  bool tsys=false );
     128
     129  // channel operation
     130  casa::CountedPtr<Scantable>
     131    arrayOperateChannel( const casa::CountedPtr<Scantable>& in,
     132                         const std::vector<float> val,
     133                         const std::string& mode, bool tsys=false );
     134
     135  // row operation
     136  casa::CountedPtr<Scantable>
     137    arrayOperateRow( const casa::CountedPtr<Scantable>& in,
     138                     const std::vector<float> val,
     139                     const std::string& mode, bool tsys=false );
     140
     141  // 2d array operation
     142  casa::CountedPtr<Scantable>
     143    array2dOperate( const casa::CountedPtr<Scantable>& in,
     144                  const std::vector< std::vector<float> > val,
     145                  const std::string& mode, bool tsys=false );
     146
    121147  casa::CountedPtr<Scantable>
    122148    binaryOperate( const casa::CountedPtr<Scantable>& left,
     
    134160  /**
    135161    * Calibrate total power scans (translated from GBTIDL)
    136     * @param calon uncalibrated Scantable with CAL noise signal
     162    * @param calon uncalibrated Scantable with CAL noise signal 
    137163    * @param caloff uncalibrated Scantable with no CAL signal
    138164    * @param tcal optional scalar Tcal, CAL temperature (K)
    139     * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
     165    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable 
    140166    * (spectrum - average of the two CAL on and off spectra;
    141167    * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
    142     */
     168    */             
    143169  casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
    144170                                            const casa::CountedPtr<Scantable>& caloff,
     
    151177    * @param smoothref optional Boxcar smooth width of the reference scans
    152178    * default: no smoothing (=1)
    153     * @param tsysv optional scalar Tsys value at the zenith, required to
    154     * set tau, as well
     179    * @param tsysv optional scalar Tsys value at the zenith, required to 
     180    * set tau, as well 
    155181    * @param tau optional scalar Tau value
    156182    * @return casa::CountedPtr<Scantable> which holds combined scans
     
    163189                                        casa::Float tau=0.0 );
    164190
    165  /**
     191  /**
    166192    * Calibrate GBT Nod scan pairs (translated from GBTIDL)
    167193    * @param s Scantable which contains Nod scans
     
    170196    * @param tsysv optional scalar Tsys value at the zenith, required to
    171197    * set tau, as well
    172     * @param tau optional scalar Tau value
     198    * @param tau optional scalar Tau value 
    173199    * @param tcal optional scalar Tcal, CAL temperature (K)
    174200    * @return casa::CountedPtr<Scantable> which holds calibrated scans
     
    199225                                    casa::Float tcal=0.0 );
    200226
     227  /**
     228   * Calibrate data with Chopper-Wheel like calibration method
     229   * which adopts position switching by antenna motion,
     230   * wobbler (nutator) switching and On-The-Fly observation.
     231   *
     232   * The method is applicable to APEX, and other telescopes other than GBT.
     233   *
     234   * @param a Scantable which contains ON and OFF scans
     235   * @param a string that indicates calibration mode
     236   * @param a string that indicates antenna name
     237   **/
     238  casa::CountedPtr<Scantable> cwcal( const casa::CountedPtr<Scantable>& s,
     239                                       const casa::String calmode,
     240                                       const casa::String antname );
     241
     242  /**
     243   * Calibrate frequency switched scans with Chopper-Wheel like
     244   * calibration method.
     245   *
     246   * The method is applicable to APEX, and other telescopes other than GBT.
     247   *
     248   * @param a Scantable which contains ON and OFF scans
     249   * @param a string that indicates antenna name
     250   **/
     251  casa::CountedPtr<Scantable> cwcalfs( const casa::CountedPtr<Scantable>& s,
     252                                       const casa::String antname );
     253
     254
     255  /**
     256   * Folding frequency-switch data
     257   * @param sig
     258   * @param ref
     259   * @param choffset
     260   **/
     261  casa::CountedPtr<Scantable> dofold( const casa::CountedPtr<Scantable> &sig,
     262                                      const casa::CountedPtr<Scantable> &ref,
     263                                      casa::Double choffset,
     264                                      casa::Double choffset = 0.0 );
     265
     266  /**
     267   * ALMA calibration
     268   **/
     269  casa::CountedPtr<Scantable> almacal( const casa::CountedPtr<Scantable>& s,
     270                                       const casa::String calmode ) ;
     271  casa::CountedPtr<Scantable> almacalfs( const casa::CountedPtr<Scantable>& s ) ;
    201272
    202273  casa::CountedPtr<Scantable>
     
    206277                               const std::vector<bool>& mask,
    207278                               const std::string& which);
     279
     280  std::vector< int > minMaxChan(const casa::CountedPtr<Scantable>& in,
     281                                const std::vector<bool>& mask,
     282                                const std::string& which);
    208283
    209284  casa::CountedPtr<Scantable> bin( const casa::CountedPtr<Scantable>& in,
     
    266341             double end, const std::string& mode="frequency");
    267342
     343  // test for average spectra with different channel/resolution
     344  casa::CountedPtr<Scantable>
     345    new_average( const std::vector<casa::CountedPtr<Scantable> >& in,
     346                 const bool& compel,
     347                 const std::vector<bool>& mask = std::vector<bool>(),
     348                 const std::string& weight = "NONE",
     349                 const std::string& avmode = "SCAN" )
     350    throw (casa::AipsError) ;
     351
    268352private:
    269353  casa::CountedPtr<Scantable>  applyToPol( const casa::CountedPtr<Scantable>& in,
     
    301385    maskedArray( const casa::Vector<casa::Float>& s,
    302386                 const casa::Vector<casa::uChar>& f );
     387  casa::MaskedArray<casa::Double>
     388    maskedArray( const casa::Vector<casa::Double>& s,
     389                 const casa::Vector<casa::uChar>& f );
    303390  casa::Vector<casa::uChar>
    304391    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
    305392
     393  vector<float> getSpectrumFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode = "before" ) ;
     394  vector<float> getTcalFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     395  vector<float> getTsysFromTime( string reftime, casa::CountedPtr<Scantable>& s, string mode="before" ) ;
     396  vector<int> getRowIdFromTime( string reftime, casa::CountedPtr<Scantable>& s ) ;
     397
     398  // Chopper-Wheel type calibration
     399  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
     400                                      casa::CountedPtr<Scantable>& off,
     401                                      casa::CountedPtr<Scantable>& sky,
     402                                      casa::CountedPtr<Scantable>& hot,
     403                                      casa::CountedPtr<Scantable>& cold,
     404                                      int index,
     405                                      string antname ) ;
     406  // Tsys * (ON-OFF)/OFF
     407  vector<float> getCalibratedSpectra( casa::CountedPtr<Scantable>& on,
     408                                      casa::CountedPtr<Scantable>& off,
     409                                      int index ) ;
     410  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
     411                                        casa::CountedPtr<Scantable>& ref,
     412                                        casa::CountedPtr<Scantable>& sky,
     413                                        casa::CountedPtr<Scantable>& hot,
     414                                        casa::CountedPtr<Scantable>& cold,
     415                                        int index ) ;
     416  vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
     417                                        casa::CountedPtr<Scantable>& ref,
     418                                        vector< casa::CountedPtr<Scantable> >& sky,
     419                                        vector< casa::CountedPtr<Scantable> >& hot,
     420                                        vector< casa::CountedPtr<Scantable> >& cold,
     421                                        int index ) ;
     422  double getMJD( string strtime ) ;
     423
    306424  bool insitu_;
    307425};
  • branches/mergetest/src/STMathWrapper.h

    r1689 r1789  
    7373  { return ScantableWrapper(STMath::unaryOperate(in.getCP(), val, mode, tsys)); }
    7474
     75  ScantableWrapper arrayOperate( const ScantableWrapper& in,
     76                                 const std::vector<float> val,
     77                                 const std::string& mode,
     78                                 bool tsys=false )
     79  { return ScantableWrapper(STMath::arrayOperateChannel(in.getCP(), val, mode, tsys)); }
     80
     81  ScantableWrapper array2dOperate( const ScantableWrapper& in,
     82                                   const std::vector< std::vector<float> > val,
     83                                   const std::string& mode, bool tsys=false )
     84  { return ScantableWrapper(STMath::array2dOperate(in.getCP(), val, mode, tsys)); }
     85
    7586  ScantableWrapper binaryOperate( const ScantableWrapper& left,
    7687                                  const ScantableWrapper& right,
     
    121132  { return STMath::statistic(in.getCP(), mask, which); }
    122133
     134  std::vector<int> minMaxChan(const ScantableWrapper& in,
     135                               const std::vector<bool>& mask,
     136                               const std::string& which)
     137  { return STMath::minMaxChan(in.getCP(), mask, which); }
     138
    123139  ScantableWrapper bin( const ScantableWrapper& in, int width=5)
    124140  { return ScantableWrapper(STMath::bin(in.getCP(), width)); }
     
    175191                                   const std::string& refTime,
    176192                                   const std::string& method  )
    177   { return ScantableWrapper(STMath::frequencyAlign(in.getCP(), refTime, method)); }
     193  { return ScantableWrapper(STMath::frequencyAlign(in.getCP())); }
    178194
    179195  ScantableWrapper convertPolarisation( const ScantableWrapper& in,
     
    191207                                            mode)); }
    192208
     209  // test for average spectra with different channel/resolution
     210  ScantableWrapper
     211    new_average( const std::vector<ScantableWrapper>& in,
     212                 const bool& compel,
     213                 const std::vector<bool>& mask,
     214                 const std::string& weight,
     215                 const std::string& avmode )
     216  {
     217    std::vector<casa::CountedPtr<Scantable> > sts;
     218    for (unsigned int i=0; i<in.size(); ++i) sts.push_back(in[i].getCP());
     219    return ScantableWrapper(STMath::new_average(sts, compel, mask, weight, avmode));
     220  }
     221
     222  // cwcal
     223  ScantableWrapper cwcal( const ScantableWrapper &in,
     224                          const std::string calmode,
     225                          const std::string antname )
     226  {
     227    casa::CountedPtr<Scantable> tab = in.getCP() ;
     228    casa::String mode( calmode ) ;
     229    casa::String name( antname ) ;
     230    return ScantableWrapper( STMath::cwcal( tab, mode, name ) ) ;
     231  }
     232  // almacal
     233  ScantableWrapper almacal( const ScantableWrapper &in,
     234                          const std::string calmode )
     235  {
     236    casa::CountedPtr<Scantable> tab = in.getCP() ;
     237    casa::String mode( calmode ) ;
     238    return ScantableWrapper( STMath::almacal( tab, mode ) ) ;
     239  }
    193240};
    194241
  • branches/mergetest/src/STMolecules.cpp

    r870 r1789  
    1414#include <tables/Tables/SetupNewTab.h>
    1515#include <tables/Tables/ScaColDesc.h>
     16#include <tables/Tables/ArrColDesc.h>
    1617#include <tables/Tables/TableRecord.h>
    1718#include <tables/Tables/TableParse.h>
    1819#include <tables/Tables/TableRow.h>
    1920#include <casa/Containers/RecordField.h>
     21
     22#include <tables/Tables/TableProxy.h>
    2023
    2124#include "STMolecules.h"
     
    5962{
    6063  // add to base class table
    61   table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
    62   table_.addColumn(ScalarColumnDesc<String>("NAME"));
    63   table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
     64  //table_.addColumn(ScalarColumnDesc<Double>("RESTFREQUENCY"));
     65  table_.addColumn(ArrayColumnDesc<Double>("RESTFREQUENCY"));
     66  //table_.addColumn(ScalarColumnDesc<String>("NAME"));
     67  table_.addColumn(ArrayColumnDesc<String>("NAME"));
     68  //table_.addColumn(ScalarColumnDesc<String>("FORMATTEDNAME"));
     69  table_.addColumn(ArrayColumnDesc<String>("FORMATTEDNAME"));
    6470  table_.rwKeywordSet().define("UNIT", String("Hz"));
    6571  // new cached columns
     
    6975}
    7076
     77/***
    7178uInt STMolecules::addEntry( Double restfreq, const String& name,
    7279                            const String& formattedname )
     
    94101  return resultid;
    95102}
    96 
     103***/
     104uInt STMolecules::addEntry( Vector<Double> restfreq, const Vector<String>& name,
     105                            const Vector<String>& formattedname )
     106{
     107// How to handle this...?
     108  Table result =
     109    table_( nelements(table_.col("RESTFREQUENCY")) == uInt (restfreq.size()) &&
     110            all(table_.col("RESTFREQUENCY")== restfreq) );
     111  uInt resultid = 0;
     112  if ( result.nrow() > 0) {
     113    ROScalarColumn<uInt> c(result, "ID");
     114    c.get(0, resultid);
     115  } else {
     116    uInt rno = table_.nrow();
     117    table_.addRow();
     118    // get last assigned _id and increment
     119    if ( rno > 0 ) {
     120      idCol_.get(rno-1, resultid);
     121      resultid++;
     122    }
     123    restfreqCol_.put(rno, restfreq);
     124    nameCol_.put(rno, name);
     125    formattednameCol_.put(rno, formattedname);
     126    idCol_.put(rno, resultid);
     127  }
     128  return resultid;
     129}
     130
     131
     132
     133
     134/***
    97135void STMolecules::getEntry( Double& restfreq, String& name,
    98136                            String& formattedname, uInt id ) const
     
    104142  ROTableRow row(t);
    105143  // get first row - there should only be one matching id
     144
    106145  const TableRecord& rec = row.get(0);
    107146  restfreq = rec.asDouble("RESTFREQUENCY");
     
    109148  formattedname = rec.asString("FORMATTEDNAME");
    110149}
     150***/
     151void STMolecules::getEntry( Vector<Double>& restfreq, Vector<String>& name,
     152                            Vector<String>& formattedname, uInt id ) const
     153{
     154  Table t = table_(table_.col("ID") == Int(id) );
     155  if (t.nrow() == 0 ) {
     156    throw(AipsError("STMolecules::getEntry - id out of range"));
     157  }
     158  ROTableRow row(t);
     159  // get first row - there should only be one matching id
     160
     161  const TableRecord& rec = row.get(0);
     162  //restfreq = rec.asDouble("RESTFREQUENCY");
     163  restfreq = rec.asArrayDouble("RESTFREQUENCY");
     164  //name = rec.asString("NAME");
     165  name = rec.asArrayString("NAME");
     166  //formattedname = rec.asString("FORMATTEDNAME");
     167  formattedname = rec.asArrayString("FORMATTEDNAME");
     168}
    111169
    112170std::vector< double > asap::STMolecules::getRestFrequencies( ) const
    113171{
    114172  std::vector<double> out;
     173  //TableProxy itsTable(table_);
     174  //Record rec;
    115175  Vector<Double> rfs = restfreqCol_.getColumn();
    116176  rfs.tovector(out);
     177  //rec = itsTable.getVarColumn("RESTFREQUENCY", 0, 1, 1);
    117178  return out;
    118179}
    119180
    120 double asap::STMolecules::getRestFrequency( uInt id ) const
    121 {
     181std::vector< double > asap::STMolecules::getRestFrequency( uInt id ) const
     182{
     183  std::vector<double> out;
    122184  Table t = table_(table_.col("ID") == Int(id) );
    123185  if (t.nrow() == 0 ) {
     
    126188  ROTableRow row(t);
    127189  const TableRecord& rec = row.get(0);
    128   return double(rec.asDouble("RESTFREQUENCY"));
    129 }
    130 
     190  //return double(rec.asDouble("RESTFREQUENCY"));
     191  Vector<Double> rfs = rec.asArrayDouble("RESTFREQUENCY");
     192  rfs.tovector(out);
     193  return out;
     194}
     195
     196int asap::STMolecules::nrow() const
     197{
     198  return int(table_.nrow());
     199}
    131200
    132201}//namespace
  • branches/mergetest/src/STMolecules.h

    r1353 r1789  
    1717#include <tables/Tables/Table.h>
    1818#include <tables/Tables/ScalarColumn.h>
     19#include <tables/Tables/ArrayColumn.h>
     20#include <casa/Arrays/Array.h>
    1921
    2022#include "STSubTable.h"
     
    3739  STMolecules& operator=(const STMolecules& other);
    3840
     41/***
    3942  casa::uInt addEntry( casa::Double restfreq, const casa::String& name="",
    4043                       const casa::String& formattedname="");
     44***/
    4145
     46  casa::uInt addEntry( casa::Vector<casa::Double> restfreq, const casa::Vector<casa::String>& name=casa::Vector<casa::String>(0),
     47                       const casa::Vector<casa::String>& formattedname=casa::Vector<casa::String>(0));
     48
     49/***
    4250  void getEntry( casa::Double& restfreq, casa::String& name,
    4351                 casa::String& formattedname, casa::uInt id) const;
     52***/
     53  void getEntry( casa::Vector<casa::Double>& restfreq, casa::Vector<casa::String>& name,
     54                 casa::Vector<casa::String>& formattedname, casa::uInt id) const;
    4455
    4556  std::vector<double> getRestFrequencies() const;
    46   double getRestFrequency( casa::uInt id ) const;
     57  std::vector<double> getRestFrequency( casa::uInt id ) const;
    4758  const casa::String& name() const { return name_; }
     59  int nrow() const;
    4860
    4961private:
     
    5264  //casa::Table table_;
    5365  //casa::ScalarColumn<casa::uInt> freqidCol_;
    54   casa::ScalarColumn<casa::Double> restfreqCol_;
    55   casa::ScalarColumn<casa::String> nameCol_;
    56   casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
     66  //casa::ScalarColumn<casa::Double> restfreqCol_;
     67  casa::ArrayColumn<casa::Double> restfreqCol_;
     68  //casa::ScalarColumn<casa::String> nameCol_;
     69  casa::ArrayColumn<casa::String> nameCol_;
     70  //casa::ScalarColumn<casa::String> formattednameCol_; // e.g. latex
     71  casa::ArrayColumn<casa::String> formattednameCol_; // e.g. latex
    5772
    5873};
  • branches/mergetest/src/STSelector.cpp

    r1542 r1789  
    8787}
    8888
     89void STSelector::setTypes( const std::vector< int >& types )
     90{
     91  setint("SRCTYPE", types);
     92}
     93
    8994void STSelector::setint(const std::string& key, const std::vector< int >& val)
    9095{
     
    116121}
    117122
     123void STSelector::setRows( const std::vector< int >& rows )
     124{
     125  rowselection_ = rows;
     126}
     127
     128// Table STSelector::apply( const Table& tab )
     129// {
     130//   if ( empty() ) {
     131//     return sort(tab);
     132//   }
     133//   TableExprNode query;
     134//   intidmap::const_iterator it;
     135//   for (it = intselections_.begin(); it != intselections_.end(); ++it) {
     136//     TableExprNode theset(Vector<Int>( (*it).second ));
     137//     if ( query.isNull() ) {
     138//       query = tab.col((*it).first).in(theset);
     139//     } else {
     140//       query = tab.col((*it).first).in(theset) && query;
     141//     }
     142//   }
     143//   stringidmap::const_iterator it1;
     144//   for (it1 = stringselections_.begin(); it1 != stringselections_.end(); ++it1) {
     145//     TableExprNode theset(mathutil::toVectorString( (*it1).second ));
     146//     if ( query.isNull() ) {
     147//       query = tab.col((*it1).first).in(theset);
     148//     } else {
     149//       query = tab.col((*it1).first).in(theset) && query;
     150//     }
     151//   }
     152//   // add taql query
     153//   if ( taql_.size() > 0 ) {
     154//     Table tmpt = tab;
     155//     std::string pytaql = "USING STYLE PYTHON " + taql_;
     156
     157//     if ( !query.isNull() ) { // taql and selection
     158//       tmpt = tableCommand(pytaql, tab(query));
     159//     } else { // taql only
     160//       tmpt = tableCommand(pytaql, tab);
     161//     }
     162//     return sort(tmpt);
     163//   } else {
     164//     if ( query.isNull() ) {
     165//       return sort(tab);
     166//     } else {
     167//       return sort(tab(query));
     168//     }
     169//   }
     170// }
    118171Table STSelector::apply( const Table& tab )
    119172{
    120173  if ( empty() ) {
    121174    return sort(tab);
     175  }
     176  Table basetab = tab;
     177  // Important!! Be sure to apply row selection first.
     178  if (rowselection_.size() > 0){
     179    //Vector<Int> intrownrs(rowselection_);
     180    Vector<uInt> rownrs( rowselection_.size() );
     181    convertArray(rownrs, Vector<Int> ( rowselection_ ));
     182    basetab = tab( rownrs );
     183    ///TableExprNode theset(Vector<Int>( rowselection_ ));
     184    ///query = tab.nodeRownr().in(theset);
    122185  }
    123186  TableExprNode query;
     
    126189    TableExprNode theset(Vector<Int>( (*it).second ));
    127190    if ( query.isNull() ) {
    128       query = tab.col((*it).first).in(theset);
     191      //query = tab.col((*it).first).in(theset);
     192      query = basetab.col((*it).first).in(theset);
    129193    } else {
    130       query = tab.col((*it).first).in(theset) && query;
     194      //query = tab.col((*it).first).in(theset) && query;
     195      query = basetab.col((*it).first).in(theset) && query;
    131196    }
    132197  }
     
    135200    TableExprNode theset(mathutil::toVectorString( (*it1).second ));
    136201    if ( query.isNull() ) {
    137       query = tab.col((*it1).first).in(theset);
     202      //query = tab.col((*it1).first).in(theset);
     203      query = basetab.col((*it1).first).in(theset);
    138204    } else {
    139       query = tab.col((*it1).first).in(theset) && query;
     205      //query = tab.col((*it1).first).in(theset) && query;
     206      query = basetab.col((*it1).first).in(theset) && query;
    140207    }
    141208  }
    142209  // add taql query
    143210  if ( taql_.size() > 0 ) {
    144     Table tmpt = tab;
     211    //Table tmpt = tab;
     212    Table tmpt = basetab;
    145213    std::string pytaql = "USING STYLE PYTHON " + taql_;
    146214
    147215    if ( !query.isNull() ) { // taql and selection
    148       tmpt = tableCommand(pytaql, tab(query));
     216      //tmpt = tableCommand(pytaql, tab(query));
     217      tmpt = tableCommand(pytaql, basetab(query));
    149218    } else { // taql only
    150       tmpt = tableCommand(pytaql, tab);
     219      //tmpt = tableCommand(pytaql, tab);
     220      tmpt = tableCommand(pytaql, basetab);
    151221    }
    152222    return sort(tmpt);
    153223  } else {
    154224    if ( query.isNull() ) {
    155       return sort(tab);
     225      //return sort(tab);
     226      return sort(basetab);
    156227    } else {
    157       return sort(tab(query));
     228      //return sort(tab(query));
     229      return sort(basetab(query));
    158230    }
    159231  }
     
    191263{
    192264  return getint("CYCLENO");
     265}
     266
     267std::vector< int > asap::STSelector::getTypes( ) const
     268{
     269  return getint("SRCTYPE") ;
    193270}
    194271
     
    227304bool asap::STSelector::empty( ) const
    228305{
    229   return (intselections_.empty() && taql_.size() == 0 );
     306  //return (intselections_.empty() && taql_.size() == 0 );
     307  return (intselections_.empty() && taql_.size() == 0 && rowselection_.size() == 0);
    230308}
    231309
  • branches/mergetest/src/STSelector.h

    r939 r1789  
    4545  void setCycles(const std::vector<int>& cycs);
    4646  void setName(const std::string&);
     47  void setTypes(const std::vector<int>& types);
    4748  virtual void setTaQL(const std::string& taql);
    4849
    4950  void setSortOrder(const std::vector<std::string>& order);
     51  void setRows(const std::vector<int>& rows);
    5052
    5153  std::vector<int> getScans() const;
     
    5456  std::vector<int> getPols() const;
    5557  std::vector<int> getCycles() const;
     58  std::vector<int> getTypes() const;
    5659  std::vector<std::string> getPolTypes() const;
    5760  std::string getTaQL() const { return taql_; }
     
    8689  casa::Block<casa::String> order_;
    8790  std::string taql_;
     91  std::vector<int> rowselection_;
    8892};
    8993
  • branches/mergetest/src/STWriter.cpp

    r1688 r1789  
    4242#include <atnf/PKSIO/PKSMS2writer.h>
    4343#include <atnf/PKSIO/PKSSDwriter.h>
     44#include <atnf/PKSIO/SrcType.h>
    4445
    4546#include <tables/Tables/Table.h>
     
    154155    havexpol(ifs[i]) = nPol(ifs[i]) > 2;
    155156  }
    156   Vector<String> obstypes(2);
    157   obstypes(0) = "TR";//on
    158   obstypes(1) = "RF TR";//off
     157//   Vector<String> obstypes(2);
     158//   obstypes(0) = "TR";//on
     159//   obstypes(1) = "RF TR";//off
    159160  const Table table = inst->table();
    160161
     
    182183    while (!beamit.pastEnd() ) {
    183184      Table btable = beamit.table();
     185      //MDirection::ScalarColumn dirCol(btable, "DIRECTION");
     186      //pksrec.direction = dirCol(0).getAngle("rad").getValue();
    184187      TableIterator cycit(btable, "CYCLENO");
    185188      ROArrayColumn<Double> srateCol(btable, "SCANRATE");
     
    188191      Vector<Float> srateflt(sratedbl.nelements());
    189192      convertArray(srateflt, sratedbl);
    190       pksrec.scanRate = srateflt;
     193      //pksrec.scanRate = srateflt;
     194      pksrec.scanRate = sratedbl;
    191195      ROArrayColumn<Double> spmCol(btable, "SRCPROPERMOTION");
    192196      spmCol.get(0, pksrec.srcPM);
     
    200204      while (!cycit.pastEnd() ) {
    201205        Table ctable = cycit.table();
    202         TableIterator ifit(ctable, "IFNO");
     206        TableIterator ifit(ctable, "IFNO", TableIterator::Ascending, TableIterator::HeapSort);
    203207        MDirection::ScalarColumn dirCol(ctable, "DIRECTION");
    204208        pksrec.direction = dirCol(0).getAngle("rad").getValue();
     
    213217          uInt nchan = specCol(0).nelements();
    214218          Double crval,crpix;
     219          //Vector<Double> restfreq;
    215220          Float tmp0,tmp1,tmp2,tmp3,tmp4;
    216           String stmp0,stmp1;
     221          String tcalt;
     222          Vector<String> stmp0, stmp1;
    217223          inst->frequencies().getEntry(crpix,crval, pksrec.freqInc,
    218224                                     rec.asuInt("FREQ_ID"));
     
    239245          pksrec.fieldName = rec.asString("FIELDNAME");
    240246          pksrec.srcName   = rec.asString("SRCNAME");
    241           pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
     247          //pksrec.obsType   = obstypes[rec.asInt("SRCTYPE")];
     248          pksrec.obsType = getObsTypes( rec.asInt("SRCTYPE") ) ;
    242249          pksrec.bandwidth = nchan * abs(pksrec.freqInc);
    243250          pksrec.azimuth   = rec.asFloat("AZIMUTH");
     
    253260          pksrec.baseSub   = 0.0f;
    254261          pksrec.xCalFctr  = 0.0;
     262          pksrec.flagrow = rec.asuInt("FLAGROW");
    255263
    256264          status = writer_->write(pksrec);
     
    348356}
    349357
    350 }
     358// get obsType string from SRCTYPE value
     359String STWriter::getObsTypes( Int srctype )
     360{
     361  String obsType ;
     362  switch( srctype ) {
     363  case Int(SrcType::PSON):
     364    obsType = "PSON" ;
     365    break ;
     366  case Int(SrcType::PSOFF):
     367    obsType = "PSOFF" ;
     368    break ;
     369  case Int(SrcType::NOD):
     370    obsType = "NOD" ;
     371    break ;
     372  case Int(SrcType::FSON):
     373    obsType = "FSON" ;
     374    break ;
     375  case Int(SrcType::FSOFF):
     376    obsType = "FSOFF" ;
     377    break ;
     378  case Int(SrcType::SKY):
     379    obsType = "SKY" ;
     380    break ;
     381  case Int(SrcType::HOT):
     382    obsType = "HOT" ;
     383    break ;
     384  case Int(SrcType::WARM):
     385    obsType = "WARM" ;
     386    break ;
     387  case Int(SrcType::COLD):
     388    obsType = "COLD" ;
     389    break ;
     390  case Int(SrcType::PONCAL):
     391    obsType = "PSON:CALON" ;
     392    break ;
     393  case Int(SrcType::POFFCAL):
     394    obsType = "PSOFF:CALON" ;
     395    break ;
     396  case Int(SrcType::NODCAL):
     397    obsType = "NOD:CALON" ;
     398    break ;
     399  case Int(SrcType::FONCAL):
     400    obsType = "FSON:CALON" ;
     401    break ;
     402  case Int(SrcType::FOFFCAL):
     403    obsType = "FSOFF:CALOFF" ;
     404    break ;
     405  case Int(SrcType::FSLO):
     406    obsType = "FSLO" ;
     407    break ;
     408  case Int(SrcType::FLOOFF):
     409    obsType = "FS:LOWER:OFF" ;
     410    break ;
     411  case Int(SrcType::FLOSKY):
     412    obsType = "FS:LOWER:SKY" ;
     413    break ;
     414  case Int(SrcType::FLOHOT):
     415    obsType = "FS:LOWER:HOT" ;
     416    break ;
     417  case Int(SrcType::FLOWARM):
     418    obsType = "FS:LOWER:WARM" ;
     419    break ;
     420  case Int(SrcType::FLOCOLD):
     421    obsType = "FS:LOWER:COLD" ;
     422    break ;
     423  case Int(SrcType::FSHI):
     424    obsType = "FSHI" ;
     425    break ;
     426  case Int(SrcType::FHIOFF):
     427    obsType = "FS:HIGHER:OFF" ;
     428    break ;
     429  case Int(SrcType::FHISKY):
     430    obsType = "FS:HIGHER:SKY" ;
     431    break ;
     432  case Int(SrcType::FHIHOT):
     433    obsType = "FS:HIGHER:HOT" ;
     434    break ;
     435  case Int(SrcType::FHIWARM):
     436    obsType = "FS:HIGHER:WARM" ;
     437    break ;
     438  case Int(SrcType::FHICOLD):
     439    obsType = "FS:HIGHER:COLD" ;
     440    break ;
     441  default:
     442    obsType = "NOTYPE" ;
     443  }
     444
     445  return obsType ;
     446}
     447
     448}
  • branches/mergetest/src/STWriter.h

    r1391 r1789  
    3636#include <casa/aips.h>
    3737#include <casa/Utilities/CountedPtr.h>
     38#include <casa/BasicSL/String.h>
    3839
    3940#include "Logger.h"
     
    8485  void replacePtTab(const casa::Table& tab, const std::string& fname);
    8586
     87  casa::String getObsTypes( casa::Int srctype ) ;
     88
    8689  std::string     format_;
    8790  PKSwriter* writer_;
  • branches/mergetest/src/Scantable.cpp

    r1743 r1789  
    1111//
    1212#include <map>
     13#include <fstream>
    1314
    1415#include <casa/aips.h>
     
    2425#include <casa/Arrays/Vector.h>
    2526#include <casa/Arrays/VectorSTLIterator.h>
     27#include <casa/Arrays/Slice.h>
    2628#include <casa/BasicMath/Math.h>
    2729#include <casa/BasicSL/Constants.h>
     
    2931#include <casa/Containers/RecordField.h>
    3032#include <casa/Utilities/GenSort.h>
     33#include <casa/Logging/LogIO.h>
    3134
    3235#include <tables/Tables/TableParse.h>
     
    106109{
    107110  initFactories();
     111
    108112  Table tab(name, Table::Update);
    109113  uInt version = tab.keywordSet().asuInt("VERSION");
     
    121125  attach();
    122126}
     127/*
     128Scantable::Scantable(const std::string& name, Table::TableType ttype) :
     129  type_(ttype)
     130{
     131  initFactories();
     132  Table tab(name, Table::Update);
     133  uInt version = tab.keywordSet().asuInt("VERSION");
     134  if (version != version_) {
     135    throw(AipsError("Unsupported version of ASAP file."));
     136  }
     137  if ( type_ == Table::Memory ) {
     138    table_ = tab.copyToMemoryTable(generateName());
     139  } else {
     140    table_ = tab;
     141  }
     142
     143  attachSubtables();
     144  originalTable_ = table_;
     145  attach();
     146}
     147*/
    123148
    124149Scantable::Scantable( const Scantable& other, bool clear )
     
    201226  td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
    202227
     228  td.addColumn(ScalarColumnDesc<uInt>("FLAGROW"));
     229
    203230  td.addColumn(ScalarColumnDesc<Double>("TIME"));
    204231  TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
     
    257284  originalTable_ = table_;
    258285}
    259 
    260286
    261287void Scantable::attach()
     
    285311  mfocusidCol_.attach(table_, "FOCUS_ID");
    286312  mmolidCol_.attach(table_, "MOLECULE_ID");
     313
     314  //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
     315  attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
     316
     317}
     318
     319template<class T, class T2>
     320void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
     321                                   const String& colName,
     322                                   const T2& defValue)
     323{
     324  try {
     325    col.attach(table_, colName);
     326  } catch (TableError& err) {
     327    String errMesg = err.getMesg();
     328    if (errMesg == "Table column " + colName + " is unknown") {
     329      table_.addColumn(ScalarColumnDesc<T>(colName));
     330      col.attach(table_, colName);
     331      col.fillColumn(static_cast<T>(defValue));
     332    } else {
     333      throw;
     334    }
     335  } catch (...) {
     336    throw;
     337  }
     338}
     339
     340template<class T, class T2>
     341void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
     342                                   const String& colName,
     343                                   const Array<T2>& defValue)
     344{
     345  try {
     346    col.attach(table_, colName);
     347  } catch (TableError& err) {
     348    String errMesg = err.getMesg();
     349    if (errMesg == "Table column " + colName + " is unknown") {
     350      table_.addColumn(ArrayColumnDesc<T>(colName));
     351      col.attach(table_, colName);
     352
     353      int size = 0;
     354      ArrayIterator<T2>& it = defValue.begin();
     355      while (it != defValue.end()) {
     356        ++size;
     357        ++it;
     358      }
     359      IPosition ip(1, size);
     360      Array<T>& arr(ip);
     361      for (int i = 0; i < size; ++i)
     362        arr[i] = static_cast<T>(defValue[i]);
     363     
     364      col.fillColumn(arr);
     365    } else {
     366      throw;
     367    }
     368  } catch (...) {
     369    throw;
     370  }
    287371}
    288372
     
    627711}
    628712
     713void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
     714{
     715  for (uInt i=0; i<table_.nrow(); ++i) {
     716    Vector<uChar> flgs = flagsCol_(i);
     717    srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
     718    flagsCol_.put(i, flgs);
     719  }
     720}
     721
     722std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
     723{
     724  Vector<uChar> flags;
     725  flagsCol_.get(uInt(whichrow), flags);
     726  srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
     727  Vector<Bool> bflag(flags.shape());
     728  convertArray(bflag, flags);
     729  //bflag = !bflag;
     730
     731  std::vector<bool> mask;
     732  bflag.tovector(mask);
     733  return mask;
     734}
     735
     736void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
     737                                   Vector<uChar> flgs)
     738{
     739    Vector<Float> spcs = specCol_(whichrow);
     740    uInt nchannel = nchan();
     741    if (spcs.nelements() != nchannel) {
     742      throw(AipsError("Data has incorrect number of channels"));
     743    }
     744    uChar userflag = 1 << 7;
     745    if (unflag) {
     746      userflag = 0 << 7;
     747    }
     748    if (clipoutside) {
     749      for (uInt j = 0; j < nchannel; ++j) {
     750        Float spc = spcs(j);
     751        if ((spc >= uthres) || (spc <= dthres)) {
     752          flgs(j) = userflag;
     753        }
     754      }
     755    } else {
     756      for (uInt j = 0; j < nchannel; ++j) {
     757        Float spc = spcs(j);
     758        if ((spc < uthres) && (spc > dthres)) {
     759          flgs(j) = userflag;
     760        }
     761      }
     762    }
     763}
     764
    629765void Scantable::flag(const std::vector<bool>& msk, bool unflag)
    630766{
     
    673809    flagsCol_.put(i, flgs);
    674810  }
     811}
     812
     813void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
     814{
     815  if ( selector_.empty() && (rows.size() == table_.nrow()) )
     816    throw(AipsError("Trying to flag whole scantable."));
     817
     818  uInt rowflag = (unflag ? 0 : 1);
     819  std::vector<uInt>::const_iterator it;
     820  for (it = rows.begin(); it != rows.end(); ++it)
     821    flagrowCol_.put(*it, rowflag);
    675822}
    676823
     
    793940  table_.keywordSet().get("FluxUnit", tmp);
    794941  oss << setw(15) << "Flux Unit:" << tmp << endl;
    795   Vector<Double> vec(moleculeTable_.getRestFrequencies());
     942  //Vector<Double> vec(moleculeTable_.getRestFrequencies());
     943  int nid = moleculeTable_.nrow();
     944  Bool firstline = True;
    796945  oss << setw(15) << "Rest Freqs:";
    797   if (vec.nelements() > 0) {
    798       oss << setprecision(10) << vec << " [Hz]" << endl;
    799   } else {
    800       oss << "none" << endl;
     946  for (int i=0; i<nid; i++) {
     947      Table t = table_(table_.col("MOLECULE_ID") == i);
     948      if (t.nrow() >  0) {
     949          Vector<Double> vec(moleculeTable_.getRestFrequency(i));
     950          if (vec.nelements() > 0) {
     951               if (firstline) {
     952                   oss << setprecision(10) << vec << " [Hz]" << endl;
     953                   firstline=False;
     954               }
     955               else{
     956                   oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
     957               }
     958          } else {
     959              oss << "none" << endl;
     960          }
     961      }
    801962  }
    802963
     
    9011062  const MDirection& md = getDirection(whichrow);
    9021063  const MEpoch& me = timeCol_(whichrow);
    903   Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     1064  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
     1065  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    9041066  return freqTable_.getSpectralCoordinate(md, mp, me, rf,
    9051067                                          mfreqidCol_(whichrow));
     
    9751137  const MDirection& md = getDirection(whichrow);
    9761138  const MEpoch& me = timeCol_(whichrow);
    977   const Double& rf = mmolidCol_(whichrow);
     1139  //const Double& rf = mmolidCol_(whichrow);
     1140  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
    9781141  SpectralCoordinate spc =
    9791142    freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
     
    9921155}
    9931156
    994 void Scantable::setRestFrequencies( double rf, const std::string& name,
     1157/**
     1158void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
    9951159                                          const std::string& unit )
     1160**/
     1161void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
     1162                                          const std::string& unit )
     1163
    9961164{
    9971165  ///@todo lookup in line table to fill in name and formattedname
    9981166  Unit u(unit);
    999   Quantum<Double> urf(rf, u);
    1000   uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
     1167  //Quantum<Double> urf(rf, u);
     1168  Quantum<Vector<Double> >urf(rf, u);
     1169  Vector<String> formattedname(0);
     1170  //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
     1171 
     1172  //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
     1173  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
    10011174  TableVector<uInt> tabvec(table_, "MOLECULE_ID");
    10021175  tabvec = id;
    10031176}
    10041177
    1005 void Scantable::setRestFrequencies( const std::string& name )
     1178/**
     1179void asap::Scantable::setRestFrequencies( const std::string& name )
    10061180{
    10071181  throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
     1182  ///@todo implement
     1183}
     1184**/
     1185void Scantable::setRestFrequencies( const vector<std::string>& name )
     1186{
     1187  throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
    10081188  ///@todo implement
    10091189}
     
    10441224void Scantable::addFit( const STFitEntry& fit, int row )
    10451225{
    1046   cout << mfitidCol_(uInt(row)) << endl;
     1226  //cout << mfitidCol_(uInt(row)) << endl;
     1227  LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
     1228  os << mfitidCol_(uInt(row)) << LogIO::POST ;
    10471229  uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
    10481230  mfitidCol_.put(uInt(row), id);
     
    10781260      Table subt = t( t.col("SCAN") == scanlist[i]+1 );
    10791261      if (subt.nrow()==0) {
    1080         cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
     1262        //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
     1263        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1264        os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
    10811265        ret = 1;
    10821266        break;
     
    10901274          Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
    10911275          if ( subt2.nrow() == 0) {
    1092             cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
     1276            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1277
     1278            //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
     1279            os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
    10931280            ret = 1;
    10941281            break;
     
    11001287          if (scan1seqn == 1 && scan2seqn == 2) {
    11011288            if (laston1 == laston2) {
    1102               cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1289              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1290              //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1291              os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
    11031292              i +=1;
    11041293            }
    11051294            else {
    1106               cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1295              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1296              //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
     1297              os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
    11071298            }
    11081299          }
    11091300          else if (scan1seqn==2 && scan2seqn == 1) {
    11101301            if (laston1 == laston2) {
    1111               cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
     1302              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1303              //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
     1304              os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
    11121305              ret = 1;
    11131306              break;
     
    11151308          }
    11161309          else {
    1117             cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
     1310            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1311            //cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
     1312            os<<LogIO::WARN<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
    11181313            ret = 1;
    11191314            break;
     
    11221317      }
    11231318      else {
    1124         cerr<<"The scan does not appear to be standard obsevation."<<endl;
     1319        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1320        //cerr<<"The scan does not appear to be standard obsevation."<<endl;
     1321        os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
    11251322      }
    11261323    //if ( i >= nscan ) break;
     
    11281325  }
    11291326  else {
    1130     cerr<<"No reference to GBT_GO table."<<endl;
     1327    LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
     1328    //cerr<<"No reference to GBT_GO table."<<endl;
     1329    os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
    11311330    ret = 1;
    11321331  }
     
    11421341}
    11431342
     1343void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
     1344  throw( casa::AipsError )
     1345{
     1346  // assumed that all rows have same nChan
     1347  Vector<Float> arr = specCol_( 0 ) ;
     1348  int nChan = arr.nelements() ;
     1349 
     1350  // if nmin < 0 or nmax < 0, nothing to do
     1351  if (  nmin < 0 ) {
     1352    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
     1353    }
     1354  if (  nmax < 0  ) {
     1355    throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
     1356  }
     1357 
     1358  // if nmin > nmax, exchange values
     1359  if ( nmin > nmax ) {
     1360    int tmp = nmax ;
     1361    nmax = nmin ;
     1362    nmin = tmp ;
     1363    LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
     1364    os << "Swap values. Applied range is ["
     1365       << nmin << ", " << nmax << "]" << LogIO::POST ;
     1366  }
     1367 
     1368  // if nmin exceeds nChan, nothing to do
     1369  if ( nmin >= nChan ) {
     1370    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
     1371  }
     1372 
     1373  // if nmax exceeds nChan, reset nmax to nChan
     1374  if ( nmax >= nChan ) {
     1375    if ( nmin == 0 ) {
     1376      // nothing to do
     1377      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
     1378      os << "Whole range is selected. Nothing to do." << LogIO::POST ;
     1379      return ;
     1380    }
     1381    else {
     1382      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
     1383      os << "Specified maximum exceeds nChan. Applied range is ["
     1384         << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
     1385      nmax = nChan - 1 ;
     1386    }
     1387  }
     1388 
     1389  // reshape specCol_ and flagCol_
     1390  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
     1391    reshapeSpectrum( nmin, nmax, irow ) ;
     1392  }
     1393
     1394  // update FREQUENCIES subtable
     1395  Double refpix ;
     1396  Double refval ;
     1397  Double increment ;
     1398  int freqnrow = freqTable_.table().nrow() ;
     1399  Vector<uInt> oldId( freqnrow ) ;
     1400  Vector<uInt> newId( freqnrow ) ;
     1401  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
     1402    freqTable_.getEntry( refpix, refval, increment, irow ) ;
     1403    /***
     1404     * need to shift refpix to nmin
     1405     * note that channel nmin in old index will be channel 0 in new one
     1406     ***/
     1407    refval = refval - ( refpix - nmin ) * increment ;
     1408    refpix = 0 ; 
     1409    freqTable_.setEntry( refpix, refval, increment, irow ) ;
     1410  }
     1411 
     1412  // update nchan
     1413  int newsize = nmax - nmin + 1 ;
     1414  table_.rwKeywordSet().define( "nChan", newsize ) ;
     1415 
     1416  // update bandwidth
     1417  // assumed all spectra in the scantable have same bandwidth
     1418  table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
     1419 
     1420  return ;
     1421}
     1422
     1423void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
     1424{
     1425  // reshape specCol_ and flagCol_
     1426  Vector<Float> oldspec = specCol_( irow ) ;
     1427  Vector<uChar> oldflag = flagsCol_( irow ) ;
     1428  uInt newsize = nmax - nmin + 1 ;
     1429  specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
     1430  flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
     1431
     1432  return ;
     1433}
     1434
     1435void asap::Scantable::regridChannel( int nChan, double dnu )
     1436{
     1437  LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
     1438  os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
     1439  // assumed that all rows have same nChan
     1440  Vector<Float> arr = specCol_( 0 ) ;
     1441  int oldsize = arr.nelements() ;
     1442
     1443  // if oldsize == nChan, nothing to do
     1444  if ( oldsize == nChan ) {
     1445    os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
     1446    return ;
     1447  }
     1448
     1449  // if oldChan < nChan, unphysical operation
     1450  if ( oldsize < nChan ) {
     1451    os << "Unphysical operation. Nothing to do." << LogIO::POST ;
     1452    return ;
     1453  }
     1454
     1455  // change channel number for specCol_ and flagCol_
     1456  Vector<Float> newspec( nChan, 0 ) ;
     1457  Vector<uChar> newflag( nChan, false ) ;
     1458  vector<string> coordinfo = getCoordInfo() ;
     1459  string oldinfo = coordinfo[0] ;
     1460  coordinfo[0] = "Hz" ;
     1461  setCoordInfo( coordinfo ) ;
     1462  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
     1463    regridChannel( nChan, dnu, irow ) ;
     1464  }
     1465  coordinfo[0] = oldinfo ;
     1466  setCoordInfo( coordinfo ) ;
     1467
     1468
     1469  // NOTE: this method does not update metadata such as
     1470  //       FREQUENCIES subtable, nChan, Bandwidth, etc.
     1471
     1472  return ;
     1473}
     1474
     1475void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
     1476{
     1477  // logging
     1478  //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
     1479  //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
     1480
     1481  Vector<Float> oldspec = specCol_( irow ) ;
     1482  Vector<uChar> oldflag = flagsCol_( irow ) ;
     1483  Vector<Float> newspec( nChan, 0 ) ;
     1484  Vector<uChar> newflag( nChan, false ) ;
     1485 
     1486  // regrid
     1487  vector<double> abcissa = getAbcissa( irow ) ;
     1488  int oldsize = abcissa.size() ;
     1489  double olddnu = abcissa[1] - abcissa[0] ;
     1490  //int refChan = 0 ;
     1491  //double frac = 0.0 ;
     1492  //double wedge = 0.0 ;
     1493  //double pile = 0.0 ;
     1494  int ichan = 0 ;
     1495  double wsum = 0.0 ;
     1496  Vector<Float> z( nChan ) ;
     1497  z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
     1498  for ( int ii = 1 ; ii < nChan ; ii++ )
     1499    z[ii] = z[ii-1] + dnu ;
     1500  Vector<Float> zi( nChan+1 ) ;
     1501  Vector<Float> yi( oldsize + 1 ) ;
     1502  zi[0] = z[0] - 0.5 * dnu ;
     1503  zi[1] = z[0] + 0.5 * dnu ;
     1504  for ( int ii = 2 ; ii < nChan ; ii++ )
     1505    zi[ii] = zi[ii-1] + dnu ;
     1506  zi[nChan] = z[nChan-1] + 0.5 * dnu ;
     1507  yi[0] = abcissa[0] - 0.5 * olddnu ;
     1508  yi[1] = abcissa[1] + 0.5 * olddnu ;
     1509  for ( int ii = 2 ; ii < oldsize ; ii++ )
     1510    yi[ii] = abcissa[ii-1] + olddnu ;
     1511  yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
     1512  if ( dnu > 0.0 ) {
     1513    for ( int ii = 0 ; ii < nChan ; ii++ ) {
     1514      double zl = zi[ii] ;
     1515      double zr = zi[ii+1] ;
     1516      for ( int j = ichan ; j < oldsize ; j++ ) {
     1517        double yl = yi[j] ;
     1518        double yr = yi[j+1] ;
     1519        if ( yl <= zl ) {
     1520          if ( yr <= zl ) {
     1521            continue ;
     1522          }
     1523          else if ( yr <= zr ) {
     1524            newspec[ii] += oldspec[j] * ( yr - zl ) ;
     1525            newflag[ii] = newflag[ii] || oldflag[j] ;
     1526            wsum += ( yr - zl ) ;
     1527          }
     1528          else {
     1529            newspec[ii] += oldspec[j] * dnu ;
     1530            newflag[ii] = newflag[ii] || oldflag[j] ;
     1531            wsum += dnu ;
     1532            ichan = j ;
     1533            break ;
     1534          }
     1535        }
     1536        else if ( yl < zr ) {
     1537          if ( yr <= zr ) {
     1538              newspec[ii] += oldspec[j] * ( yr - yl ) ;
     1539              newflag[ii] = newflag[ii] || oldflag[j] ;
     1540              wsum += ( yr - yl ) ;
     1541          }
     1542          else {
     1543            newspec[ii] += oldspec[j] * ( zr - yl ) ;
     1544            newflag[ii] = newflag[ii] || oldflag[j] ;
     1545            wsum += ( zr - yl ) ;
     1546            ichan = j ;
     1547            break ;
     1548          }
     1549        }
     1550        else {
     1551          ichan = j - 1 ;
     1552          break ;
     1553        }
     1554      }
     1555      newspec[ii] /= wsum ;
     1556      wsum = 0.0 ;
     1557    }
     1558  }
     1559  else if ( dnu < 0.0 ) {
     1560    for ( int ii = 0 ; ii < nChan ; ii++ ) {
     1561      double zl = zi[ii] ;
     1562      double zr = zi[ii+1] ;
     1563      for ( int j = ichan ; j < oldsize ; j++ ) {
     1564        double yl = yi[j] ;
     1565        double yr = yi[j+1] ;
     1566        if ( yl >= zl ) {
     1567          if ( yr >= zl ) {
     1568            continue ;
     1569          }
     1570          else if ( yr >= zr ) {
     1571            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
     1572            newflag[ii] = newflag[ii] || oldflag[j] ;
     1573            wsum += abs( yr - zl ) ;
     1574          }
     1575          else {
     1576            newspec[ii] += oldspec[j] * abs( dnu ) ;
     1577            newflag[ii] = newflag[ii] || oldflag[j] ;
     1578            wsum += abs( dnu ) ;
     1579            ichan = j ;
     1580            break ;
     1581          }
     1582        }
     1583        else if ( yl > zr ) {
     1584          if ( yr >= zr ) {
     1585            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
     1586            newflag[ii] = newflag[ii] || oldflag[j] ;
     1587            wsum += abs( yr - yl ) ;
     1588          }
     1589          else {
     1590            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
     1591            newflag[ii] = newflag[ii] || oldflag[j] ;
     1592            wsum += abs( zr - yl ) ;
     1593            ichan = j ;
     1594            break ;
     1595          }
     1596        }
     1597        else {
     1598          ichan = j - 1 ;
     1599          break ;
     1600        }
     1601      }
     1602      newspec[ii] /= wsum ;
     1603      wsum = 0.0 ;
     1604    }
     1605  }
     1606//    * ichan = 0
     1607//    ***/
     1608//   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
     1609//   pile += dnu ;
     1610//   wedge = olddnu * ( refChan + 1 ) ;
     1611//   while ( wedge < pile ) {
     1612//     newspec[0] += olddnu * oldspec[refChan] ;
     1613//     newflag[0] = newflag[0] || oldflag[refChan] ;
     1614//     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
     1615//     refChan++ ;
     1616//     wedge += olddnu ;
     1617//     wsum += olddnu ;
     1618//     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1619//   }
     1620//   frac = ( wedge - pile ) / olddnu ;
     1621//   wsum += ( 1.0 - frac ) * olddnu ;
     1622//   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1623//   newflag[0] = newflag[0] || oldflag[refChan] ;
     1624//   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
     1625//   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
     1626//   newspec[0] /= wsum ;
     1627//   //ofs << "newspec[0] = " << newspec[0] << endl ;
     1628//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1629
     1630//   /***
     1631//    * ichan = 1 - nChan-2
     1632//    ***/
     1633//   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
     1634//     pile += dnu ;
     1635//     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
     1636//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1637//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
     1638//     refChan++ ;
     1639//     wedge += olddnu ;
     1640//     wsum = frac * olddnu ;
     1641//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1642//     while ( wedge < pile ) {
     1643//       newspec[ichan] += olddnu * oldspec[refChan] ;
     1644//       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1645//       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
     1646//       refChan++ ;
     1647//       wedge += olddnu ;
     1648//       wsum += olddnu ;
     1649//       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1650//     }
     1651//     frac = ( wedge - pile ) / olddnu ;
     1652//     wsum += ( 1.0 - frac ) * olddnu ;
     1653//     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
     1654//     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
     1655//     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
     1656//     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1657//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
     1658//     newspec[ichan] /= wsum ;
     1659//     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
     1660//   }
     1661
     1662//   /***
     1663//    * ichan = nChan-1
     1664//    ***/
     1665//   // NOTE: Assumed that all spectra have the same bandwidth
     1666//   pile += dnu ;
     1667//   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
     1668//   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
     1669//   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1670//   refChan++ ;
     1671//   wedge += olddnu ;
     1672//   wsum = frac * olddnu ;
     1673//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1674//   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
     1675//     newspec[nChan-1] += olddnu * oldspec[jchan] ;
     1676//     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
     1677//     wsum += olddnu ;
     1678//     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
     1679//     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1680//   }
     1681//   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
     1682//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
     1683//   newspec[nChan-1] /= wsum ;
     1684//   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
     1685 
     1686//   specCol_.put( irow, newspec ) ;
     1687//   flagsCol_.put( irow, newflag ) ;
     1688
     1689//   // ofs.close() ;
     1690
     1691
     1692  return ;
     1693}
     1694
    11441695std::vector<float> Scantable::getWeather(int whichrow) const
    11451696{
     
    11541705
    11551706}
    1156  //namespace asap
     1707//namespace asap
  • branches/mergetest/src/Scantable.h

    r1730 r1789  
    2929
    3030#include <coordinates/Coordinates/SpectralCoordinate.h>
     31
     32#include <casa/Arrays/Vector.h>
     33#include <casa/Quanta/Quantum.h>
     34
     35#include <casa/Exceptions/Error.h>
    3136
    3237#include "Logger.h"
     
    226231
    227232  /**
     233   * Flag the data in a row-based manner. (CAS-1433 Wataru Kawasaki)
     234   * param[in] rows    list of row numbers to be flagged
     235   */
     236  void flagRow( const std::vector<casa::uInt>& rows = std::vector<casa::uInt>(), bool unflag=false);
     237
     238  /**
     239   * Get flagRow info at the specified row. If true, the whole data
     240   * at the row should be flagged.
     241   */
     242  bool getFlagRow(int whichrow) const
     243    { return (flagrowCol_(whichrow) > 0); }
     244
     245  /**
     246   * Flag the data outside a specified range (in a channel-based manner).
     247   * (CAS-1807 Wataru Kawasaki)
     248   */
     249  void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
     250
     251  /**
     252   * Return a list of booleans with the size of nchan for a specified row, to get info
     253   * about which channel is clipped.
     254   */
     255  std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag);
     256  void srchChannelsToClip(casa::uInt whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag,
     257                          casa::Vector<casa::uChar> flgs);
     258
     259  /**
    228260   * Return a list of row numbers with respect to the original table.
    229261   * @return a list of unsigned ints
     
    276308  std::vector<uint> getScanNos() const { return getNumbers(scanCol_); }
    277309  int getScan(int whichrow) const { return scanCol_(whichrow); }
     310
     311  //TT addition
     312  std::vector<uint> getMolNos() {return getNumbers(mmolidCol_); }
    278313
    279314  /**
     
    300335    { return azCol_(whichrow); }
    301336  float getParAngle(int whichrow) const
    302   { return focus().getParAngle(mfocusidCol_(whichrow)); }
     337    { return focus().getParAngle(mfocusidCol_(whichrow)); }
     338  int getTcalId(int whichrow) const
     339    { return mtcalidCol_(whichrow); }
    303340
    304341  std::string getSourceName(int whichrow) const
     
    353390  std::vector<double> getRestFrequencies() const
    354391    { return moleculeTable_.getRestFrequencies(); }
    355 
     392  std::vector<double> getRestFrequency(int id) const
     393    { return moleculeTable_.getRestFrequency(id); }
     394
     395  /**
    356396  void setRestFrequencies(double rf, const std::string& name = "",
    357397                          const std::string& = "Hz");
    358   void setRestFrequencies(const std::string& name);
     398  **/
     399  // Modified by Takeshi Nakazato 05/09/2008
     400  /***
     401  void setRestFrequencies(vector<double> rf, const vector<std::string>& name = "",
     402                          const std::string& = "Hz");
     403  ***/
     404  void setRestFrequencies(vector<double> rf,
     405                          const vector<std::string>& name = vector<std::string>(1,""),
     406                          const std::string& = "Hz");
     407
     408  //void setRestFrequencies(const std::string& name);
     409  void setRestFrequencies(const vector<std::string>& name);
    359410
    360411  void shift(int npix);
     
    415466    { focus().setParallactify(flag); }
    416467
     468  /**
     469   * Reshape spectrum
     470   * @param[in] nmin, nmax minimum and maximum channel
     471   * @param[in] irow       row number
     472   *
     473   * 30/07/2008 Takeshi Nakazato 
     474   **/
     475  void reshapeSpectrum( int nmin, int nmax ) throw( casa::AipsError );
     476  void reshapeSpectrum( int nmin, int nmax, int irow ) ;
     477
     478  /**
     479   * Change channel number under fixed bandwidth
     480   * @param[in] nchan, dnu new channel number and spectral resolution
     481   * @param[in] irow       row number
     482   *
     483   * 27/08/2008 Takeshi Nakazato
     484   **/
     485  void regridChannel( int nchan, double dnu ) ;
     486  void regridChannel( int nchan, double dnu, int irow ) ;
     487
     488 
    417489private:
    418490
     
    489561  casa::ScalarColumn<casa::Float> elCol_;
    490562  casa::ScalarColumn<casa::String> srcnCol_, fldnCol_;
    491   casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_;
     563  casa::ScalarColumn<casa::uInt> scanCol_, beamCol_, ifCol_, polCol_, cycleCol_, flagrowCol_;
    492564  casa::ScalarColumn<casa::Int> rbeamCol_, srctCol_;
    493565  casa::ArrayColumn<casa::Float> specCol_, tsysCol_;
     
    510582  void initFactories();
    511583
     584  /**
     585   * Add an auxiliary column to the main table and attach it to a
     586   * cached column. Use for adding new columns that the original asap2
     587   * tables do not have.
     588   * @param[in] col      reference to the cached column to be attached
     589   * @param[in] colName  column name in asap table
     590   * @param[in] defValue default value to fill in the column
     591   *
     592   * 25/10/2009 Wataru Kawasaki
     593   */
     594  template<class T, class T2> void attachAuxColumnDef(casa::ScalarColumn<T>&,
     595                                                       const casa::String&,
     596                                                       const T2&);
     597  template<class T, class T2> void attachAuxColumnDef(casa::ArrayColumn<T>&,
     598                                                      const casa::String&,
     599                                                      const casa::Array<T2>&);
    512600};
    513601
  • branches/mergetest/src/ScantableWrapper.h

    r1730 r1789  
    109109    { table_->flag(msk, unflag); }
    110110
     111  void flagRow(const std::vector<casa::uInt>& rows=std::vector<casa::uInt>(), bool unflag=false)
     112    { table_->flagRow(rows, unflag); }
     113
     114  bool getFlagRow(int whichrow=0) const
     115    { return table_->getFlagRow(whichrow); }
     116
     117  void clip(const casa::Float uthres, const casa::Float dthres, bool clipoutside=true, bool unflag=false)
     118    { table_->clip(uthres, dthres, clipoutside, unflag); }
     119
     120  std::vector<bool> getClipMask(int whichrow, const casa::Float uthres, const casa::Float dthres, bool clipoutside, bool unflag) const
     121    { return table_->getClipMask(whichrow, uthres, dthres, clipoutside, unflag); }
     122
    111123  std::string getSourceName(int whichrow=0) const
    112124    { return table_->getSourceName(whichrow); }
     
    134146  std::vector<uint> getScanNos() { return table_->getScanNos(); }
    135147  int getScan(int whichrow) const {return table_->getScan(whichrow);}
     148  std::vector<uint> getMolNos() { return table_->getMolNos();}
    136149
    137150  STSelector getSelection() const { return table_->getSelection(); }
     
    158171  { table_->shift(npix); }
    159172
     173/**
     174  commented out by TT
    160175  void setRestFrequencies(double rf, const std::string& name,
    161176                          const std::string& unit)
    162177    { table_->setRestFrequencies(rf, name, unit); }
     178**/
     179  void setRestFrequencies(vector<double> rf, const vector<std::string>& name,
     180                          const std::string& unit)
     181    { table_->setRestFrequencies(rf, name, unit); }
     182
    163183/*
    164184  void setRestFrequencies(const std::string& name) {
     
    167187*/
    168188
     189/*
    169190  std::vector<double> getRestFrequencies() const
    170191    { return table_->getRestFrequencies(); }
     192*/
     193  std::vector<double> getRestFrequency(int id) const
     194    { return table_->getRestFrequency(id); }
    171195
    172196  void setCoordInfo(std::vector<string> theinfo) {
     
    209233  int checkScanInfo(const vector<int>& scanlist) const
    210234    { return table_->checkScanInfo(scanlist); }
    211 
     235 
    212236  std::vector<double> getDirectionVector(int whichrow) const
    213237    { return table_->getDirectionVector(whichrow); }
     
    222246  std::vector<float> getWeather(int whichrow) const
    223247    { return table_->getWeather(whichrow); }
     248
     249  void reshapeSpectrum( int nmin, int nmax )
     250  { table_->reshapeSpectrum( nmin, nmax ); }
    224251
    225252private:
     
    229256} // namespace
    230257#endif
     258
  • branches/mergetest/src/python_STMath.cpp

    r1391 r1789  
    4949        .def("_averagebeams", &STMathWrapper::averageBeams)
    5050        .def("_unaryop", &STMathWrapper::unaryOperate)
     51        .def("_arrayop", &STMathWrapper::arrayOperate)
     52        //.def("_array2dop", &STMathWrapper::array2dOperate)
    5153        .def("_binaryop", &STMathWrapper::binaryOperate)
    5254        .def("_auto_quotient", &STMathWrapper::autoQuotient)
     
    5759        .def("_dofs", &STMathWrapper::dofs)
    5860        .def("_stats", &STMathWrapper::statistic)
     61        .def("_minmaxchan", &STMathWrapper::minMaxChan)
    5962        .def("_freqswitch", &STMathWrapper::freqSwitch)
    6063        .def("_bin", &STMathWrapper::bin)
     
    7376        .def("_mx_extract", &STMathWrapper::mxExtract)
    7477        .def("_lag_flag", &STMathWrapper::lagFlag)
     78        // testing average spectra with different channel/resolution
     79        .def("_new_average", &STMathWrapper::new_average)
     80        // cwcal
     81        .def("cwcal", &STMathWrapper::cwcal)
     82        .def("almacal", &STMathWrapper::almacal)
    7583          ;
    7684    };
  • branches/mergetest/src/python_STSelector.cpp

    r939 r1789  
    2929        .def("_getscans", &STSelector::getScans)
    3030        .def("_getcycles", &STSelector::getCycles)
     31        .def("_gettypes", &STSelector::getTypes)
    3132        .def("_gettaql", &STSelector::getTaQL)
    3233        .def("_getorder", &STSelector::getSortOrder)
     
    4142        .def("_settaql", &STSelector::setTaQL)
    4243        .def("_setorder", &STSelector::setSortOrder)
     44        .def("_setrows", &STSelector::setRows)
     45        .def("_settypes", &STSelector::setTypes)
    4346        .def("_empty", &STSelector::empty)
    4447      ;
  • branches/mergetest/src/python_Scantable.cpp

    r1730 r1789  
    5858    .def("getscannos", &ScantableWrapper::getScanNos)
    5959    .def("getcycle", &ScantableWrapper::getCycle)
     60    .def("getmolnos", &ScantableWrapper::getMolNos)
    6061    .def("nif", &ScantableWrapper::nif,
    6162         (boost::python::arg("scanno")=-1) )
     
    8788    .def("_getmask", &ScantableWrapper::getMask,
    8889         (boost::python::arg("whichrow")=0) )
     90    .def("_getclipmask", &ScantableWrapper::getClipMask,
     91         (boost::python::arg("whichrow")=0) )
    8992    .def("_gettsys", &ScantableWrapper::getTsys)
    9093    .def("_getsourcename", &ScantableWrapper::getSourceName,
     
    103106         (boost::python::arg("whichrow")=0) )
    104107    .def("get_antennaname", &ScantableWrapper::getAntennaName)
    105     .def("_flag", &ScantableWrapper::flag,
    106                   (boost::python::arg("unflag") = false) )
     108    .def("_flag", &ScantableWrapper::flag)
     109    .def("_flag_row", &ScantableWrapper::flagRow)
     110    .def("_getflagrow", &ScantableWrapper::getFlagRow,
     111         (boost::python::arg("whichrow")=0) )
     112    .def("_clip", &ScantableWrapper::clip,
     113         (boost::python::arg("clipoutside")=true,
     114          boost::python::arg("unflag")=false) )
    107115    .def("_save",  &ScantableWrapper::makePersistent)
    108116    .def("_summary",  &ScantableWrapper::summary,
    109117         (boost::python::arg("verbose")=true) )
    110     .def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
     118    //.def("_getrestfreqs",  &ScantableWrapper::getRestFrequencies)
     119    .def("_getrestfreqs",  &ScantableWrapper::getRestFrequency)
    111120    .def("_setrestfreqs",  &ScantableWrapper::setRestFrequencies)
    112121    .def("shift_refpix", &ScantableWrapper::shift)
     
    127136    .def("get_coordinate", &ScantableWrapper::getCoordinate)
    128137    .def("_get_weather", &ScantableWrapper::getWeather)
     138    .def("_reshape", &ScantableWrapper::reshapeSpectrum,
     139         (boost::python::arg("nmin")=-1,
     140          boost::python::arg("nmax")=-1) )
    129141  ;
    130142};
  • branches/mergetest/src/python_asap.cpp

    r1715 r1789  
    7777  asap::python::python_STCoordinate();
    7878  asap::python::python_STAtmosphere();
     79  asap::python::python_SrcType();
    7980
    8081#ifndef HAVE_LIBPYRAP
  • branches/mergetest/src/python_asap.h

    r1715 r1789  
    4949    void python_STCoordinate();
    5050    void python_STAtmosphere();
     51    void python_SrcType();
    5152
    5253  } // python
Note: See TracChangeset for help on using the changeset viewer.