source: trunk/src/Scantable.cpp @ 2463

Last change on this file since 2463 was 2463, checked in by Kana Sugimoto, 12 years ago

New Development: No (a bug fix)

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): scantable._regrid_specchan

Description: a bug fix in Scantable::regridSpecChannel.

Fixed the calculation of regridded reference channel.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 116.1 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13
14#include <atnf/PKSIO/SrcType.h>
15
16#include <casa/aips.h>
17#include <casa/iomanip.h>
18#include <casa/iostream.h>
19#include <casa/OS/File.h>
20#include <casa/OS/Path.h>
21#include <casa/Arrays/Array.h>
22#include <casa/Arrays/ArrayAccessor.h>
23#include <casa/Arrays/ArrayLogical.h>
24#include <casa/Arrays/ArrayMath.h>
25#include <casa/Arrays/MaskArrMath.h>
26#include <casa/Arrays/Slice.h>
27#include <casa/Arrays/Vector.h>
28#include <casa/Arrays/VectorSTLIterator.h>
29#include <casa/BasicMath/Math.h>
30#include <casa/BasicSL/Constants.h>
31#include <casa/Containers/RecordField.h>
32#include <casa/Logging/LogIO.h>
33#include <casa/Quanta/MVAngle.h>
34#include <casa/Quanta/MVTime.h>
35#include <casa/Utilities/GenSort.h>
36
37#include <coordinates/Coordinates/CoordinateUtil.h>
38
39// needed to avoid error in .tcc
40#include <measures/Measures/MCDirection.h>
41//
42#include <measures/Measures/MDirection.h>
43#include <measures/Measures/MEpoch.h>
44#include <measures/Measures/MFrequency.h>
45#include <measures/Measures/MeasRef.h>
46#include <measures/Measures/MeasTable.h>
47#include <measures/TableMeasures/ScalarMeasColumn.h>
48#include <measures/TableMeasures/TableMeasDesc.h>
49#include <measures/TableMeasures/TableMeasRefDesc.h>
50#include <measures/TableMeasures/TableMeasValueDesc.h>
51
52#include <tables/Tables/ArrColDesc.h>
53#include <tables/Tables/ExprNode.h>
54#include <tables/Tables/ScaColDesc.h>
55#include <tables/Tables/SetupNewTab.h>
56#include <tables/Tables/TableCopy.h>
57#include <tables/Tables/TableDesc.h>
58#include <tables/Tables/TableIter.h>
59#include <tables/Tables/TableParse.h>
60#include <tables/Tables/TableRecord.h>
61#include <tables/Tables/TableRow.h>
62#include <tables/Tables/TableVector.h>
63
64#include "MathUtils.h"
65#include "STAttr.h"
66#include "STLineFinder.h"
67#include "STPolCircular.h"
68#include "STPolLinear.h"
69#include "STPolStokes.h"
70#include "STUpgrade.h"
71#include "Scantable.h"
72
73#define debug 1
74
75using namespace casa;
76
77namespace asap {
78
79std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
80
81void Scantable::initFactories() {
82  if ( factories_.empty() ) {
83    Scantable::factories_["linear"] = &STPolLinear::myFactory;
84    Scantable::factories_["circular"] = &STPolCircular::myFactory;
85    Scantable::factories_["stokes"] = &STPolStokes::myFactory;
86  }
87}
88
89Scantable::Scantable(Table::TableType ttype) :
90  type_(ttype)
91{
92  initFactories();
93  setupMainTable();
94  freqTable_ = STFrequencies(*this);
95  table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
96  weatherTable_ = STWeather(*this);
97  table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
98  focusTable_ = STFocus(*this);
99  table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
100  tcalTable_ = STTcal(*this);
101  table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
102  moleculeTable_ = STMolecules(*this);
103  table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
104  historyTable_ = STHistory(*this);
105  table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
106  fitTable_ = STFit(*this);
107  table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
108  table_.tableInfo().setType( "Scantable" ) ;
109  originalTable_ = table_;
110  attach();
111}
112
113Scantable::Scantable(const std::string& name, Table::TableType ttype) :
114  type_(ttype)
115{
116  initFactories();
117
118  Table tab(name, Table::Update);
119  uInt version = tab.keywordSet().asuInt("VERSION");
120  if (version != version_) {
121      STUpgrade upgrader(version_);
122      LogIO os( LogOrigin( "Scantable" ) ) ;
123      os << LogIO::WARN
124         << name << " data format version " << version
125         << " is deprecated" << endl
126         << "Running upgrade."<< endl 
127         << LogIO::POST ; 
128      std::string outname = upgrader.upgrade(name);
129      if ( outname != name ) {
130        os << LogIO::WARN
131           << "Data will be loaded from " << outname << " instead of "
132           << name << LogIO::POST ;
133        tab = Table(outname, Table::Update ) ;
134      }
135  }
136  if ( type_ == Table::Memory ) {
137    table_ = tab.copyToMemoryTable(generateName());
138  } else {
139    table_ = tab;
140  }
141  table_.tableInfo().setType( "Scantable" ) ;
142
143  attachSubtables();
144  originalTable_ = table_;
145  attach();
146}
147/*
148Scantable::Scantable(const std::string& name, Table::TableType ttype) :
149  type_(ttype)
150{
151  initFactories();
152  Table tab(name, Table::Update);
153  uInt version = tab.keywordSet().asuInt("VERSION");
154  if (version != version_) {
155    throw(AipsError("Unsupported version of ASAP file."));
156  }
157  if ( type_ == Table::Memory ) {
158    table_ = tab.copyToMemoryTable(generateName());
159  } else {
160    table_ = tab;
161  }
162
163  attachSubtables();
164  originalTable_ = table_;
165  attach();
166}
167*/
168
169Scantable::Scantable( const Scantable& other, bool clear ):
170  Logger()
171{
172  // with or without data
173  String newname = String(generateName());
174  type_ = other.table_.tableType();
175  if ( other.table_.tableType() == Table::Memory ) {
176      if ( clear ) {
177        table_ = TableCopy::makeEmptyMemoryTable(newname,
178                                                 other.table_, True);
179      } else
180        table_ = other.table_.copyToMemoryTable(newname);
181  } else {
182      other.table_.deepCopy(newname, Table::New, False,
183                            other.table_.endianFormat(),
184                            Bool(clear));
185      table_ = Table(newname, Table::Update);
186      table_.markForDelete();
187  }
188  table_.tableInfo().setType( "Scantable" ) ;
189  /// @todo reindex SCANNO, recompute nbeam, nif, npol
190  if ( clear ) copySubtables(other);
191  attachSubtables();
192  originalTable_ = table_;
193  attach();
194}
195
196void Scantable::copySubtables(const Scantable& other) {
197  Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
198  TableCopy::copyRows(t, other.freqTable_.table());
199  t = table_.rwKeywordSet().asTable("FOCUS");
200  TableCopy::copyRows(t, other.focusTable_.table());
201  t = table_.rwKeywordSet().asTable("WEATHER");
202  TableCopy::copyRows(t, other.weatherTable_.table());
203  t = table_.rwKeywordSet().asTable("TCAL");
204  TableCopy::copyRows(t, other.tcalTable_.table());
205  t = table_.rwKeywordSet().asTable("MOLECULES");
206  TableCopy::copyRows(t, other.moleculeTable_.table());
207  t = table_.rwKeywordSet().asTable("HISTORY");
208  TableCopy::copyRows(t, other.historyTable_.table());
209  t = table_.rwKeywordSet().asTable("FIT");
210  TableCopy::copyRows(t, other.fitTable_.table());
211}
212
213void Scantable::attachSubtables()
214{
215  freqTable_ = STFrequencies(table_);
216  focusTable_ = STFocus(table_);
217  weatherTable_ = STWeather(table_);
218  tcalTable_ = STTcal(table_);
219  moleculeTable_ = STMolecules(table_);
220  historyTable_ = STHistory(table_);
221  fitTable_ = STFit(table_);
222}
223
224Scantable::~Scantable()
225{
226}
227
228void Scantable::setupMainTable()
229{
230  TableDesc td("", "1", TableDesc::Scratch);
231  td.comment() = "An ASAP Scantable";
232  td.rwKeywordSet().define("VERSION", uInt(version_));
233
234  // n Cycles
235  td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
236  // new index every nBeam x nIF x nPol
237  td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
238
239  td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
240  td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
241  // linear, circular, stokes
242  td.rwKeywordSet().define("POLTYPE", String("linear"));
243  td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
244
245  td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
246  td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
247
248  ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
249  refbeamnoColumn.setDefault(Int(-1));
250  td.addColumn(refbeamnoColumn);
251
252  ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
253  flagrowColumn.setDefault(uInt(0));
254  td.addColumn(flagrowColumn);
255
256  td.addColumn(ScalarColumnDesc<Double>("TIME"));
257  TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
258  TableMeasValueDesc measVal(td, "TIME");
259  TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
260  mepochCol.write(td);
261
262  td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
263
264  td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
265  // Type of source (on=0, off=1, other=-1)
266  ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
267  stypeColumn.setDefault(Int(-1));
268  td.addColumn(stypeColumn);
269  td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
270
271  //The actual Data Vectors
272  td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
273  td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
274  td.addColumn(ArrayColumnDesc<Float>("TSYS"));
275
276  td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
277                                       IPosition(1,2),
278                                       ColumnDesc::Direct));
279  TableMeasRefDesc mdirRef(MDirection::J2000); // default
280  TableMeasValueDesc tmvdMDir(td, "DIRECTION");
281  // the TableMeasDesc gives the column a type
282  TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
283  // a uder set table type e.g. GALCTIC, B1950 ...
284  td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
285  // writing create the measure column
286  mdirCol.write(td);
287  td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
288  td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
289  td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
290
291  td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
292  ScalarColumnDesc<Int> fitColumn("FIT_ID");
293  fitColumn.setDefault(Int(-1));
294  td.addColumn(fitColumn);
295
296  td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
297  td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
298
299  // columns which just get dragged along, as they aren't used in asap
300  td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
301  td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
302  td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
303  td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
304
305  td.rwKeywordSet().define("OBSMODE", String(""));
306
307  // Now create Table SetUp from the description.
308  SetupNewTable aNewTab(generateName(), td, Table::Scratch);
309  table_ = Table(aNewTab, type_, 0);
310  originalTable_ = table_;
311}
312
313void Scantable::attach()
314{
315  timeCol_.attach(table_, "TIME");
316  srcnCol_.attach(table_, "SRCNAME");
317  srctCol_.attach(table_, "SRCTYPE");
318  specCol_.attach(table_, "SPECTRA");
319  flagsCol_.attach(table_, "FLAGTRA");
320  tsysCol_.attach(table_, "TSYS");
321  cycleCol_.attach(table_,"CYCLENO");
322  scanCol_.attach(table_, "SCANNO");
323  beamCol_.attach(table_, "BEAMNO");
324  ifCol_.attach(table_, "IFNO");
325  polCol_.attach(table_, "POLNO");
326  integrCol_.attach(table_, "INTERVAL");
327  azCol_.attach(table_, "AZIMUTH");
328  elCol_.attach(table_, "ELEVATION");
329  dirCol_.attach(table_, "DIRECTION");
330  fldnCol_.attach(table_, "FIELDNAME");
331  rbeamCol_.attach(table_, "REFBEAMNO");
332
333  mweatheridCol_.attach(table_,"WEATHER_ID");
334  mfitidCol_.attach(table_,"FIT_ID");
335  mfreqidCol_.attach(table_, "FREQ_ID");
336  mtcalidCol_.attach(table_, "TCAL_ID");
337  mfocusidCol_.attach(table_, "FOCUS_ID");
338  mmolidCol_.attach(table_, "MOLECULE_ID");
339
340  //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
341  attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
342
343}
344
345template<class T, class T2>
346void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
347                                   const String& colName,
348                                   const T2& defValue)
349{
350  try {
351    col.attach(table_, colName);
352  } catch (TableError& err) {
353    String errMesg = err.getMesg();
354    if (errMesg == "Table column " + colName + " is unknown") {
355      table_.addColumn(ScalarColumnDesc<T>(colName));
356      col.attach(table_, colName);
357      col.fillColumn(static_cast<T>(defValue));
358    } else {
359      throw;
360    }
361  } catch (...) {
362    throw;
363  }
364}
365
366template<class T, class T2>
367void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
368                                   const String& colName,
369                                   const Array<T2>& defValue)
370{
371  try {
372    col.attach(table_, colName);
373  } catch (TableError& err) {
374    String errMesg = err.getMesg();
375    if (errMesg == "Table column " + colName + " is unknown") {
376      table_.addColumn(ArrayColumnDesc<T>(colName));
377      col.attach(table_, colName);
378
379      int size = 0;
380      ArrayIterator<T2>& it = defValue.begin();
381      while (it != defValue.end()) {
382        ++size;
383        ++it;
384      }
385      IPosition ip(1, size);
386      Array<T>& arr(ip);
387      for (int i = 0; i < size; ++i)
388        arr[i] = static_cast<T>(defValue[i]);
389
390      col.fillColumn(arr);
391    } else {
392      throw;
393    }
394  } catch (...) {
395    throw;
396  }
397}
398
399void Scantable::setHeader(const STHeader& sdh)
400{
401  table_.rwKeywordSet().define("nIF", sdh.nif);
402  table_.rwKeywordSet().define("nBeam", sdh.nbeam);
403  table_.rwKeywordSet().define("nPol", sdh.npol);
404  table_.rwKeywordSet().define("nChan", sdh.nchan);
405  table_.rwKeywordSet().define("Observer", sdh.observer);
406  table_.rwKeywordSet().define("Project", sdh.project);
407  table_.rwKeywordSet().define("Obstype", sdh.obstype);
408  table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
409  table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
410  table_.rwKeywordSet().define("Equinox", sdh.equinox);
411  table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
412  table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
413  table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
414  table_.rwKeywordSet().define("UTC", sdh.utc);
415  table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
416  table_.rwKeywordSet().define("Epoch", sdh.epoch);
417  table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
418}
419
420STHeader Scantable::getHeader() const
421{
422  STHeader sdh;
423  table_.keywordSet().get("nBeam",sdh.nbeam);
424  table_.keywordSet().get("nIF",sdh.nif);
425  table_.keywordSet().get("nPol",sdh.npol);
426  table_.keywordSet().get("nChan",sdh.nchan);
427  table_.keywordSet().get("Observer", sdh.observer);
428  table_.keywordSet().get("Project", sdh.project);
429  table_.keywordSet().get("Obstype", sdh.obstype);
430  table_.keywordSet().get("AntennaName", sdh.antennaname);
431  table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
432  table_.keywordSet().get("Equinox", sdh.equinox);
433  table_.keywordSet().get("FreqRefFrame", sdh.freqref);
434  table_.keywordSet().get("FreqRefVal", sdh.reffreq);
435  table_.keywordSet().get("Bandwidth", sdh.bandwidth);
436  table_.keywordSet().get("UTC", sdh.utc);
437  table_.keywordSet().get("FluxUnit", sdh.fluxunit);
438  table_.keywordSet().get("Epoch", sdh.epoch);
439  table_.keywordSet().get("POLTYPE", sdh.poltype);
440  return sdh;
441}
442
443void Scantable::setSourceType( int stype )
444{
445  if ( stype < 0 || stype > 1 )
446    throw(AipsError("Illegal sourcetype."));
447  TableVector<Int> tabvec(table_, "SRCTYPE");
448  tabvec = Int(stype);
449}
450
451bool Scantable::conformant( const Scantable& other )
452{
453  return this->getHeader().conformant(other.getHeader());
454}
455
456
457
458std::string Scantable::formatSec(Double x) const
459{
460  Double xcop = x;
461  MVTime mvt(xcop/24./3600.);  // make days
462
463  if (x < 59.95)
464    return  String("      ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
465  else if (x < 3599.95)
466    return String("   ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
467  else {
468    ostringstream oss;
469    oss << setw(2) << std::right << setprecision(1) << mvt.hour();
470    oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
471    return String(oss);
472  }
473};
474
475std::string Scantable::formatDirection(const MDirection& md) const
476{
477  Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
478  Int prec = 7;
479
480  MVAngle mvLon(t[0]);
481  String sLon = mvLon.string(MVAngle::TIME,prec);
482  uInt tp = md.getRef().getType();
483  if (tp == MDirection::GALACTIC ||
484      tp == MDirection::SUPERGAL ) {
485    sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
486  }
487  MVAngle mvLat(t[1]);
488  String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
489  return sLon + String(" ") + sLat;
490}
491
492
493std::string Scantable::getFluxUnit() const
494{
495  return table_.keywordSet().asString("FluxUnit");
496}
497
498void Scantable::setFluxUnit(const std::string& unit)
499{
500  String tmp(unit);
501  Unit tU(tmp);
502  if (tU==Unit("K") || tU==Unit("Jy")) {
503     table_.rwKeywordSet().define(String("FluxUnit"), tmp);
504  } else {
505     throw AipsError("Illegal unit - must be compatible with Jy or K");
506  }
507}
508
509void Scantable::setInstrument(const std::string& name)
510{
511  bool throwIt = true;
512  // create an Instrument to see if this is valid
513  STAttr::convertInstrument(name, throwIt);
514  String nameU(name);
515  nameU.upcase();
516  table_.rwKeywordSet().define(String("AntennaName"), nameU);
517}
518
519void Scantable::setFeedType(const std::string& feedtype)
520{
521  if ( Scantable::factories_.find(feedtype) ==  Scantable::factories_.end() ) {
522    std::string msg = "Illegal feed type "+ feedtype;
523    throw(casa::AipsError(msg));
524  }
525  table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
526}
527
528MPosition Scantable::getAntennaPosition() const
529{
530  Vector<Double> antpos;
531  table_.keywordSet().get("AntennaPosition", antpos);
532  MVPosition mvpos(antpos(0),antpos(1),antpos(2));
533  return MPosition(mvpos);
534}
535
536void Scantable::makePersistent(const std::string& filename)
537{
538  String inname(filename);
539  Path path(inname);
540  /// @todo reindex SCANNO, recompute nbeam, nif, npol
541  inname = path.expandedName();
542  // 2011/03/04 TN
543  // We can comment out this workaround since the essential bug is
544  // fixed in casacore (r20889 in google code).
545  table_.deepCopy(inname, Table::New);
546//   // WORKAROUND !!! for Table bug
547//   // Remove when fixed in casacore
548//   if ( table_.tableType() == Table::Memory  && !selector_.empty() ) {
549//     Table tab = table_.copyToMemoryTable(generateName());
550//     tab.deepCopy(inname, Table::New);
551//     tab.markForDelete();
552//
553//   } else {
554//     table_.deepCopy(inname, Table::New);
555//   }
556}
557
558int Scantable::nbeam( int scanno ) const
559{
560  if ( scanno < 0 ) {
561    Int n;
562    table_.keywordSet().get("nBeam",n);
563    return int(n);
564  } else {
565    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
566    Table t = table_(table_.col("SCANNO") == scanno);
567    ROTableRow row(t);
568    const TableRecord& rec = row.get(0);
569    Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
570                    && t.col("POLNO") == Int(rec.asuInt("POLNO"))
571                    && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
572    ROTableVector<uInt> v(subt, "BEAMNO");
573    return int(v.nelements());
574  }
575  return 0;
576}
577
578int Scantable::nif( int scanno ) const
579{
580  if ( scanno < 0 ) {
581    Int n;
582    table_.keywordSet().get("nIF",n);
583    return int(n);
584  } else {
585    // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
586    Table t = table_(table_.col("SCANNO") == scanno);
587    ROTableRow row(t);
588    const TableRecord& rec = row.get(0);
589    Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
590                    && t.col("POLNO") == Int(rec.asuInt("POLNO"))
591                    && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
592    if ( subt.nrow() == 0 ) return 0;
593    ROTableVector<uInt> v(subt, "IFNO");
594    return int(v.nelements());
595  }
596  return 0;
597}
598
599int Scantable::npol( int scanno ) const
600{
601  if ( scanno < 0 ) {
602    Int n;
603    table_.keywordSet().get("nPol",n);
604    return n;
605  } else {
606    // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
607    Table t = table_(table_.col("SCANNO") == scanno);
608    ROTableRow row(t);
609    const TableRecord& rec = row.get(0);
610    Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
611                    && t.col("IFNO") == Int(rec.asuInt("IFNO"))
612                    && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
613    if ( subt.nrow() == 0 ) return 0;
614    ROTableVector<uInt> v(subt, "POLNO");
615    return int(v.nelements());
616  }
617  return 0;
618}
619
620int Scantable::ncycle( int scanno ) const
621{
622  if ( scanno < 0 ) {
623    Block<String> cols(2);
624    cols[0] = "SCANNO";
625    cols[1] = "CYCLENO";
626    TableIterator it(table_, cols);
627    int n = 0;
628    while ( !it.pastEnd() ) {
629      ++n;
630      ++it;
631    }
632    return n;
633  } else {
634    Table t = table_(table_.col("SCANNO") == scanno);
635    ROTableRow row(t);
636    const TableRecord& rec = row.get(0);
637    Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
638                    && t.col("POLNO") == Int(rec.asuInt("POLNO"))
639                    && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
640    if ( subt.nrow() == 0 ) return 0;
641    return int(subt.nrow());
642  }
643  return 0;
644}
645
646
647int Scantable::nrow( int scanno ) const
648{
649  return int(table_.nrow());
650}
651
652int Scantable::nchan( int ifno ) const
653{
654  if ( ifno < 0 ) {
655    Int n;
656    table_.keywordSet().get("nChan",n);
657    return int(n);
658  } else {
659    // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
660    // vary with these
661    Table t = table_(table_.col("IFNO") == ifno, 1);
662    if ( t.nrow() == 0 ) return 0;
663    ROArrayColumn<Float> v(t, "SPECTRA");
664    return v.shape(0)(0);
665  }
666  return 0;
667}
668
669int Scantable::nscan() const {
670  Vector<uInt> scannos(scanCol_.getColumn());
671  uInt nout = genSort( scannos, Sort::Ascending,
672                       Sort::QuickSort|Sort::NoDuplicates );
673  return int(nout);
674}
675
676int Scantable::getChannels(int whichrow) const
677{
678  return specCol_.shape(whichrow)(0);
679}
680
681int Scantable::getBeam(int whichrow) const
682{
683  return beamCol_(whichrow);
684}
685
686std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
687{
688  Vector<uInt> nos(col.getColumn());
689  uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
690  nos.resize(n, True);
691  std::vector<uint> stlout;
692  nos.tovector(stlout);
693  return stlout;
694}
695
696int Scantable::getIF(int whichrow) const
697{
698  return ifCol_(whichrow);
699}
700
701int Scantable::getPol(int whichrow) const
702{
703  return polCol_(whichrow);
704}
705
706std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
707{
708  return formatTime(me, showdate, 0);
709}
710
711std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
712{
713  MVTime mvt(me.getValue());
714  if (showdate)
715    //mvt.setFormat(MVTime::YMD);
716    mvt.setFormat(MVTime::YMD, prec);
717  else
718    //mvt.setFormat(MVTime::TIME);
719    mvt.setFormat(MVTime::TIME, prec);
720  ostringstream oss;
721  oss << mvt;
722  return String(oss);
723}
724
725void Scantable::calculateAZEL()
726{
727  MPosition mp = getAntennaPosition();
728  MEpoch::ROScalarColumn timeCol(table_, "TIME");
729  ostringstream oss;
730  oss << "Computed azimuth/elevation using " << endl
731      << mp << endl;
732  for (Int i=0; i<nrow(); ++i) {
733    MEpoch me = timeCol(i);
734    MDirection md = getDirection(i);
735    oss  << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
736         << endl << "     => ";
737    MeasFrame frame(mp, me);
738    Vector<Double> azel =
739        MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
740                                                frame)
741                            )().getAngle("rad").getValue();
742    azCol_.put(i,Float(azel[0]));
743    elCol_.put(i,Float(azel[1]));
744    oss << "azel: " << azel[0]/C::pi*180.0 << " "
745        << azel[1]/C::pi*180.0 << " (deg)" << endl;
746  }
747  pushLog(String(oss));
748}
749
750void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
751{
752  for (uInt i=0; i<table_.nrow(); ++i) {
753    Vector<uChar> flgs = flagsCol_(i);
754    srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
755    flagsCol_.put(i, flgs);
756  }
757}
758
759std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
760{
761  Vector<uChar> flags;
762  flagsCol_.get(uInt(whichrow), flags);
763  srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
764  Vector<Bool> bflag(flags.shape());
765  convertArray(bflag, flags);
766  //bflag = !bflag;
767
768  std::vector<bool> mask;
769  bflag.tovector(mask);
770  return mask;
771}
772
773void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
774                                   Vector<uChar> flgs)
775{
776    Vector<Float> spcs = specCol_(whichrow);
777    uInt nchannel = spcs.nelements();
778    if (spcs.nelements() != nchannel) {
779      throw(AipsError("Data has incorrect number of channels"));
780    }
781    uChar userflag = 1 << 7;
782    if (unflag) {
783      userflag = 0 << 7;
784    }
785    if (clipoutside) {
786      for (uInt j = 0; j < nchannel; ++j) {
787        Float spc = spcs(j);
788        if ((spc >= uthres) || (spc <= dthres)) {
789          flgs(j) = userflag;
790        }
791      }
792    } else {
793      for (uInt j = 0; j < nchannel; ++j) {
794        Float spc = spcs(j);
795        if ((spc < uthres) && (spc > dthres)) {
796          flgs(j) = userflag;
797        }
798      }
799    }
800}
801
802
803void Scantable::flag( int whichrow, const std::vector<bool>& msk, bool unflag ) {
804  std::vector<bool>::const_iterator it;
805  uInt ntrue = 0;
806  if (whichrow >= int(table_.nrow()) ) {
807    throw(AipsError("Invalid row number"));
808  }
809  for (it = msk.begin(); it != msk.end(); ++it) {
810    if ( *it ) {
811      ntrue++;
812    }
813  }
814  //if ( selector_.empty()  && (msk.size() == 0 || msk.size() == ntrue) )
815  if ( whichrow == -1 && !unflag && selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
816    throw(AipsError("Trying to flag whole scantable."));
817  uChar userflag = 1 << 7;
818  if ( unflag ) {
819    userflag = 0 << 7;
820  }
821  if (whichrow > -1 ) {
822    applyChanFlag(uInt(whichrow), msk, userflag);
823  } else {
824    for ( uInt i=0; i<table_.nrow(); ++i) {
825      applyChanFlag(i, msk, userflag);
826    }
827  }
828}
829
830void Scantable::applyChanFlag( uInt whichrow, const std::vector<bool>& msk, uChar flagval )
831{
832  if (whichrow >= table_.nrow() ) {
833    throw( casa::indexError<int>( whichrow, "asap::Scantable::applyChanFlag: Invalid row number" ) );
834  }
835  Vector<uChar> flgs = flagsCol_(whichrow);
836  if ( msk.size() == 0 ) {
837    flgs = flagval;
838    flagsCol_.put(whichrow, flgs);
839    return;
840  }
841  if ( int(msk.size()) != nchan( getIF(whichrow) ) ) {
842    throw(AipsError("Mask has incorrect number of channels."));
843  }
844  if ( flgs.nelements() != msk.size() ) {
845    throw(AipsError("Mask has incorrect number of channels."
846                    " Probably varying with IF. Please flag per IF"));
847  }
848  std::vector<bool>::const_iterator it;
849  uInt j = 0;
850  for (it = msk.begin(); it != msk.end(); ++it) {
851    if ( *it ) {
852      flgs(j) = flagval;
853    }
854    ++j;
855  }
856  flagsCol_.put(whichrow, flgs);
857}
858
859void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
860{
861  if ( selector_.empty() && (rows.size() == table_.nrow()) )
862    throw(AipsError("Trying to flag whole scantable."));
863
864  uInt rowflag = (unflag ? 0 : 1);
865  std::vector<uInt>::const_iterator it;
866  for (it = rows.begin(); it != rows.end(); ++it)
867    flagrowCol_.put(*it, rowflag);
868}
869
870std::vector<bool> Scantable::getMask(int whichrow) const
871{
872  Vector<uChar> flags;
873  flagsCol_.get(uInt(whichrow), flags);
874  Vector<Bool> bflag(flags.shape());
875  convertArray(bflag, flags);
876  bflag = !bflag;
877  std::vector<bool> mask;
878  bflag.tovector(mask);
879  return mask;
880}
881
882std::vector<float> Scantable::getSpectrum( int whichrow,
883                                           const std::string& poltype ) const
884{
885  String ptype = poltype;
886  if (poltype == "" ) ptype = getPolType();
887  if ( whichrow  < 0 || whichrow >= nrow() )
888    throw(AipsError("Illegal row number."));
889  std::vector<float> out;
890  Vector<Float> arr;
891  uInt requestedpol = polCol_(whichrow);
892  String basetype = getPolType();
893  if ( ptype == basetype ) {
894    specCol_.get(whichrow, arr);
895  } else {
896    CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
897                                               basetype));
898    uInt row = uInt(whichrow);
899    stpol->setSpectra(getPolMatrix(row));
900    Float fang,fhand;
901    fang = focusTable_.getTotalAngle(mfocusidCol_(row));
902    fhand = focusTable_.getFeedHand(mfocusidCol_(row));
903    stpol->setPhaseCorrections(fang, fhand);
904    arr = stpol->getSpectrum(requestedpol, ptype);
905  }
906  if ( arr.nelements() == 0 )
907    pushLog("Not enough polarisations present to do the conversion.");
908  arr.tovector(out);
909  return out;
910}
911
912void Scantable::setSpectrum( const std::vector<float>& spec,
913                                   int whichrow )
914{
915  Vector<Float> spectrum(spec);
916  Vector<Float> arr;
917  specCol_.get(whichrow, arr);
918  if ( spectrum.nelements() != arr.nelements() )
919    throw AipsError("The spectrum has incorrect number of channels.");
920  specCol_.put(whichrow, spectrum);
921}
922
923
924String Scantable::generateName()
925{
926  return (File::newUniqueName("./","temp")).baseName();
927}
928
929const casa::Table& Scantable::table( ) const
930{
931  return table_;
932}
933
934casa::Table& Scantable::table( )
935{
936  return table_;
937}
938
939std::string Scantable::getPolType() const
940{
941  return table_.keywordSet().asString("POLTYPE");
942}
943
944void Scantable::unsetSelection()
945{
946  table_ = originalTable_;
947  attach();
948  selector_.reset();
949}
950
951void Scantable::setSelection( const STSelector& selection )
952{
953  Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
954  if ( tab.nrow() == 0 ) {
955    throw(AipsError("Selection contains no data. Not applying it."));
956  }
957  table_ = tab;
958  attach();
959//   tab.rwKeywordSet().define("nBeam",(Int)(getBeamNos().size())) ;
960//   vector<uint> selectedIFs = getIFNos() ;
961//   Int newnIF = selectedIFs.size() ;
962//   tab.rwKeywordSet().define("nIF",newnIF) ;
963//   if ( newnIF != 0 ) {
964//     Int newnChan = 0 ;
965//     for ( Int i = 0 ; i < newnIF ; i++ ) {
966//       Int nChan = nchan( selectedIFs[i] ) ;
967//       if ( newnChan > nChan )
968//         newnChan = nChan ;
969//     }
970//     tab.rwKeywordSet().define("nChan",newnChan) ;
971//   }
972//   tab.rwKeywordSet().define("nPol",(Int)(getPolNos().size())) ;
973  selector_ = selection;
974}
975
976
977std::string Scantable::headerSummary()
978{
979  // Format header info
980//   STHeader sdh;
981//   sdh = getHeader();
982//   sdh.print();
983  ostringstream oss;
984  oss.flags(std::ios_base::left);
985  String tmp;
986  // Project
987  table_.keywordSet().get("Project", tmp);
988  oss << setw(15) << "Project:" << tmp << endl;
989  // Observation date
990  oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
991  // Observer
992  oss << setw(15) << "Observer:"
993      << table_.keywordSet().asString("Observer") << endl;
994  // Antenna Name
995  table_.keywordSet().get("AntennaName", tmp);
996  oss << setw(15) << "Antenna Name:" << tmp << endl;
997  // Obs type
998  table_.keywordSet().get("Obstype", tmp);
999  // Records (nrow)
1000  oss << setw(15) << "Data Records:" << table_.nrow() << " rows" << endl;
1001  oss << setw(15) << "Obs. Type:" << tmp << endl;
1002  // Beams, IFs, Polarizations, and Channels
1003  oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1004      << setw(15) << "IFs:" << setw(4) << nif() << endl
1005      << setw(15) << "Polarisations:" << setw(4) << npol()
1006      << "(" << getPolType() << ")" << endl
1007      << setw(15) << "Channels:" << nchan() << endl;
1008  // Flux unit
1009  table_.keywordSet().get("FluxUnit", tmp);
1010  oss << setw(15) << "Flux Unit:" << tmp << endl;
1011  // Abscissa Unit
1012  oss << setw(15) << "Abscissa:" << getAbcissaLabel(0) << endl;
1013  // Selection
1014  oss << selector_.print() << endl;
1015
1016  return String(oss);
1017}
1018
1019void Scantable::summary( const std::string& filename )
1020{
1021  ostringstream oss;
1022  ofstream ofs;
1023  LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1024
1025  if (filename != "")
1026    ofs.open( filename.c_str(),  ios::out );
1027
1028  oss << endl;
1029  oss << asap::SEPERATOR << endl;
1030  oss << " Scan Table Summary" << endl;
1031  oss << asap::SEPERATOR << endl;
1032
1033  // Format header info
1034  oss << headerSummary();
1035  oss << endl;
1036
1037  if (table_.nrow() <= 0){
1038    oss << asap::SEPERATOR << endl;
1039    oss << "The MAIN table is empty: there are no data!!!" << endl;
1040    oss << asap::SEPERATOR << endl;
1041
1042    ols << String(oss) << LogIO::POST;
1043    if (ofs) {
1044      ofs << String(oss) << flush;
1045      ofs.close();
1046    }
1047    return;
1048  }
1049
1050
1051
1052  // main table
1053  String dirtype = "Position ("
1054                  + getDirectionRefString()
1055                  + ")";
1056  oss.flags(std::ios_base::left);
1057  oss << setw(5) << "Scan"
1058      << setw(15) << "Source"
1059      << setw(35) << "Time range"
1060      << setw(2) << "" << setw(7) << "Int[s]"
1061      << setw(7) << "Record"
1062      << setw(8) << "SrcType"
1063      << setw(8) << "FreqIDs"
1064      << setw(7) << "MolIDs" << endl;
1065  oss << setw(7)<< "" << setw(6) << "Beam"
1066      << setw(23) << dirtype << endl;
1067
1068  oss << asap::SEPERATOR << endl;
1069
1070  // Flush summary and clear up the string
1071  ols << String(oss) << LogIO::POST;
1072  if (ofs) ofs << String(oss) << flush;
1073  oss.str("");
1074  oss.clear();
1075
1076
1077  // Get Freq_ID map
1078  ROScalarColumn<uInt> ftabIds(frequencies().table(), "ID");
1079  Int nfid = ftabIds.nrow();
1080  if (nfid <= 0){
1081    oss << "FREQUENCIES subtable is empty: there are no data!!!" << endl;
1082    oss << asap::SEPERATOR << endl;
1083
1084    ols << String(oss) << LogIO::POST;
1085    if (ofs) {
1086      ofs << String(oss) << flush;
1087      ofs.close();
1088    }
1089    return;
1090  }
1091  // Storages of overall IFNO, POLNO, and nchan per FREQ_ID
1092  // the orders are identical to ID in FREQ subtable
1093  Block< Vector<uInt> > ifNos(nfid), polNos(nfid);
1094  Vector<Int> fIdchans(nfid,-1);
1095  map<uInt, Int> fidMap;  // (FREQ_ID, row # in FREQ subtable) pair
1096  for (Int i=0; i < nfid; i++){
1097   // fidMap[freqId] returns row number in FREQ subtable
1098   fidMap.insert(pair<uInt, Int>(ftabIds(i),i));
1099   ifNos[i] = Vector<uInt>();
1100   polNos[i] = Vector<uInt>();
1101  }
1102
1103  TableIterator iter(table_, "SCANNO");
1104
1105  // Vars for keeping track of time, freqids, molIds in a SCANNO
1106  Vector<uInt> freqids;
1107  Vector<uInt> molids;
1108  Vector<uInt> beamids(1,0);
1109  Vector<MDirection> beamDirs;
1110  Vector<Int> stypeids(1,0);
1111  Vector<String> stypestrs;
1112  Int nfreq(1);
1113  Int nmol(1);
1114  uInt nbeam(1);
1115  uInt nstype(1);
1116
1117  Double btime(0.0), etime(0.0);
1118  Double meanIntTim(0.0);
1119
1120  uInt currFreqId(0), ftabRow(0);
1121  Int iflen(0), pollen(0);
1122
1123  while (!iter.pastEnd()) {
1124    Table subt = iter.table();
1125    uInt snrow = subt.nrow();
1126    ROTableRow row(subt);
1127    const TableRecord& rec = row.get(0);
1128
1129    // relevant columns
1130    ROScalarColumn<Double> mjdCol(subt,"TIME");
1131    ROScalarColumn<Double> intervalCol(subt,"INTERVAL");
1132    MDirection::ROScalarColumn dirCol(subt,"DIRECTION");
1133
1134    ScalarColumn<uInt> freqIdCol(subt,"FREQ_ID");
1135    ScalarColumn<uInt> molIdCol(subt,"MOLECULE_ID");
1136    ROScalarColumn<uInt> beamCol(subt,"BEAMNO");
1137    ROScalarColumn<Int> stypeCol(subt,"SRCTYPE");
1138
1139    ROScalarColumn<uInt> ifNoCol(subt,"IFNO");
1140    ROScalarColumn<uInt> polNoCol(subt,"POLNO");
1141
1142
1143    // Times
1144    meanIntTim = sum(intervalCol.getColumn()) / (double) snrow;
1145    minMax(btime, etime, mjdCol.getColumn());
1146    etime += meanIntTim/C::day;
1147
1148    // MOLECULE_ID and FREQ_ID
1149    molids = getNumbers(molIdCol);
1150    molids.shape(nmol);
1151
1152    freqids = getNumbers(freqIdCol);
1153    freqids.shape(nfreq);
1154
1155    // Add first beamid, and srcNames
1156    beamids.resize(1,False);
1157    beamDirs.resize(1,False);
1158    beamids(0)=beamCol(0);
1159    beamDirs(0)=dirCol(0);
1160    nbeam = 1;
1161
1162    stypeids.resize(1,False);
1163    stypeids(0)=stypeCol(0);
1164    nstype = 1;
1165
1166    // Global listings of nchan/IFNO/POLNO per FREQ_ID
1167    currFreqId=freqIdCol(0);
1168    ftabRow = fidMap[currFreqId];
1169    // Assumes an identical number of channels per FREQ_ID
1170    if (fIdchans(ftabRow) < 0 ) {
1171      RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1172      fIdchans(ftabRow)=(*spec).shape()(0);
1173    }
1174    // Should keep ifNos and polNos form the previous SCANNO
1175    if ( !anyEQ(ifNos[ftabRow],ifNoCol(0)) ) {
1176      ifNos[ftabRow].shape(iflen);
1177      iflen++;
1178      ifNos[ftabRow].resize(iflen,True);
1179      ifNos[ftabRow](iflen-1) = ifNoCol(0);
1180    }
1181    if ( !anyEQ(polNos[ftabRow],polNoCol(0)) ) {
1182      polNos[ftabRow].shape(pollen);
1183      pollen++;
1184      polNos[ftabRow].resize(pollen,True);
1185      polNos[ftabRow](pollen-1) = polNoCol(0);
1186    }
1187
1188    for (uInt i=1; i < snrow; i++){
1189      // Need to list BEAMNO and DIRECTION in the same order
1190      if ( !anyEQ(beamids,beamCol(i)) ) {
1191        nbeam++;
1192        beamids.resize(nbeam,True);
1193        beamids(nbeam-1)=beamCol(i);
1194        beamDirs.resize(nbeam,True);
1195        beamDirs(nbeam-1)=dirCol(i);
1196      }
1197
1198      // SRCTYPE is Int (getNumber takes only uInt)
1199      if ( !anyEQ(stypeids,stypeCol(i)) ) {
1200        nstype++;
1201        stypeids.resize(nstype,True);
1202        stypeids(nstype-1)=stypeCol(i);
1203      }
1204
1205      // Global listings of nchan/IFNO/POLNO per FREQ_ID
1206      currFreqId=freqIdCol(i);
1207      ftabRow = fidMap[currFreqId];
1208      if (fIdchans(ftabRow) < 0 ) {
1209        const TableRecord& rec = row.get(i);
1210        RORecordFieldPtr< Array<Float> > spec(rec, "SPECTRA");
1211        fIdchans(ftabRow) = (*spec).shape()(0);
1212      }
1213      if ( !anyEQ(ifNos[ftabRow],ifNoCol(i)) ) {
1214        ifNos[ftabRow].shape(iflen);
1215        iflen++;
1216        ifNos[ftabRow].resize(iflen,True);
1217        ifNos[ftabRow](iflen-1) = ifNoCol(i);
1218      }
1219      if ( !anyEQ(polNos[ftabRow],polNoCol(i)) ) {
1220        polNos[ftabRow].shape(pollen);
1221        pollen++;
1222        polNos[ftabRow].resize(pollen,True);
1223        polNos[ftabRow](pollen-1) = polNoCol(i);
1224      }
1225    } // end of row iteration
1226
1227    stypestrs.resize(nstype,False);
1228    for (uInt j=0; j < nstype; j++)
1229      stypestrs(j) = SrcType::getName(stypeids(j));
1230
1231    // Format Scan summary
1232    oss << setw(4) << std::right << rec.asuInt("SCANNO")
1233        << std::left << setw(1) << ""
1234        << setw(15) << rec.asString("SRCNAME")
1235        << setw(21) << MVTime(btime).string(MVTime::YMD,7)
1236        << setw(3) << " - " << MVTime(etime).string(MVTime::TIME,7)
1237        << setw(3) << "" << setw(6) << meanIntTim << setw(1) << ""
1238        << std::right << setw(5) << snrow << setw(2) << ""
1239        << std::left << stypestrs << setw(1) << ""
1240        << freqids << setw(1) << ""
1241        << molids  << endl;
1242    // Format Beam summary
1243    for (uInt j=0; j < nbeam; j++) {
1244      oss << setw(7) << "" << setw(6) << beamids(j) << setw(1) << ""
1245          << formatDirection(beamDirs(j)) << endl;
1246    }
1247    // Flush summary every scan and clear up the string
1248    ols << String(oss) << LogIO::POST;
1249    if (ofs) ofs << String(oss) << flush;
1250    oss.str("");
1251    oss.clear();
1252
1253    ++iter;
1254  } // end of scan iteration
1255  oss << asap::SEPERATOR << endl;
1256 
1257  // List FRECUENCIES Table (using STFrequencies.print may be slow)
1258  oss << "FREQUENCIES: " << nfreq << endl;
1259  oss << std::right << setw(5) << "ID" << setw(2) << ""
1260      << std::left  << setw(5) << "IFNO" << setw(2) << ""
1261      << setw(8) << "Frame"
1262      << setw(16) << "RefVal"
1263      << setw(7) << "RefPix"
1264      << setw(15) << "Increment"
1265      << setw(9) << "Channels"
1266      << setw(6) << "POLNOs" << endl;
1267  Int tmplen;
1268  for (Int i=0; i < nfid; i++){
1269    // List row=i of FREQUENCIES subtable
1270    ifNos[i].shape(tmplen);
1271    if (tmplen == 1) {
1272      oss << std::right << setw(5) << ftabIds(i) << setw(2) << ""
1273          << setw(3) << ifNos[i](0) << setw(1) << ""
1274          << std::left << setw(46) << frequencies().print(ftabIds(i))
1275          << setw(2) << ""
1276          << std::right << setw(8) << fIdchans[i] << setw(2) << ""
1277          << std::left << polNos[i] << endl;
1278    } else if (tmplen > 0 ) {
1279      // You shouldn't come here
1280      oss << std::left
1281          << "Multiple IFNOs in FREQ_ID = " << ftabIds(i)
1282          << " !!!" << endl;
1283    }
1284  }
1285  oss << asap::SEPERATOR << endl;
1286
1287  // List MOLECULES Table (currently lists all rows)
1288  oss << "MOLECULES: " << endl;
1289  if (molecules().nrow() <= 0) {
1290    oss << "   MOLECULES subtable is empty: there are no data" << endl;
1291  } else {
1292    ROTableRow row(molecules().table());
1293    oss << std::right << setw(5) << "ID"
1294        << std::left << setw(3) << ""
1295        << setw(18) << "RestFreq"
1296        << setw(15) << "Name" << endl;
1297    for (Int i=0; i < molecules().nrow(); i++){
1298      const TableRecord& rec=row.get(i);
1299      oss << std::right << setw(5) << rec.asuInt("ID")
1300          << std::left << setw(3) << ""
1301          << rec.asArrayDouble("RESTFREQUENCY") << setw(1) << ""
1302          << rec.asArrayString("NAME") << endl;
1303    }
1304  }
1305  oss << asap::SEPERATOR << endl;
1306  ols << String(oss) << LogIO::POST;
1307  if (ofs) {
1308    ofs << String(oss) << flush;
1309    ofs.close();
1310  }
1311  //  return String(oss);
1312}
1313
1314
1315std::string Scantable::oldheaderSummary()
1316{
1317  // Format header info
1318//   STHeader sdh;
1319//   sdh = getHeader();
1320//   sdh.print();
1321  ostringstream oss;
1322  oss.flags(std::ios_base::left);
1323  oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
1324      << setw(15) << "IFs:" << setw(4) << nif() << endl
1325      << setw(15) << "Polarisations:" << setw(4) << npol()
1326      << "(" << getPolType() << ")" << endl
1327      << setw(15) << "Channels:" << nchan() << endl;
1328  String tmp;
1329  oss << setw(15) << "Observer:"
1330      << table_.keywordSet().asString("Observer") << endl;
1331  oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
1332  table_.keywordSet().get("Project", tmp);
1333  oss << setw(15) << "Project:" << tmp << endl;
1334  table_.keywordSet().get("Obstype", tmp);
1335  oss << setw(15) << "Obs. Type:" << tmp << endl;
1336  table_.keywordSet().get("AntennaName", tmp);
1337  oss << setw(15) << "Antenna Name:" << tmp << endl;
1338  table_.keywordSet().get("FluxUnit", tmp);
1339  oss << setw(15) << "Flux Unit:" << tmp << endl;
1340  int nid = moleculeTable_.nrow();
1341  Bool firstline = True;
1342  oss << setw(15) << "Rest Freqs:";
1343  for (int i=0; i<nid; i++) {
1344    Table t = table_(table_.col("MOLECULE_ID") == i, 1);
1345      if (t.nrow() >  0) {
1346          Vector<Double> vec(moleculeTable_.getRestFrequency(i));
1347          if (vec.nelements() > 0) {
1348               if (firstline) {
1349                   oss << setprecision(10) << vec << " [Hz]" << endl;
1350                   firstline=False;
1351               }
1352               else{
1353                   oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
1354               }
1355          } else {
1356              oss << "none" << endl;
1357          }
1358      }
1359  }
1360
1361  oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
1362  oss << selector_.print() << endl;
1363  return String(oss);
1364}
1365
1366  //std::string Scantable::summary( const std::string& filename )
1367void Scantable::oldsummary( const std::string& filename )
1368{
1369  ostringstream oss;
1370  ofstream ofs;
1371  LogIO ols(LogOrigin("Scantable", "summary", WHERE));
1372
1373  if (filename != "")
1374    ofs.open( filename.c_str(),  ios::out );
1375
1376  oss << endl;
1377  oss << asap::SEPERATOR << endl;
1378  oss << " Scan Table Summary" << endl;
1379  oss << asap::SEPERATOR << endl;
1380
1381  // Format header info
1382  oss << oldheaderSummary();
1383  oss << endl;
1384
1385  // main table
1386  String dirtype = "Position ("
1387                  + getDirectionRefString()
1388                  + ")";
1389  oss.flags(std::ios_base::left);
1390  oss << setw(5) << "Scan" << setw(15) << "Source"
1391      << setw(10) << "Time" << setw(18) << "Integration"
1392      << setw(15) << "Source Type" << endl;
1393  oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
1394  oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
1395      << setw(8) << "Frame" << setw(16)
1396      << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
1397      << setw(7) << "Channels"
1398      << endl;
1399  oss << asap::SEPERATOR << endl;
1400
1401  // Flush summary and clear up the string
1402  ols << String(oss) << LogIO::POST;
1403  if (ofs) ofs << String(oss) << flush;
1404  oss.str("");
1405  oss.clear();
1406
1407  TableIterator iter(table_, "SCANNO");
1408  while (!iter.pastEnd()) {
1409    Table subt = iter.table();
1410    ROTableRow row(subt);
1411    MEpoch::ROScalarColumn timeCol(subt,"TIME");
1412    const TableRecord& rec = row.get(0);
1413    oss << setw(4) << std::right << rec.asuInt("SCANNO")
1414        << std::left << setw(1) << ""
1415        << setw(15) << rec.asString("SRCNAME")
1416        << setw(10) << formatTime(timeCol(0), false);
1417    // count the cycles in the scan
1418    TableIterator cyciter(subt, "CYCLENO");
1419    int nint = 0;
1420    while (!cyciter.pastEnd()) {
1421      ++nint;
1422      ++cyciter;
1423    }
1424    oss << setw(3) << std::right << nint  << setw(3) << " x " << std::left
1425        << setw(11) <<  formatSec(rec.asFloat("INTERVAL")) << setw(1) << ""
1426        << setw(15) << SrcType::getName(rec.asInt("SRCTYPE")) << endl;
1427
1428    TableIterator biter(subt, "BEAMNO");
1429    while (!biter.pastEnd()) {
1430      Table bsubt = biter.table();
1431      ROTableRow brow(bsubt);
1432      const TableRecord& brec = brow.get(0);
1433      uInt row0 = bsubt.rowNumbers(table_)[0];
1434      oss << setw(5) << "" <<  setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
1435      oss  << setw(4) << ""  << formatDirection(getDirection(row0)) << endl;
1436      TableIterator iiter(bsubt, "IFNO");
1437      while (!iiter.pastEnd()) {
1438        Table isubt = iiter.table();
1439        ROTableRow irow(isubt);
1440        const TableRecord& irec = irow.get(0);
1441        oss << setw(9) << "";
1442        oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
1443            << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1444            << setw(3) << "" << nchan(irec.asuInt("IFNO"))
1445            << endl;
1446
1447        ++iiter;
1448      }
1449      ++biter;
1450    }
1451    // Flush summary every scan and clear up the string
1452    ols << String(oss) << LogIO::POST;
1453    if (ofs) ofs << String(oss) << flush;
1454    oss.str("");
1455    oss.clear();
1456
1457    ++iter;
1458  }
1459  oss << asap::SEPERATOR << endl;
1460  ols << String(oss) << LogIO::POST;
1461  if (ofs) {
1462    ofs << String(oss) << flush;
1463    ofs.close();
1464  }
1465  //  return String(oss);
1466}
1467
1468// std::string Scantable::getTime(int whichrow, bool showdate) const
1469// {
1470//   MEpoch::ROScalarColumn timeCol(table_, "TIME");
1471//   MEpoch me;
1472//   if (whichrow > -1) {
1473//     me = timeCol(uInt(whichrow));
1474//   } else {
1475//     Double tm;
1476//     table_.keywordSet().get("UTC",tm);
1477//     me = MEpoch(MVEpoch(tm));
1478//   }
1479//   return formatTime(me, showdate);
1480// }
1481
1482std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
1483{
1484  MEpoch me;
1485  me = getEpoch(whichrow);
1486  return formatTime(me, showdate, prec);
1487}
1488
1489MEpoch Scantable::getEpoch(int whichrow) const
1490{
1491  if (whichrow > -1) {
1492    return timeCol_(uInt(whichrow));
1493  } else {
1494    Double tm;
1495    table_.keywordSet().get("UTC",tm);
1496    return MEpoch(MVEpoch(tm));
1497  }
1498}
1499
1500std::string Scantable::getDirectionString(int whichrow) const
1501{
1502  return formatDirection(getDirection(uInt(whichrow)));
1503}
1504
1505
1506SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1507  const MPosition& mp = getAntennaPosition();
1508  const MDirection& md = getDirection(whichrow);
1509  const MEpoch& me = timeCol_(whichrow);
1510  //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1511  Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1512  return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1513                                          mfreqidCol_(whichrow));
1514}
1515
1516std::vector< double > Scantable::getAbcissa( int whichrow ) const
1517{
1518  if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
1519  std::vector<double> stlout;
1520  int nchan = specCol_(whichrow).nelements();
1521  String us = freqTable_.getUnitString();
1522  if ( us == "" || us == "pixel" || us == "channel" ) {
1523    for (int i=0; i<nchan; ++i) {
1524      stlout.push_back(double(i));
1525    }
1526    return stlout;
1527  }
1528  SpectralCoordinate spc = getSpectralCoordinate(whichrow);
1529  Vector<Double> pixel(nchan);
1530  Vector<Double> world;
1531  indgen(pixel);
1532  if ( Unit(us) == Unit("Hz") ) {
1533    for ( int i=0; i < nchan; ++i) {
1534      Double world;
1535      spc.toWorld(world, pixel[i]);
1536      stlout.push_back(double(world));
1537    }
1538  } else if ( Unit(us) == Unit("km/s") ) {
1539    Vector<Double> world;
1540    spc.pixelToVelocity(world, pixel);
1541    world.tovector(stlout);
1542  }
1543  return stlout;
1544}
1545void Scantable::setDirectionRefString( const std::string & refstr )
1546{
1547  MDirection::Types mdt;
1548  if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1549    throw(AipsError("Illegal Direction frame."));
1550  }
1551  if ( refstr == "" ) {
1552    String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1553    table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1554  } else {
1555    table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1556  }
1557}
1558
1559std::string Scantable::getDirectionRefString( ) const
1560{
1561  return table_.keywordSet().asString("DIRECTIONREF");
1562}
1563
1564MDirection Scantable::getDirection(int whichrow ) const
1565{
1566  String usertype = table_.keywordSet().asString("DIRECTIONREF");
1567  String type = MDirection::showType(dirCol_.getMeasRef().getType());
1568  if ( usertype != type ) {
1569    MDirection::Types mdt;
1570    if (!MDirection::getType(mdt, usertype)) {
1571      throw(AipsError("Illegal Direction frame."));
1572    }
1573    return dirCol_.convert(uInt(whichrow), mdt);
1574  } else {
1575    return dirCol_(uInt(whichrow));
1576  }
1577}
1578
1579std::string Scantable::getAbcissaLabel( int whichrow ) const
1580{
1581  if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
1582  const MPosition& mp = getAntennaPosition();
1583  const MDirection& md = getDirection(whichrow);
1584  const MEpoch& me = timeCol_(whichrow);
1585  //const Double& rf = mmolidCol_(whichrow);
1586  const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1587  SpectralCoordinate spc =
1588    freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1589
1590  String s = "Channel";
1591  Unit u = Unit(freqTable_.getUnitString());
1592  if (u == Unit("km/s")) {
1593    s = CoordinateUtil::axisLabel(spc, 0, True,True,  True);
1594  } else if (u == Unit("Hz")) {
1595    Vector<String> wau(1);wau = u.getName();
1596    spc.setWorldAxisUnits(wau);
1597    s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
1598  }
1599  return s;
1600
1601}
1602
1603/**
1604void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1605                                          const std::string& unit )
1606**/
1607void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1608                                          const std::string& unit )
1609
1610{
1611  ///@todo lookup in line table to fill in name and formattedname
1612  Unit u(unit);
1613  //Quantum<Double> urf(rf, u);
1614  Quantum<Vector<Double> >urf(rf, u);
1615  Vector<String> formattedname(0);
1616  //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1617
1618  //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1619  uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1620  TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1621  tabvec = id;
1622}
1623
1624/**
1625void asap::Scantable::setRestFrequencies( const std::string& name )
1626{
1627  throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1628  ///@todo implement
1629}
1630**/
1631
1632void Scantable::setRestFrequencies( const vector<std::string>& name )
1633{
1634  (void) name; // suppress unused warning
1635  throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1636  ///@todo implement
1637}
1638
1639std::vector< unsigned int > Scantable::rownumbers( ) const
1640{
1641  std::vector<unsigned int> stlout;
1642  Vector<uInt> vec = table_.rowNumbers();
1643  vec.tovector(stlout);
1644  return stlout;
1645}
1646
1647
1648Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1649{
1650  ROTableRow row(table_);
1651  const TableRecord& rec = row.get(whichrow);
1652  Table t =
1653    originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1654                    && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1655                    && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1656                    && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1657  ROArrayColumn<Float> speccol(t, "SPECTRA");
1658  return speccol.getColumn();
1659}
1660
1661std::vector< std::string > Scantable::columnNames( ) const
1662{
1663  Vector<String> vec = table_.tableDesc().columnNames();
1664  return mathutil::tovectorstring(vec);
1665}
1666
1667MEpoch::Types Scantable::getTimeReference( ) const
1668{
1669  return MEpoch::castType(timeCol_.getMeasRef().getType());
1670}
1671
1672void Scantable::addFit( const STFitEntry& fit, int row )
1673{
1674  //cout << mfitidCol_(uInt(row)) << endl;
1675  LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1676  os << mfitidCol_(uInt(row)) << LogIO::POST ;
1677  uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1678  mfitidCol_.put(uInt(row), id);
1679}
1680
1681void Scantable::shift(int npix)
1682{
1683  Vector<uInt> fids(mfreqidCol_.getColumn());
1684  genSort( fids, Sort::Ascending,
1685           Sort::QuickSort|Sort::NoDuplicates );
1686  for (uInt i=0; i<fids.nelements(); ++i) {
1687    frequencies().shiftRefPix(npix, fids[i]);
1688  }
1689}
1690
1691String Scantable::getAntennaName() const
1692{
1693  String out;
1694  table_.keywordSet().get("AntennaName", out);
1695  String::size_type pos1 = out.find("@") ;
1696  String::size_type pos2 = out.find("//") ;
1697  if ( pos2 != String::npos )
1698    out = out.substr(pos2+2,pos1-pos2-2) ;
1699  else if ( pos1 != String::npos )
1700    out = out.substr(0,pos1) ;
1701  return out;
1702}
1703
1704int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1705{
1706  String tbpath;
1707  int ret = 0;
1708  if ( table_.keywordSet().isDefined("GBT_GO") ) {
1709    table_.keywordSet().get("GBT_GO", tbpath);
1710    Table t(tbpath,Table::Old);
1711    // check each scan if other scan of the pair exist
1712    int nscan = scanlist.size();
1713    for (int i = 0; i < nscan; i++) {
1714      Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1715      if (subt.nrow()==0) {
1716        //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1717        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1718        os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1719        ret = 1;
1720        break;
1721      }
1722      ROTableRow row(subt);
1723      const TableRecord& rec = row.get(0);
1724      int scan1seqn = rec.asuInt("PROCSEQN");
1725      int laston1 = rec.asuInt("LASTON");
1726      if ( rec.asuInt("PROCSIZE")==2 ) {
1727        if ( i < nscan-1 ) {
1728          Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1729          if ( subt2.nrow() == 0) {
1730            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1731
1732            //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1733            os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1734            ret = 1;
1735            break;
1736          }
1737          ROTableRow row2(subt2);
1738          const TableRecord& rec2 = row2.get(0);
1739          int scan2seqn = rec2.asuInt("PROCSEQN");
1740          int laston2 = rec2.asuInt("LASTON");
1741          if (scan1seqn == 1 && scan2seqn == 2) {
1742            if (laston1 == laston2) {
1743              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1744              //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1745              os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1746              i +=1;
1747            }
1748            else {
1749              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1750              //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1751              os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1752            }
1753          }
1754          else if (scan1seqn==2 && scan2seqn == 1) {
1755            if (laston1 == laston2) {
1756              LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1757              //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1758              os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1759              ret = 1;
1760              break;
1761            }
1762          }
1763          else {
1764            LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1765            //cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1766            os<<LogIO::WARN<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1767            ret = 1;
1768            break;
1769          }
1770        }
1771      }
1772      else {
1773        LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1774        //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1775        os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1776      }
1777    //if ( i >= nscan ) break;
1778    }
1779  }
1780  else {
1781    LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1782    //cerr<<"No reference to GBT_GO table."<<endl;
1783    os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1784    ret = 1;
1785  }
1786  return ret;
1787}
1788
1789std::vector<double> Scantable::getDirectionVector(int whichrow) const
1790{
1791  Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1792  std::vector<double> dir;
1793  Dir.tovector(dir);
1794  return dir;
1795}
1796
1797void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1798  throw( casa::AipsError )
1799{
1800  // assumed that all rows have same nChan
1801  Vector<Float> arr = specCol_( 0 ) ;
1802  int nChan = arr.nelements() ;
1803
1804  // if nmin < 0 or nmax < 0, nothing to do
1805  if (  nmin < 0 ) {
1806    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1807    }
1808  if (  nmax < 0  ) {
1809    throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1810  }
1811
1812  // if nmin > nmax, exchange values
1813  if ( nmin > nmax ) {
1814    int tmp = nmax ;
1815    nmax = nmin ;
1816    nmin = tmp ;
1817    LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1818    os << "Swap values. Applied range is ["
1819       << nmin << ", " << nmax << "]" << LogIO::POST ;
1820  }
1821
1822  // if nmin exceeds nChan, nothing to do
1823  if ( nmin >= nChan ) {
1824    throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1825  }
1826
1827  // if nmax exceeds nChan, reset nmax to nChan
1828  if ( nmax >= nChan ) {
1829    if ( nmin == 0 ) {
1830      // nothing to do
1831      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1832      os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1833      return ;
1834    }
1835    else {
1836      LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1837      os << "Specified maximum exceeds nChan. Applied range is ["
1838         << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1839      nmax = nChan - 1 ;
1840    }
1841  }
1842
1843  // reshape specCol_ and flagCol_
1844  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1845    reshapeSpectrum( nmin, nmax, irow ) ;
1846  }
1847
1848  // update FREQUENCIES subtable
1849  Double refpix ;
1850  Double refval ;
1851  Double increment ;
1852  int freqnrow = freqTable_.table().nrow() ;
1853  Vector<uInt> oldId( freqnrow ) ;
1854  Vector<uInt> newId( freqnrow ) ;
1855  for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1856    freqTable_.getEntry( refpix, refval, increment, irow ) ;
1857    /***
1858     * need to shift refpix to nmin
1859     * note that channel nmin in old index will be channel 0 in new one
1860     ***/
1861    refval = refval - ( refpix - nmin ) * increment ;
1862    refpix = 0 ;
1863    freqTable_.setEntry( refpix, refval, increment, irow ) ;
1864  }
1865
1866  // update nchan
1867  int newsize = nmax - nmin + 1 ;
1868  table_.rwKeywordSet().define( "nChan", newsize ) ;
1869
1870  // update bandwidth
1871  // assumed all spectra in the scantable have same bandwidth
1872  table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1873
1874  return ;
1875}
1876
1877void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1878{
1879  // reshape specCol_ and flagCol_
1880  Vector<Float> oldspec = specCol_( irow ) ;
1881  Vector<uChar> oldflag = flagsCol_( irow ) ;
1882  uInt newsize = nmax - nmin + 1 ;
1883  specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1884  flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1885
1886  return ;
1887}
1888
1889void asap::Scantable::regridSpecChannel( double dnu, int nChan )
1890{
1891  LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1892  os << "Regrid abcissa with spectral resoultion " << dnu << " " << freqTable_.getUnitString() << " with channel number " << ((nChan>0)? String(nChan) : "covering band width")<< LogIO::POST ;
1893  int freqnrow = freqTable_.table().nrow() ;
1894  Vector<bool> firstTime( freqnrow, true ) ;
1895  double oldincr, factor;
1896  uInt currId;
1897  Double refpix ;
1898  Double refval ;
1899  Double increment ;
1900  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1901    currId = mfreqidCol_(irow);
1902    vector<double> abcissa = getAbcissa( irow ) ;
1903    if (nChan < 0) {
1904      int oldsize = abcissa.size() ;
1905      double bw = (abcissa[oldsize-1]-abcissa[0]) +                     \
1906        0.5 * (abcissa[1]-abcissa[0] + abcissa[oldsize-1]-abcissa[oldsize-2]) ;
1907      nChan = int( ceil( abs(bw/dnu) ) ) ;
1908    }
1909    // actual regridding
1910    regridChannel( nChan, dnu, irow ) ;
1911
1912    // update FREQUENCIES subtable
1913    if (firstTime[currId]) {
1914      oldincr = abcissa[1]-abcissa[0] ;
1915      factor = dnu/oldincr ;
1916      firstTime[currId] = false ;
1917      freqTable_.getEntry( refpix, refval, increment, currId ) ;
1918
1919      //refval = refval - ( refpix + 0.5 * (1 - factor) ) * increment ;
1920      if (factor > 0 ) {
1921        refpix = (refpix + 0.5)/factor - 0.5;
1922      } else {
1923        refpix = (abcissa.size() - 0.5 - refpix)/abs(factor) - 0.5;
1924      }
1925      freqTable_.setEntry( refpix, refval, increment*factor, currId ) ;
1926      //os << "ID" << currId << ": channel width (Orig) = " << oldincr << " [" << freqTable_.getUnitString() << "], scale factor = " << factor << LogIO::POST ;
1927      //os << "     frequency increment (Orig) = " << increment << "-> (New) " << increment*factor << LogIO::POST ;
1928    }
1929  }
1930}
1931
1932void asap::Scantable::regridChannel( int nChan, double dnu )
1933{
1934  LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1935  os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1936  // assumed that all rows have same nChan
1937  Vector<Float> arr = specCol_( 0 ) ;
1938  int oldsize = arr.nelements() ;
1939
1940  // if oldsize == nChan, nothing to do
1941  if ( oldsize == nChan ) {
1942    os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1943    return ;
1944  }
1945
1946  // if oldChan < nChan, unphysical operation
1947  if ( oldsize < nChan ) {
1948    os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1949    return ;
1950  }
1951
1952  // change channel number for specCol_, flagCol_, and tsysCol_ (if necessary)
1953  vector<string> coordinfo = getCoordInfo() ;
1954  string oldinfo = coordinfo[0] ;
1955  coordinfo[0] = "Hz" ;
1956  setCoordInfo( coordinfo ) ;
1957  for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1958    regridChannel( nChan, dnu, irow ) ;
1959  }
1960  coordinfo[0] = oldinfo ;
1961  setCoordInfo( coordinfo ) ;
1962
1963
1964  // NOTE: this method does not update metadata such as
1965  //       FREQUENCIES subtable, nChan, Bandwidth, etc.
1966
1967  return ;
1968}
1969
1970void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1971{
1972  // logging
1973  //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1974  //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1975
1976  Vector<Float> oldspec = specCol_( irow ) ;
1977  Vector<uChar> oldflag = flagsCol_( irow ) ;
1978  Vector<Float> oldtsys = tsysCol_( irow ) ;
1979  Vector<Float> newspec( nChan, 0 ) ;
1980  Vector<uChar> newflag( nChan, true ) ;
1981  Vector<Float> newtsys ;
1982  bool regridTsys = false ;
1983  if (oldtsys.size() == oldspec.size()) {
1984    regridTsys = true ;
1985    newtsys.resize(nChan,false) ;
1986    newtsys = 0 ;
1987  }
1988
1989  // regrid
1990  vector<double> abcissa = getAbcissa( irow ) ;
1991  int oldsize = abcissa.size() ;
1992  double olddnu = abcissa[1] - abcissa[0] ;
1993  //int ichan = 0 ;
1994  double wsum = 0.0 ;
1995  Vector<double> zi( nChan+1 ) ;
1996  Vector<double> yi( oldsize + 1 ) ;
1997  yi[0] = abcissa[0] - 0.5 * olddnu ;
1998  for ( int ii = 1 ; ii < oldsize ; ii++ )
1999    yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
2000  yi[oldsize] = abcissa[oldsize-1] \
2001    + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
2002  //zi[0] = abcissa[0] - 0.5 * olddnu ;
2003  zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
2004  for ( int ii = 1 ; ii < nChan ; ii++ )
2005    zi[ii] = zi[0] + dnu * ii ;
2006  zi[nChan] = zi[nChan-1] + dnu ;
2007  // Access zi and yi in ascending order
2008  int izs = ((dnu > 0) ? 0 : nChan ) ;
2009  int ize = ((dnu > 0) ? nChan : 0 ) ;
2010  int izincr = ((dnu > 0) ? 1 : -1 ) ;
2011  int ichan =  ((olddnu > 0) ? 0 : oldsize ) ;
2012  int iye = ((olddnu > 0) ? oldsize : 0 ) ;
2013  int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
2014  //for ( int ii = izs ; ii != ize ; ii+=izincr ){
2015  int ii = izs ;
2016  while (ii != ize) {
2017    // always zl < zr
2018    double zl = zi[ii] ;
2019    double zr = zi[ii+izincr] ;
2020    // Need to access smaller index for the new spec, flag, and tsys.
2021    // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
2022    int i = min(ii, ii+izincr) ;
2023    //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
2024    int jj = ichan ;
2025    while (jj != iye) {
2026      // always yl < yr
2027      double yl = yi[jj] ;
2028      double yr = yi[jj+iyincr] ;
2029      // Need to access smaller index for the original spec, flag, and tsys.
2030      // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
2031      int j = min(jj, jj+iyincr) ;
2032      if ( yr <= zl ) {
2033        jj += iyincr ;
2034        continue ;
2035      }
2036      else if ( yl <= zl ) {
2037        if ( yr < zr ) {
2038          if (!oldflag[j]) {
2039            newspec[i] += oldspec[j] * ( yr - zl ) ;
2040            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
2041            wsum += ( yr - zl ) ;
2042          }
2043          newflag[i] = newflag[i] && oldflag[j] ;
2044        }
2045        else {
2046          if (!oldflag[j]) {
2047            newspec[i] += oldspec[j] * abs(dnu) ;
2048            if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
2049            wsum += abs(dnu) ;
2050          }
2051          newflag[i] = newflag[i] && oldflag[j] ;
2052          ichan = jj ;
2053          break ;
2054        }
2055      }
2056      else if ( yl < zr ) {
2057        if ( yr <= zr ) {
2058          if (!oldflag[j]) {
2059            newspec[i] += oldspec[j] * ( yr - yl ) ;
2060            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
2061            wsum += ( yr - yl ) ;
2062          }
2063          newflag[i] = newflag[i] && oldflag[j] ;
2064        }
2065        else {
2066          if (!oldflag[j]) {
2067            newspec[i] += oldspec[j] * ( zr - yl ) ;
2068            if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
2069            wsum += ( zr - yl ) ;
2070          }
2071          newflag[i] = newflag[i] && oldflag[j] ;
2072          ichan = jj ;
2073          break ;
2074        }
2075      }
2076      else {
2077        ichan = jj - iyincr ;
2078        break ;
2079      }
2080      jj += iyincr ;
2081    }
2082    if ( wsum != 0.0 ) {
2083      newspec[i] /= wsum ;
2084      if (regridTsys) newtsys[i] /= wsum ;
2085    }
2086    wsum = 0.0 ;
2087    ii += izincr ;
2088  }
2089//   if ( dnu > 0.0 ) {
2090//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
2091//       double zl = zi[ii] ;
2092//       double zr = zi[ii+1] ;
2093//       for ( int j = ichan ; j < oldsize ; j++ ) {
2094//         double yl = yi[j] ;
2095//         double yr = yi[j+1] ;
2096//         if ( yl <= zl ) {
2097//           if ( yr <= zl ) {
2098//             continue ;
2099//           }
2100//           else if ( yr <= zr ) {
2101//          if (!oldflag[j]) {
2102//            newspec[ii] += oldspec[j] * ( yr - zl ) ;
2103//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
2104//            wsum += ( yr - zl ) ;
2105//          }
2106//          newflag[ii] = newflag[ii] && oldflag[j] ;
2107//           }
2108//           else {
2109//          if (!oldflag[j]) {
2110//            newspec[ii] += oldspec[j] * dnu ;
2111//            if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
2112//            wsum += dnu ;
2113//          }
2114//          newflag[ii] = newflag[ii] && oldflag[j] ;
2115//             ichan = j ;
2116//             break ;
2117//           }
2118//         }
2119//         else if ( yl < zr ) {
2120//           if ( yr <= zr ) {
2121//          if (!oldflag[j]) {
2122//            newspec[ii] += oldspec[j] * ( yr - yl ) ;
2123//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
2124//               wsum += ( yr - yl ) ;
2125//          }
2126//          newflag[ii] = newflag[ii] && oldflag[j] ;
2127//           }
2128//           else {
2129//          if (!oldflag[j]) {
2130//            newspec[ii] += oldspec[j] * ( zr - yl ) ;
2131//            if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
2132//            wsum += ( zr - yl ) ;
2133//          }
2134//          newflag[ii] = newflag[ii] && oldflag[j] ;
2135//             ichan = j ;
2136//             break ;
2137//           }
2138//         }
2139//         else {
2140//           ichan = j - 1 ;
2141//           break ;
2142//         }
2143//       }
2144//       if ( wsum != 0.0 ) {
2145//         newspec[ii] /= wsum ;
2146//      if (regridTsys) newtsys[ii] /= wsum ;
2147//       }
2148//       wsum = 0.0 ;
2149//     }
2150//   }
2151//   else if ( dnu < 0.0 ) {
2152//     for ( int ii = 0 ; ii < nChan ; ii++ ) {
2153//       double zl = zi[ii] ;
2154//       double zr = zi[ii+1] ;
2155//       for ( int j = ichan ; j < oldsize ; j++ ) {
2156//         double yl = yi[j] ;
2157//         double yr = yi[j+1] ;
2158//         if ( yl >= zl ) {
2159//           if ( yr >= zl ) {
2160//             continue ;
2161//           }
2162//           else if ( yr >= zr ) {
2163//          if (!oldflag[j]) {
2164//            newspec[ii] += oldspec[j] * abs( yr - zl ) ;
2165//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
2166//            wsum += abs( yr - zl ) ;
2167//          }
2168//          newflag[ii] = newflag[ii] && oldflag[j] ;
2169//           }
2170//           else {
2171//          if (!oldflag[j]) {
2172//            newspec[ii] += oldspec[j] * abs( dnu ) ;
2173//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
2174//            wsum += abs( dnu ) ;
2175//          }
2176//          newflag[ii] = newflag[ii] && oldflag[j] ;
2177//             ichan = j ;
2178//             break ;
2179//           }
2180//         }
2181//         else if ( yl > zr ) {
2182//           if ( yr >= zr ) {
2183//          if (!oldflag[j]) {
2184//            newspec[ii] += oldspec[j] * abs( yr - yl ) ;
2185//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
2186//            wsum += abs( yr - yl ) ;
2187//          }
2188//          newflag[ii] = newflag[ii] && oldflag[j] ;
2189//           }
2190//           else {
2191//          if (!oldflag[j]) {
2192//            newspec[ii] += oldspec[j] * abs( zr - yl ) ;
2193//            if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
2194//            wsum += abs( zr - yl ) ;
2195//          }
2196//          newflag[ii] = newflag[ii] && oldflag[j] ;
2197//             ichan = j ;
2198//             break ;
2199//           }
2200//         }
2201//         else {
2202//           ichan = j - 1 ;
2203//           break ;
2204//         }
2205//       }
2206//       if ( wsum != 0.0 ) {
2207//         newspec[ii] /= wsum ;
2208//      if (regridTsys) newtsys[ii] /= wsum ;
2209//       }
2210//       wsum = 0.0 ;
2211//     }
2212//   }
2213// //   //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
2214// //   pile += dnu ;
2215// //   wedge = olddnu * ( refChan + 1 ) ;
2216// //   while ( wedge < pile ) {
2217// //     newspec[0] += olddnu * oldspec[refChan] ;
2218// //     newflag[0] = newflag[0] || oldflag[refChan] ;
2219// //     //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
2220// //     refChan++ ;
2221// //     wedge += olddnu ;
2222// //     wsum += olddnu ;
2223// //     //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2224// //   }
2225// //   frac = ( wedge - pile ) / olddnu ;
2226// //   wsum += ( 1.0 - frac ) * olddnu ;
2227// //   newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2228// //   newflag[0] = newflag[0] || oldflag[refChan] ;
2229// //   //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
2230// //   //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2231// //   newspec[0] /= wsum ;
2232// //   //ofs << "newspec[0] = " << newspec[0] << endl ;
2233// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2234
2235// //   /***
2236// //    * ichan = 1 - nChan-2
2237// //    ***/
2238// //   for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
2239// //     pile += dnu ;
2240// //     newspec[ichan] += frac * olddnu * oldspec[refChan] ;
2241// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2242// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
2243// //     refChan++ ;
2244// //     wedge += olddnu ;
2245// //     wsum = frac * olddnu ;
2246// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2247// //     while ( wedge < pile ) {
2248// //       newspec[ichan] += olddnu * oldspec[refChan] ;
2249// //       newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2250// //       //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
2251// //       refChan++ ;
2252// //       wedge += olddnu ;
2253// //       wsum += olddnu ;
2254// //       //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2255// //     }
2256// //     frac = ( wedge - pile ) / olddnu ;
2257// //     wsum += ( 1.0 - frac ) * olddnu ;
2258// //     newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2259// //     newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2260// //     //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
2261// //     //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2262// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2263// //     newspec[ichan] /= wsum ;
2264// //     //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
2265// //   }
2266
2267// //   /***
2268// //    * ichan = nChan-1
2269// //    ***/
2270// //   // NOTE: Assumed that all spectra have the same bandwidth
2271// //   pile += dnu ;
2272// //   newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
2273// //   newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
2274// //   //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2275// //   refChan++ ;
2276// //   wedge += olddnu ;
2277// //   wsum = frac * olddnu ;
2278// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2279// //   for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
2280// //     newspec[nChan-1] += olddnu * oldspec[jchan] ;
2281// //     newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
2282// //     wsum += olddnu ;
2283// //     //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2284// //     //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2285// //   }
2286// //   //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2287// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2288// //   newspec[nChan-1] /= wsum ;
2289// //   //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
2290
2291// //   // ofs.close() ;
2292
2293  specCol_.put( irow, newspec ) ;
2294  flagsCol_.put( irow, newflag ) ;
2295  if (regridTsys) tsysCol_.put( irow, newtsys );
2296
2297  return ;
2298}
2299
2300std::vector<float> Scantable::getWeather(int whichrow) const
2301{
2302  std::vector<float> out(5);
2303  //Float temperature, pressure, humidity, windspeed, windaz;
2304  weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
2305                         mweatheridCol_(uInt(whichrow)));
2306
2307
2308  return out;
2309}
2310
2311bool Scantable::getFlagtraFast(uInt whichrow)
2312{
2313  uChar flag;
2314  Vector<uChar> flags;
2315  flagsCol_.get(whichrow, flags);
2316  flag = flags[0];
2317  for (uInt i = 1; i < flags.size(); ++i) {
2318    flag &= flags[i];
2319  }
2320  return ((flag >> 7) == 1);
2321}
2322
2323void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2324{
2325  try {
2326    ofstream ofs;
2327    String coordInfo = "";
2328    bool hasSameNchan = true;
2329    bool outTextFile = false;
2330
2331    if (blfile != "") {
2332      ofs.open(blfile.c_str(), ios::out | ios::app);
2333      if (ofs) outTextFile = true;
2334    }
2335
2336    if (outLogger || outTextFile) {
2337      coordInfo = getCoordInfo()[0];
2338      if (coordInfo == "") coordInfo = "channel";
2339      hasSameNchan = hasSameNchanOverIFs();
2340    }
2341
2342    Fitter fitter = Fitter();
2343    fitter.setExpression("poly", order);
2344    //fitter.setIterClipping(thresClip, nIterClip);
2345
2346    int nRow = nrow();
2347    std::vector<bool> chanMask;
2348    bool showProgress;
2349    int minNRow;
2350    parseProgressInfo(progressInfo, showProgress, minNRow);
2351
2352    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2353      chanMask = getCompositeChanMask(whichrow, mask);
2354      fitBaseline(chanMask, whichrow, fitter);
2355      setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2356      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
2357      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2358    }
2359
2360    if (outTextFile) ofs.close();
2361
2362  } catch (...) {
2363    throw;
2364  }
2365}
2366
2367void Scantable::autoPolyBaseline(const std::vector<bool>& mask, int order, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2368{
2369  try {
2370    ofstream ofs;
2371    String coordInfo = "";
2372    bool hasSameNchan = true;
2373    bool outTextFile = false;
2374
2375    if (blfile != "") {
2376      ofs.open(blfile.c_str(), ios::out | ios::app);
2377      if (ofs) outTextFile = true;
2378    }
2379
2380    if (outLogger || outTextFile) {
2381      coordInfo = getCoordInfo()[0];
2382      if (coordInfo == "") coordInfo = "channel";
2383      hasSameNchan = hasSameNchanOverIFs();
2384    }
2385
2386    Fitter fitter = Fitter();
2387    fitter.setExpression("poly", order);
2388    //fitter.setIterClipping(thresClip, nIterClip);
2389
2390    int nRow = nrow();
2391    std::vector<bool> chanMask;
2392    int minEdgeSize = getIFNos().size()*2;
2393    STLineFinder lineFinder = STLineFinder();
2394    lineFinder.setOptions(threshold, 3, chanAvgLimit);
2395
2396    bool showProgress;
2397    int minNRow;
2398    parseProgressInfo(progressInfo, showProgress, minNRow);
2399
2400    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2401
2402      //-------------------------------------------------------
2403      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2404      //-------------------------------------------------------
2405      int edgeSize = edge.size();
2406      std::vector<int> currentEdge;
2407      if (edgeSize >= 2) {
2408        int idx = 0;
2409        if (edgeSize > 2) {
2410          if (edgeSize < minEdgeSize) {
2411            throw(AipsError("Length of edge element info is less than that of IFs"));
2412          }
2413          idx = 2 * getIF(whichrow);
2414        }
2415        currentEdge.push_back(edge[idx]);
2416        currentEdge.push_back(edge[idx+1]);
2417      } else {
2418        throw(AipsError("Wrong length of edge element"));
2419      }
2420      lineFinder.setData(getSpectrum(whichrow));
2421      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2422      chanMask = lineFinder.getMask();
2423      //-------------------------------------------------------
2424
2425      fitBaseline(chanMask, whichrow, fitter);
2426      setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2427
2428      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
2429      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2430    }
2431
2432    if (outTextFile) ofs.close();
2433
2434  } catch (...) {
2435    throw;
2436  }
2437}
2438
2439void Scantable::cubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2440{
2441  try {
2442    ofstream ofs;
2443    String coordInfo = "";
2444    bool hasSameNchan = true;
2445    bool outTextFile = false;
2446
2447    if (blfile != "") {
2448      ofs.open(blfile.c_str(), ios::out | ios::app);
2449      if (ofs) outTextFile = true;
2450    }
2451
2452    if (outLogger || outTextFile) {
2453      coordInfo = getCoordInfo()[0];
2454      if (coordInfo == "") coordInfo = "channel";
2455      hasSameNchan = hasSameNchanOverIFs();
2456    }
2457
2458    //Fitter fitter = Fitter();
2459    //fitter.setExpression("cspline", nPiece);
2460    //fitter.setIterClipping(thresClip, nIterClip);
2461
2462    bool showProgress;
2463    int minNRow;
2464    parseProgressInfo(progressInfo, showProgress, minNRow);
2465
2466    int nRow = nrow();
2467    std::vector<bool> chanMask;
2468
2469    //--------------------------------
2470    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2471      chanMask = getCompositeChanMask(whichrow, mask);
2472      //fitBaseline(chanMask, whichrow, fitter);
2473      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2474      std::vector<int> pieceEdges(nPiece+1);
2475      std::vector<float> params(nPiece*4);
2476      int nClipped = 0;
2477      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2478      setSpectrum(res, whichrow);
2479      //
2480
2481      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
2482      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2483    }
2484    //--------------------------------
2485   
2486    if (outTextFile) ofs.close();
2487
2488  } catch (...) {
2489    throw;
2490  }
2491}
2492
2493void Scantable::autoCubicSplineBaseline(const std::vector<bool>& mask, int nPiece, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2494{
2495  try {
2496    ofstream ofs;
2497    String coordInfo = "";
2498    bool hasSameNchan = true;
2499    bool outTextFile = false;
2500
2501    if (blfile != "") {
2502      ofs.open(blfile.c_str(), ios::out | ios::app);
2503      if (ofs) outTextFile = true;
2504    }
2505
2506    if (outLogger || outTextFile) {
2507      coordInfo = getCoordInfo()[0];
2508      if (coordInfo == "") coordInfo = "channel";
2509      hasSameNchan = hasSameNchanOverIFs();
2510    }
2511
2512    //Fitter fitter = Fitter();
2513    //fitter.setExpression("cspline", nPiece);
2514    //fitter.setIterClipping(thresClip, nIterClip);
2515
2516    int nRow = nrow();
2517    std::vector<bool> chanMask;
2518    int minEdgeSize = getIFNos().size()*2;
2519    STLineFinder lineFinder = STLineFinder();
2520    lineFinder.setOptions(threshold, 3, chanAvgLimit);
2521
2522    bool showProgress;
2523    int minNRow;
2524    parseProgressInfo(progressInfo, showProgress, minNRow);
2525
2526    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2527
2528      //-------------------------------------------------------
2529      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2530      //-------------------------------------------------------
2531      int edgeSize = edge.size();
2532      std::vector<int> currentEdge;
2533      if (edgeSize >= 2) {
2534        int idx = 0;
2535        if (edgeSize > 2) {
2536          if (edgeSize < minEdgeSize) {
2537            throw(AipsError("Length of edge element info is less than that of IFs"));
2538          }
2539          idx = 2 * getIF(whichrow);
2540        }
2541        currentEdge.push_back(edge[idx]);
2542        currentEdge.push_back(edge[idx+1]);
2543      } else {
2544        throw(AipsError("Wrong length of edge element"));
2545      }
2546      lineFinder.setData(getSpectrum(whichrow));
2547      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2548      chanMask = lineFinder.getMask();
2549      //-------------------------------------------------------
2550
2551
2552      //fitBaseline(chanMask, whichrow, fitter);
2553      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2554      std::vector<int> pieceEdges(nPiece+1);
2555      std::vector<float> params(nPiece*4);
2556      int nClipped = 0;
2557      std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2558      setSpectrum(res, whichrow);
2559      //
2560
2561      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
2562      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2563    }
2564
2565    if (outTextFile) ofs.close();
2566
2567  } catch (...) {
2568    throw;
2569  }
2570}
2571
2572std::vector<float> Scantable::doCubicSplineFitting(const std::vector<float>& data, const std::vector<bool>& mask, int nPiece, std::vector<int>& idxEdge, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
2573{
2574  if (data.size() != mask.size()) {
2575    throw(AipsError("data and mask sizes are not identical"));
2576  }
2577  if (nPiece < 1) {
2578    throw(AipsError("number of the sections must be one or more"));
2579  }
2580
2581  int nChan = data.size();
2582  std::vector<int> maskArray(nChan);
2583  std::vector<int> x(nChan);
2584  int j = 0;
2585  for (int i = 0; i < nChan; ++i) {
2586    maskArray[i] = mask[i] ? 1 : 0;
2587    if (mask[i]) {
2588      x[j] = i;
2589      j++;
2590    }
2591  }
2592  int initNData = j;
2593
2594  if (initNData < nPiece) {
2595    throw(AipsError("too few non-flagged channels"));
2596  }
2597
2598  int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
2599  std::vector<double> invEdge(nPiece-1);
2600  idxEdge[0] = x[0];
2601  for (int i = 1; i < nPiece; ++i) {
2602    int valX = x[nElement*i];
2603    idxEdge[i] = valX;
2604    invEdge[i-1] = 1.0/(double)valX;
2605  }
2606  idxEdge[nPiece] = x[initNData-1]+1;
2607
2608  int nData = initNData;
2609  int nDOF = nPiece + 3;  //number of parameters to solve, namely, 4+(nPiece-1).
2610
2611  std::vector<double> x1(nChan), x2(nChan), x3(nChan);
2612  std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
2613  std::vector<double> r1(nChan), residual(nChan);
2614  for (int i = 0; i < nChan; ++i) {
2615    double di = (double)i;
2616    double dD = (double)data[i];
2617    x1[i]   = di;
2618    x2[i]   = di*di;
2619    x3[i]   = di*di*di;
2620    z1[i]   = dD;
2621    x1z1[i] = dD*di;
2622    x2z1[i] = dD*di*di;
2623    x3z1[i] = dD*di*di*di;
2624    r1[i]   = 0.0;
2625    residual[i] = 0.0;
2626  }
2627
2628  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
2629    // xMatrix : horizontal concatenation of
2630    //           the least-sq. matrix (left) and an
2631    //           identity matrix (right).
2632    // the right part is used to calculate the inverse matrix of the left part.
2633    double xMatrix[nDOF][2*nDOF];
2634    double zMatrix[nDOF];
2635    for (int i = 0; i < nDOF; ++i) {
2636      for (int j = 0; j < 2*nDOF; ++j) {
2637        xMatrix[i][j] = 0.0;
2638      }
2639      xMatrix[i][nDOF+i] = 1.0;
2640      zMatrix[i] = 0.0;
2641    }
2642
2643    for (int n = 0; n < nPiece; ++n) {
2644      int nUseDataInPiece = 0;
2645      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2646
2647        if (maskArray[i] == 0) continue;
2648
2649        xMatrix[0][0] += 1.0;
2650        xMatrix[0][1] += x1[i];
2651        xMatrix[0][2] += x2[i];
2652        xMatrix[0][3] += x3[i];
2653        xMatrix[1][1] += x2[i];
2654        xMatrix[1][2] += x3[i];
2655        xMatrix[1][3] += x2[i]*x2[i];
2656        xMatrix[2][2] += x2[i]*x2[i];
2657        xMatrix[2][3] += x3[i]*x2[i];
2658        xMatrix[3][3] += x3[i]*x3[i];
2659        zMatrix[0] += z1[i];
2660        zMatrix[1] += x1z1[i];
2661        zMatrix[2] += x2z1[i];
2662        zMatrix[3] += x3z1[i];
2663
2664        for (int j = 0; j < n; ++j) {
2665          double q = 1.0 - x1[i]*invEdge[j];
2666          q = q*q*q;
2667          xMatrix[0][j+4] += q;
2668          xMatrix[1][j+4] += q*x1[i];
2669          xMatrix[2][j+4] += q*x2[i];
2670          xMatrix[3][j+4] += q*x3[i];
2671          for (int k = 0; k < j; ++k) {
2672            double r = 1.0 - x1[i]*invEdge[k];
2673            r = r*r*r;
2674            xMatrix[k+4][j+4] += r*q;
2675          }
2676          xMatrix[j+4][j+4] += q*q;
2677          zMatrix[j+4] += q*z1[i];
2678        }
2679
2680        nUseDataInPiece++;
2681      }
2682
2683      if (nUseDataInPiece < 1) {
2684        std::vector<string> suffixOfPieceNumber(4);
2685        suffixOfPieceNumber[0] = "th";
2686        suffixOfPieceNumber[1] = "st";
2687        suffixOfPieceNumber[2] = "nd";
2688        suffixOfPieceNumber[3] = "rd";
2689        int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
2690        ostringstream oss;
2691        oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
2692        oss << " piece of the spectrum. can't execute fitting anymore.";
2693        throw(AipsError(String(oss)));
2694      }
2695    }
2696
2697    for (int i = 0; i < nDOF; ++i) {
2698      for (int j = 0; j < i; ++j) {
2699        xMatrix[i][j] = xMatrix[j][i];
2700      }
2701    }
2702
2703    std::vector<double> invDiag(nDOF);
2704    for (int i = 0; i < nDOF; ++i) {
2705      invDiag[i] = 1.0/xMatrix[i][i];
2706      for (int j = 0; j < nDOF; ++j) {
2707        xMatrix[i][j] *= invDiag[i];
2708      }
2709    }
2710
2711    for (int k = 0; k < nDOF; ++k) {
2712      for (int i = 0; i < nDOF; ++i) {
2713        if (i != k) {
2714          double factor1 = xMatrix[k][k];
2715          double factor2 = xMatrix[i][k];
2716          for (int j = k; j < 2*nDOF; ++j) {
2717            xMatrix[i][j] *= factor1;
2718            xMatrix[i][j] -= xMatrix[k][j]*factor2;
2719            xMatrix[i][j] /= factor1;
2720          }
2721        }
2722      }
2723      double xDiag = xMatrix[k][k];
2724      for (int j = k; j < 2*nDOF; ++j) {
2725        xMatrix[k][j] /= xDiag;
2726      }
2727    }
2728   
2729    for (int i = 0; i < nDOF; ++i) {
2730      for (int j = 0; j < nDOF; ++j) {
2731        xMatrix[i][nDOF+j] *= invDiag[j];
2732      }
2733    }
2734    //compute a vector y which consists of the coefficients of the best-fit spline curves
2735    //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
2736    //cubic terms for the other pieces (in case nPiece>1).
2737    std::vector<double> y(nDOF);
2738    for (int i = 0; i < nDOF; ++i) {
2739      y[i] = 0.0;
2740      for (int j = 0; j < nDOF; ++j) {
2741        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
2742      }
2743    }
2744
2745    double a0 = y[0];
2746    double a1 = y[1];
2747    double a2 = y[2];
2748    double a3 = y[3];
2749
2750    int j = 0;
2751    for (int n = 0; n < nPiece; ++n) {
2752      for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2753        r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
2754      }
2755      params[j]   = a0;
2756      params[j+1] = a1;
2757      params[j+2] = a2;
2758      params[j+3] = a3;
2759      j += 4;
2760
2761      if (n == nPiece-1) break;
2762
2763      double d = y[4+n];
2764      double iE = invEdge[n];
2765      a0 +=     d;
2766      a1 -= 3.0*d*iE;
2767      a2 += 3.0*d*iE*iE;
2768      a3 -=     d*iE*iE*iE;
2769    }
2770
2771    //subtract constant value for masked regions at the edge of spectrum
2772    if (idxEdge[0] > 0) {
2773      int n = idxEdge[0];
2774      for (int i = 0; i < idxEdge[0]; ++i) {
2775        //--cubic extrapolate--
2776        //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
2777        //--linear extrapolate--
2778        //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
2779        //--constant--
2780        r1[i] = r1[n];
2781      }
2782    }
2783    if (idxEdge[nPiece] < nChan) {
2784      int n = idxEdge[nPiece]-1;
2785      for (int i = idxEdge[nPiece]; i < nChan; ++i) {
2786        //--cubic extrapolate--
2787        //int m = 4*(nPiece-1);
2788        //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
2789        //--linear extrapolate--
2790        //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
2791        //--constant--
2792        r1[i] = r1[n];
2793      }
2794    }
2795
2796    for (int i = 0; i < nChan; ++i) {
2797      residual[i] = z1[i] - r1[i];
2798    }
2799
2800    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
2801      break;
2802    } else {
2803      double stdDev = 0.0;
2804      for (int i = 0; i < nChan; ++i) {
2805        stdDev += residual[i]*residual[i]*(double)maskArray[i];
2806      }
2807      stdDev = sqrt(stdDev/(double)nData);
2808     
2809      double thres = stdDev * thresClip;
2810      int newNData = 0;
2811      for (int i = 0; i < nChan; ++i) {
2812        if (abs(residual[i]) >= thres) {
2813          maskArray[i] = 0;
2814        }
2815        if (maskArray[i] > 0) {
2816          newNData++;
2817        }
2818      }
2819      if (newNData == nData) {
2820        break; //no more flag to add. iteration stops.
2821      } else {
2822        nData = newNData;
2823      }
2824    }
2825  }
2826
2827  nClipped = initNData - nData;
2828
2829  std::vector<float> result(nChan);
2830  if (getResidual) {
2831    for (int i = 0; i < nChan; ++i) {
2832      result[i] = (float)residual[i];
2833    }
2834  } else {
2835    for (int i = 0; i < nChan; ++i) {
2836      result[i] = (float)r1[i];
2837    }
2838  }
2839
2840  return result;
2841}
2842
2843void Scantable::selectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
2844{
2845  nWaves.clear();
2846
2847  if (applyFFT) {
2848    string fftThAttr;
2849    float fftThSigma;
2850    int fftThTop;
2851    parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
2852    doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
2853  }
2854
2855  addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
2856}
2857
2858void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
2859{
2860  uInt idxSigma = fftThresh.find("sigma");
2861  uInt idxTop   = fftThresh.find("top");
2862
2863  if (idxSigma == fftThresh.size() - 5) {
2864    std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
2865    is >> fftThSigma;
2866    fftThAttr = "sigma";
2867  } else if (idxTop == 0) {
2868    std::istringstream is(fftThresh.substr(3));
2869    is >> fftThTop;
2870    fftThAttr = "top";
2871  } else {
2872    bool isNumber = true;
2873    for (uInt i = 0; i < fftThresh.size()-1; ++i) {
2874      char ch = (fftThresh.substr(i, 1).c_str())[0];
2875      if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
2876        isNumber = false;
2877        break;
2878      }
2879    }
2880    if (isNumber) {
2881      std::istringstream is(fftThresh);
2882      is >> fftThSigma;
2883      fftThAttr = "sigma";
2884    } else {
2885      throw(AipsError("fftthresh has a wrong value"));
2886    }
2887  }
2888}
2889
2890void Scantable::doSelectWaveNumbers(const int whichrow, const std::vector<bool>& chanMask, const std::string& fftMethod, const float fftThSigma, const int fftThTop, const std::string& fftThAttr, std::vector<int>& nWaves)
2891{
2892  std::vector<float> fspec;
2893  if (fftMethod == "fft") {
2894    fspec = execFFT(whichrow, chanMask, false, true);
2895  //} else if (fftMethod == "lsp") {
2896  //  fspec = lombScarglePeriodogram(whichrow);
2897  }
2898
2899  if (fftThAttr == "sigma") {
2900    float mean  = 0.0;
2901    float mean2 = 0.0;
2902    for (uInt i = 0; i < fspec.size(); ++i) {
2903      mean  += fspec[i];
2904      mean2 += fspec[i]*fspec[i];
2905    }
2906    mean  /= float(fspec.size());
2907    mean2 /= float(fspec.size());
2908    float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
2909
2910    for (uInt i = 0; i < fspec.size(); ++i) {
2911      if (fspec[i] >= thres) {
2912        nWaves.push_back(i);
2913      }
2914    }
2915
2916  } else if (fftThAttr == "top") {
2917    for (int i = 0; i < fftThTop; ++i) {
2918      float max = 0.0;
2919      int maxIdx = 0;
2920      for (uInt j = 0; j < fspec.size(); ++j) {
2921        if (fspec[j] > max) {
2922          max = fspec[j];
2923          maxIdx = j;
2924        }
2925      }
2926      nWaves.push_back(maxIdx);
2927      fspec[maxIdx] = 0.0;
2928    }
2929
2930  }
2931
2932  if (nWaves.size() > 1) {
2933    sort(nWaves.begin(), nWaves.end());
2934  }
2935}
2936
2937void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
2938{
2939  std::vector<int> tempAddNWaves, tempRejectNWaves;
2940  for (uInt i = 0; i < addNWaves.size(); ++i) {
2941    tempAddNWaves.push_back(addNWaves[i]);
2942  }
2943  if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
2944    setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
2945  }
2946
2947  for (uInt i = 0; i < rejectNWaves.size(); ++i) {
2948    tempRejectNWaves.push_back(rejectNWaves[i]);
2949  }
2950  if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
2951    setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
2952  }
2953
2954  for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
2955    bool found = false;
2956    for (uInt j = 0; j < nWaves.size(); ++j) {
2957      if (nWaves[j] == tempAddNWaves[i]) {
2958        found = true;
2959        break;
2960      }
2961    }
2962    if (!found) nWaves.push_back(tempAddNWaves[i]);
2963  }
2964
2965  for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
2966    for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
2967      if (*j == tempRejectNWaves[i]) {
2968        j = nWaves.erase(j);
2969      } else {
2970        ++j;
2971      }
2972    }
2973  }
2974
2975  if (nWaves.size() > 1) {
2976    sort(nWaves.begin(), nWaves.end());
2977    unique(nWaves.begin(), nWaves.end());
2978  }
2979}
2980
2981void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
2982{
2983  if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
2984    int val = nWaves[0];
2985    int nyquistFreq = nchan(getIF(whichrow))/2+1;
2986    nWaves.clear();
2987    if (val > nyquistFreq) {  // for safety, at least nWaves contains a constant; CAS-3759
2988      nWaves.push_back(0);
2989    }
2990    while (val <= nyquistFreq) {
2991      nWaves.push_back(val);
2992      val++;
2993    }
2994  }
2995}
2996
2997void Scantable::sinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2998{
2999  try {
3000    ofstream ofs;
3001    String coordInfo = "";
3002    bool hasSameNchan = true;
3003    bool outTextFile = false;
3004
3005    if (blfile != "") {
3006      ofs.open(blfile.c_str(), ios::out | ios::app);
3007      if (ofs) outTextFile = true;
3008    }
3009
3010    if (outLogger || outTextFile) {
3011      coordInfo = getCoordInfo()[0];
3012      if (coordInfo == "") coordInfo = "channel";
3013      hasSameNchan = hasSameNchanOverIFs();
3014    }
3015
3016    //Fitter fitter = Fitter();
3017    //fitter.setExpression("sinusoid", nWaves);
3018    //fitter.setIterClipping(thresClip, nIterClip);
3019
3020    int nRow = nrow();
3021    std::vector<bool> chanMask;
3022    std::vector<int> nWaves;
3023
3024    bool showProgress;
3025    int minNRow;
3026    parseProgressInfo(progressInfo, showProgress, minNRow);
3027
3028    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3029      chanMask = getCompositeChanMask(whichrow, mask);
3030      selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
3031
3032      //FOR DEBUGGING------------
3033      /*
3034      if (whichrow < 0) {// == nRow -1) {
3035        cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
3036        if (applyFFT) {
3037          cout << "[ ";
3038          for (uInt j = 0; j < nWaves.size(); ++j) {
3039            cout << nWaves[j] << ", ";
3040          }
3041          cout << " ]    " << endl;
3042        }
3043        cout << flush;
3044      }
3045      */
3046      //-------------------------
3047
3048      //fitBaseline(chanMask, whichrow, fitter);
3049      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
3050      std::vector<float> params;
3051      int nClipped = 0;
3052      std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
3053      setSpectrum(res, whichrow);
3054      //
3055
3056      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
3057      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3058    }
3059
3060    if (outTextFile) ofs.close();
3061
3062  } catch (...) {
3063    throw;
3064  }
3065}
3066
3067void Scantable::autoSinusoidBaseline(const std::vector<bool>& mask, const bool applyFFT, const std::string& fftMethod, const std::string& fftThresh, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, float thresClip, int nIterClip, const std::vector<int>& edge, float threshold, int chanAvgLimit, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
3068{
3069  try {
3070    ofstream ofs;
3071    String coordInfo = "";
3072    bool hasSameNchan = true;
3073    bool outTextFile = false;
3074
3075    if (blfile != "") {
3076      ofs.open(blfile.c_str(), ios::out | ios::app);
3077      if (ofs) outTextFile = true;
3078    }
3079
3080    if (outLogger || outTextFile) {
3081      coordInfo = getCoordInfo()[0];
3082      if (coordInfo == "") coordInfo = "channel";
3083      hasSameNchan = hasSameNchanOverIFs();
3084    }
3085
3086    //Fitter fitter = Fitter();
3087    //fitter.setExpression("sinusoid", nWaves);
3088    //fitter.setIterClipping(thresClip, nIterClip);
3089
3090    int nRow = nrow();
3091    std::vector<bool> chanMask;
3092    std::vector<int> nWaves;
3093
3094    int minEdgeSize = getIFNos().size()*2;
3095    STLineFinder lineFinder = STLineFinder();
3096    lineFinder.setOptions(threshold, 3, chanAvgLimit);
3097
3098    bool showProgress;
3099    int minNRow;
3100    parseProgressInfo(progressInfo, showProgress, minNRow);
3101
3102    for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3103
3104      //-------------------------------------------------------
3105      //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
3106      //-------------------------------------------------------
3107      int edgeSize = edge.size();
3108      std::vector<int> currentEdge;
3109      if (edgeSize >= 2) {
3110        int idx = 0;
3111        if (edgeSize > 2) {
3112          if (edgeSize < minEdgeSize) {
3113            throw(AipsError("Length of edge element info is less than that of IFs"));
3114          }
3115          idx = 2 * getIF(whichrow);
3116        }
3117        currentEdge.push_back(edge[idx]);
3118        currentEdge.push_back(edge[idx+1]);
3119      } else {
3120        throw(AipsError("Wrong length of edge element"));
3121      }
3122      lineFinder.setData(getSpectrum(whichrow));
3123      lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
3124      chanMask = lineFinder.getMask();
3125      //-------------------------------------------------------
3126
3127      selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
3128
3129      //fitBaseline(chanMask, whichrow, fitter);
3130      //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
3131      std::vector<float> params;
3132      int nClipped = 0;
3133      std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
3134      setSpectrum(res, whichrow);
3135      //
3136
3137      outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
3138      showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3139    }
3140
3141    if (outTextFile) ofs.close();
3142
3143  } catch (...) {
3144    throw;
3145  }
3146}
3147
3148std::vector<float> Scantable::doSinusoidFitting(const std::vector<float>& data, const std::vector<bool>& mask, const std::vector<int>& waveNumbers, std::vector<float>& params, int& nClipped, float thresClip, int nIterClip, bool getResidual)
3149{
3150  if (data.size() != mask.size()) {
3151    throw(AipsError("data and mask sizes are not identical"));
3152  }
3153  if (data.size() < 2) {
3154    throw(AipsError("data size is too short"));
3155  }
3156  if (waveNumbers.size() == 0) {
3157    throw(AipsError("no wave numbers given"));
3158  }
3159  std::vector<int> nWaves;  // sorted and uniqued array of wave numbers
3160  nWaves.reserve(waveNumbers.size());
3161  copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
3162  sort(nWaves.begin(), nWaves.end());
3163  std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
3164  nWaves.erase(end_it, nWaves.end());
3165
3166  int minNWaves = nWaves[0];
3167  if (minNWaves < 0) {
3168    throw(AipsError("wave number must be positive or zero (i.e. constant)"));
3169  }
3170  bool hasConstantTerm = (minNWaves == 0);
3171
3172  int nChan = data.size();
3173  std::vector<int> maskArray;
3174  std::vector<int> x;
3175  for (int i = 0; i < nChan; ++i) {
3176    maskArray.push_back(mask[i] ? 1 : 0);
3177    if (mask[i]) {
3178      x.push_back(i);
3179    }
3180  }
3181
3182  int initNData = x.size();
3183
3184  int nData = initNData;
3185  int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0);  //number of parameters to solve.
3186
3187  const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
3188  double baseXFactor = 2.0*PI/(double)(nChan-1);  //the denominator (nChan-1) should be changed to (xdata[nChan-1]-xdata[0]) for accepting x-values given in velocity or frequency when this function is moved to fitter. (2011/03/30 WK)
3189
3190  // xArray : contains elemental values for computing the least-square matrix.
3191  //          xArray.size() is nDOF and xArray[*].size() is nChan.
3192  //          Each xArray element are as follows:
3193  //          xArray[0]    = {1.0, 1.0, 1.0, ..., 1.0},
3194  //          xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
3195  //          xArray[2n]   = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
3196  //          where (1 <= n <= nMaxWavesInSW),
3197  //          or,
3198  //          xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
3199  //          xArray[2n]   = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
3200  //          where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
3201  std::vector<std::vector<double> > xArray;
3202  if (hasConstantTerm) {
3203    std::vector<double> xu;
3204    for (int j = 0; j < nChan; ++j) {
3205      xu.push_back(1.0);
3206    }
3207    xArray.push_back(xu);
3208  }
3209  for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
3210    double xFactor = baseXFactor*(double)nWaves[i];
3211    std::vector<double> xs, xc;
3212    xs.clear();
3213    xc.clear();
3214    for (int j = 0; j < nChan; ++j) {
3215      xs.push_back(sin(xFactor*(double)j));
3216      xc.push_back(cos(xFactor*(double)j));
3217    }
3218    xArray.push_back(xs);
3219    xArray.push_back(xc);
3220  }
3221
3222  std::vector<double> z1, r1, residual;
3223  for (int i = 0; i < nChan; ++i) {
3224    z1.push_back((double)data[i]);
3225    r1.push_back(0.0);
3226    residual.push_back(0.0);
3227  }
3228
3229  for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
3230    // xMatrix : horizontal concatenation of
3231    //           the least-sq. matrix (left) and an
3232    //           identity matrix (right).
3233    // the right part is used to calculate the inverse matrix of the left part.
3234    double xMatrix[nDOF][2*nDOF];
3235    double zMatrix[nDOF];
3236    for (int i = 0; i < nDOF; ++i) {
3237      for (int j = 0; j < 2*nDOF; ++j) {
3238        xMatrix[i][j] = 0.0;
3239      }
3240      xMatrix[i][nDOF+i] = 1.0;
3241      zMatrix[i] = 0.0;
3242    }
3243
3244    int nUseData = 0;
3245    for (int k = 0; k < nChan; ++k) {
3246      if (maskArray[k] == 0) continue;
3247
3248      for (int i = 0; i < nDOF; ++i) {
3249        for (int j = i; j < nDOF; ++j) {
3250          xMatrix[i][j] += xArray[i][k] * xArray[j][k];
3251        }
3252        zMatrix[i] += z1[k] * xArray[i][k];
3253      }
3254
3255      nUseData++;
3256    }
3257
3258    if (nUseData < 1) {
3259        throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));     
3260    }
3261
3262    for (int i = 0; i < nDOF; ++i) {
3263      for (int j = 0; j < i; ++j) {
3264        xMatrix[i][j] = xMatrix[j][i];
3265      }
3266    }
3267
3268    std::vector<double> invDiag;
3269    for (int i = 0; i < nDOF; ++i) {
3270      invDiag.push_back(1.0/xMatrix[i][i]);
3271      for (int j = 0; j < nDOF; ++j) {
3272        xMatrix[i][j] *= invDiag[i];
3273      }
3274    }
3275
3276    for (int k = 0; k < nDOF; ++k) {
3277      for (int i = 0; i < nDOF; ++i) {
3278        if (i != k) {
3279          double factor1 = xMatrix[k][k];
3280          double factor2 = xMatrix[i][k];
3281          for (int j = k; j < 2*nDOF; ++j) {
3282            xMatrix[i][j] *= factor1;
3283            xMatrix[i][j] -= xMatrix[k][j]*factor2;
3284            xMatrix[i][j] /= factor1;
3285          }
3286        }
3287      }
3288      double xDiag = xMatrix[k][k];
3289      for (int j = k; j < 2*nDOF; ++j) {
3290        xMatrix[k][j] /= xDiag;
3291      }
3292    }
3293   
3294    for (int i = 0; i < nDOF; ++i) {
3295      for (int j = 0; j < nDOF; ++j) {
3296        xMatrix[i][nDOF+j] *= invDiag[j];
3297      }
3298    }
3299    //compute a vector y which consists of the coefficients of the sinusoids forming the
3300    //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
3301    //and cosine functions, respectively.
3302    std::vector<double> y;
3303    params.clear();
3304    for (int i = 0; i < nDOF; ++i) {
3305      y.push_back(0.0);
3306      for (int j = 0; j < nDOF; ++j) {
3307        y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
3308      }
3309      params.push_back(y[i]);
3310    }
3311
3312    for (int i = 0; i < nChan; ++i) {
3313      r1[i] = y[0];
3314      for (int j = 1; j < nDOF; ++j) {
3315        r1[i] += y[j]*xArray[j][i];
3316      }
3317      residual[i] = z1[i] - r1[i];
3318    }
3319
3320    if ((nClip == nIterClip) || (thresClip <= 0.0)) {
3321      break;
3322    } else {
3323      double stdDev = 0.0;
3324      for (int i = 0; i < nChan; ++i) {
3325        stdDev += residual[i]*residual[i]*(double)maskArray[i];
3326      }
3327      stdDev = sqrt(stdDev/(double)nData);
3328     
3329      double thres = stdDev * thresClip;
3330      int newNData = 0;
3331      for (int i = 0; i < nChan; ++i) {
3332        if (abs(residual[i]) >= thres) {
3333          maskArray[i] = 0;
3334        }
3335        if (maskArray[i] > 0) {
3336          newNData++;
3337        }
3338      }
3339      if (newNData == nData) {
3340        break; //no more flag to add. iteration stops.
3341      } else {
3342        nData = newNData;
3343      }
3344    }
3345  }
3346
3347  nClipped = initNData - nData;
3348
3349  std::vector<float> result;
3350  if (getResidual) {
3351    for (int i = 0; i < nChan; ++i) {
3352      result.push_back((float)residual[i]);
3353    }
3354  } else {
3355    for (int i = 0; i < nChan; ++i) {
3356      result.push_back((float)r1[i]);
3357    }
3358  }
3359
3360  return result;
3361}
3362
3363void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
3364{
3365  std::vector<double> dAbcissa = getAbcissa(whichrow);
3366  std::vector<float> abcissa;
3367  for (uInt i = 0; i < dAbcissa.size(); ++i) {
3368    abcissa.push_back((float)dAbcissa[i]);
3369  }
3370  std::vector<float> spec = getSpectrum(whichrow);
3371
3372  fitter.setData(abcissa, spec, mask);
3373  fitter.lfit();
3374}
3375
3376std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
3377{
3378  std::vector<bool> mask = getMask(whichrow);
3379  uInt maskSize = mask.size();
3380  if (inMask.size() != 0) {
3381    if (maskSize != inMask.size()) {
3382      throw(AipsError("mask sizes are not the same."));
3383    }
3384    for (uInt i = 0; i < maskSize; ++i) {
3385      mask[i] = mask[i] && inMask[i];
3386    }
3387  }
3388
3389  return mask;
3390}
3391
3392/*
3393std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
3394{
3395  int edgeSize = edge.size();
3396  std::vector<int> currentEdge;
3397  if (edgeSize >= 2) {
3398      int idx = 0;
3399      if (edgeSize > 2) {
3400        if (edgeSize < minEdgeSize) {
3401          throw(AipsError("Length of edge element info is less than that of IFs"));
3402        }
3403        idx = 2 * getIF(whichrow);
3404      }
3405      currentEdge.push_back(edge[idx]);
3406      currentEdge.push_back(edge[idx+1]);
3407  } else {
3408    throw(AipsError("Wrong length of edge element"));
3409  }
3410
3411  lineFinder.setData(getSpectrum(whichrow));
3412  lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
3413
3414  return lineFinder.getMask();
3415}
3416*/
3417
3418/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK)  */
3419void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, Fitter& fitter)
3420{
3421  if (outLogger || outTextFile) {
3422    std::vector<float> params = fitter.getParameters();
3423    std::vector<bool>  fixed  = fitter.getFixedParameters();
3424    float rms = getRms(chanMask, whichrow);
3425    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3426
3427    if (outLogger) {
3428      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3429      ols << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, false) << LogIO::POST ;
3430    }
3431    if (outTextFile) {
3432      ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true) << flush;
3433    }
3434  }
3435}
3436
3437/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
3438void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<int>& edge, const std::vector<float>& params, const int nClipped)
3439{
3440  if (outLogger || outTextFile) {
3441    float rms = getRms(chanMask, whichrow);
3442    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3443    std::vector<bool> fixed;
3444    fixed.clear();
3445
3446    if (outLogger) {
3447      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3448      ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3449    }
3450    if (outTextFile) {
3451      ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3452    }
3453  }
3454}
3455
3456/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
3457void Scantable::outputFittingResult(bool outLogger, bool outTextFile, const std::vector<bool>& chanMask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, ofstream& ofs, const casa::String& funcName, const std::vector<float>& params, const int nClipped)
3458{
3459  if (outLogger || outTextFile) {
3460    float rms = getRms(chanMask, whichrow);
3461    String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3462    std::vector<bool> fixed;
3463    fixed.clear();
3464
3465    if (outLogger) {
3466      LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3467      ols << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3468    }
3469    if (outTextFile) {
3470      ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3471    }
3472  }
3473}
3474
3475void Scantable::parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow)
3476{
3477  int idxDelimiter = progressInfo.find(",");
3478  if (idxDelimiter < 0) {
3479    throw(AipsError("wrong value in 'showprogress' parameter")) ;
3480  }
3481  showProgress = (progressInfo.substr(0, idxDelimiter) == "true");
3482  std::istringstream is(progressInfo.substr(idxDelimiter+1));
3483  is >> minNRow;
3484}
3485
3486void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress, const int nTotalThreshold)
3487{
3488  if (showProgress && (nTotal >= nTotalThreshold)) {
3489    int nInterval = int(floor(double(nTotal)/100.0));
3490    if (nInterval == 0) nInterval++;
3491
3492    if (nProcessed % nInterval == 0) {
3493      printf("\r");                          //go to the head of line
3494      printf("\x1b[31m\x1b[1m");             //set red color, highlighted
3495      printf("[%3d%%]", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
3496      printf("\x1b[39m\x1b[0m");             //set default attributes
3497      fflush(NULL);
3498    }
3499
3500    if (nProcessed == nTotal - 1) {
3501      printf("\r\x1b[K");                    //clear
3502      fflush(NULL);
3503    }
3504
3505  }
3506}
3507
3508std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
3509{
3510  std::vector<bool>  mask = getMask(whichrow);
3511
3512  if (inMask.size() > 0) {
3513    uInt maskSize = mask.size();
3514    if (maskSize != inMask.size()) {
3515      throw(AipsError("mask sizes are not the same."));
3516    }
3517    for (uInt i = 0; i < maskSize; ++i) {
3518      mask[i] = mask[i] && inMask[i];
3519    }
3520  }
3521
3522  Vector<Float> spec = getSpectrum(whichrow);
3523  mathutil::doZeroOrderInterpolation(spec, mask);
3524
3525  FFTServer<Float,Complex> ffts;
3526  Vector<Complex> fftres;
3527  ffts.fft0(fftres, spec);
3528
3529  std::vector<float> res;
3530  float norm = float(2.0/double(spec.size()));
3531
3532  if (getRealImag) {
3533    for (uInt i = 0; i < fftres.size(); ++i) {
3534      res.push_back(real(fftres[i])*norm);
3535      res.push_back(imag(fftres[i])*norm);
3536    }
3537  } else {
3538    for (uInt i = 0; i < fftres.size(); ++i) {
3539      res.push_back(abs(fftres[i])*norm);
3540      if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
3541    }
3542  }
3543
3544  return res;
3545}
3546
3547
3548float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
3549{
3550  Vector<Float> spec;
3551  specCol_.get(whichrow, spec);
3552
3553  float mean = 0.0;
3554  float smean = 0.0;
3555  int n = 0;
3556  for (uInt i = 0; i < spec.nelements(); ++i) {
3557    if (mask[i]) {
3558      mean += spec[i];
3559      smean += spec[i]*spec[i];
3560      n++;
3561    }
3562  }
3563
3564  mean /= (float)n;
3565  smean /= (float)n;
3566
3567  return sqrt(smean - mean*mean);
3568}
3569
3570
3571std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
3572{
3573  ostringstream oss;
3574
3575  if (verbose) {
3576    oss <<  " Scan[" << getScan(whichrow)  << "]";
3577    oss <<  " Beam[" << getBeam(whichrow)  << "]";
3578    oss <<    " IF[" << getIF(whichrow)    << "]";
3579    oss <<   " Pol[" << getPol(whichrow)   << "]";
3580    oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
3581    oss << "Fitter range = " << masklist << endl;
3582    oss << "Baseline parameters" << endl;
3583    oss << flush;
3584  }
3585
3586  return String(oss);
3587}
3588
3589std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const
3590{
3591  ostringstream oss;
3592
3593  if (verbose) {
3594    oss << "Results of baseline fit" << endl;
3595    oss << "  rms = " << setprecision(6) << rms << endl;
3596    if (nClipped >= 0) {
3597      oss << "  Number of clipped channels = " << nClipped << endl;
3598    }
3599    for (int i = 0; i < 60; ++i) {
3600      oss << "-";
3601    }
3602    oss << endl;
3603    oss << flush;
3604  }
3605
3606  return String(oss);
3607}
3608
3609std::string Scantable::formatBaselineParams(const std::vector<float>& params,
3610                                            const std::vector<bool>& fixed,
3611                                            float rms,
3612                                            int nClipped,
3613                                            const std::string& masklist,
3614                                            int whichrow,
3615                                            bool verbose,
3616                                            int start, int count,
3617                                            bool resetparamid) const
3618{
3619  int nParam = (int)(params.size());
3620
3621  if (nParam < 1) {
3622    return("  Not fitted");
3623  } else {
3624
3625    ostringstream oss;
3626    oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3627
3628    if (start < 0) start = 0;
3629    if (count < 0) count = nParam;
3630    int end = start + count;
3631    if (end > nParam) end = nParam;
3632    int paramidoffset = (resetparamid) ? (-start) : 0;
3633
3634    for (int i = start; i < end; ++i) {
3635      if (i > start) {
3636        oss << ",";
3637      }
3638      std::string sFix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
3639      oss << "  p" << (i+paramidoffset) << sFix << "= " << right << setw(13) << setprecision(6) << params[i];
3640    }
3641
3642    oss << endl;
3643    oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3644
3645    return String(oss);
3646  }
3647
3648}
3649
3650  std::string Scantable::formatPiecewiseBaselineParams(const std::vector<int>& ranges, const std::vector<float>& params, const std::vector<bool>& fixed, float rms, int nClipped, const std::string& masklist, int whichrow, bool verbose) const
3651{
3652  int nOutParam = (int)(params.size());
3653  int nPiece = (int)(ranges.size()) - 1;
3654
3655  if (nOutParam < 1) {
3656    return("  Not fitted");
3657  } else if (nPiece < 0) {
3658    return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose);
3659  } else if (nPiece < 1) {
3660    return("  Bad count of the piece edge info");
3661  } else if (nOutParam % nPiece != 0) {
3662    return("  Bad count of the output baseline parameters");
3663  } else {
3664
3665    int nParam = nOutParam / nPiece;
3666
3667    ostringstream oss;
3668    oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3669
3670    stringstream ss;
3671    ss << ranges[nPiece] << flush;
3672    int wRange = ss.str().size() * 2 + 5;
3673
3674    for (int i = 0; i < nPiece; ++i) {
3675      ss.str("");
3676      ss << "  [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
3677      oss << left << setw(wRange) << ss.str();
3678      oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, i*nParam, nParam, true);
3679    }
3680
3681    oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3682
3683    return String(oss);
3684  }
3685
3686}
3687
3688bool Scantable::hasSameNchanOverIFs()
3689{
3690  int nIF = nif(-1);
3691  int nCh;
3692  int totalPositiveNChan = 0;
3693  int nPositiveNChan = 0;
3694
3695  for (int i = 0; i < nIF; ++i) {
3696    nCh = nchan(i);
3697    if (nCh > 0) {
3698      totalPositiveNChan += nCh;
3699      nPositiveNChan++;
3700    }
3701  }
3702
3703  return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
3704}
3705
3706std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
3707{
3708  if (mask.size() <= 0) {
3709    throw(AipsError("The mask elements should be > 0"));
3710  }
3711  int IF = getIF(whichrow);
3712  if (mask.size() != (uInt)nchan(IF)) {
3713    throw(AipsError("Number of channels in scantable != number of mask elements"));
3714  }
3715
3716  if (verbose) {
3717    LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
3718    logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
3719    if (!hasSameNchan) {
3720      logOs << endl << "This mask is only valid for IF=" << IF;
3721    }
3722    logOs << LogIO::POST;
3723  }
3724
3725  std::vector<double> abcissa = getAbcissa(whichrow);
3726  std::vector<int> edge = getMaskEdgeIndices(mask);
3727
3728  ostringstream oss;
3729  oss.setf(ios::fixed);
3730  oss << setprecision(1) << "[";
3731  for (uInt i = 0; i < edge.size(); i+=2) {
3732    if (i > 0) oss << ",";
3733    oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
3734  }
3735  oss << "]" << flush;
3736
3737  return String(oss);
3738}
3739
3740std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
3741{
3742  if (mask.size() <= 0) {
3743    throw(AipsError("The mask elements should be > 0"));
3744  }
3745
3746  std::vector<int> out, startIndices, endIndices;
3747  int maskSize = mask.size();
3748
3749  startIndices.clear();
3750  endIndices.clear();
3751
3752  if (mask[0]) {
3753    startIndices.push_back(0);
3754  }
3755  for (int i = 1; i < maskSize; ++i) {
3756    if ((!mask[i-1]) && mask[i]) {
3757      startIndices.push_back(i);
3758    } else if (mask[i-1] && (!mask[i])) {
3759      endIndices.push_back(i-1);
3760    }
3761  }
3762  if (mask[maskSize-1]) {
3763    endIndices.push_back(maskSize-1);
3764  }
3765
3766  if (startIndices.size() != endIndices.size()) {
3767    throw(AipsError("Inconsistent Mask Size: bad data?"));
3768  }
3769  for (uInt i = 0; i < startIndices.size(); ++i) {
3770    if (startIndices[i] > endIndices[i]) {
3771      throw(AipsError("Mask start index > mask end index"));
3772    }
3773  }
3774
3775  out.clear();
3776  for (uInt i = 0; i < startIndices.size(); ++i) {
3777    out.push_back(startIndices[i]);
3778    out.push_back(endIndices[i]);
3779  }
3780
3781  return out;
3782}
3783
3784vector<float> Scantable::getTsysSpectrum( int whichrow ) const
3785{
3786  Vector<Float> tsys( tsysCol_(whichrow) ) ;
3787  vector<float> stlTsys ;
3788  tsys.tovector( stlTsys ) ;
3789  return stlTsys ;
3790}
3791
3792
3793}
3794//namespace asap
Note: See TracBrowser for help on using the repository browser.