source: trunk/src/Scantable.cpp@ 2537

Last change on this file since 2537 was 2531, checked in by Malte Marquarding, 12 years ago

Issue #265: Support multiple chains in IF in summary

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 116.2 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];
1278 if (tmplen > 1) {
1279 oss << " (" << tmplen << " chains)";
1280 }
1281 oss << endl;
1282 }
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 Vector<Float> oldtsys = tsysCol_( irow ) ;
1883 uInt newsize = nmax - nmin + 1 ;
1884 Slice slice( nmin, newsize, 1 ) ;
1885 specCol_.put( irow, oldspec( slice ) ) ;
1886 flagsCol_.put( irow, oldflag( slice ) ) ;
1887 if ( oldspec.size() == oldtsys.size() )
1888 tsysCol_.put( irow, oldtsys( slice ) ) ;
1889
1890 return ;
1891}
1892
1893void asap::Scantable::regridSpecChannel( double dnu, int nChan )
1894{
1895 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1896 os << "Regrid abcissa with spectral resoultion " << dnu << " " << freqTable_.getUnitString() << " with channel number " << ((nChan>0)? String(nChan) : "covering band width")<< LogIO::POST ;
1897 int freqnrow = freqTable_.table().nrow() ;
1898 Vector<bool> firstTime( freqnrow, true ) ;
1899 double oldincr, factor;
1900 uInt currId;
1901 Double refpix ;
1902 Double refval ;
1903 Double increment ;
1904 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1905 currId = mfreqidCol_(irow);
1906 vector<double> abcissa = getAbcissa( irow ) ;
1907 if (nChan < 0) {
1908 int oldsize = abcissa.size() ;
1909 double bw = (abcissa[oldsize-1]-abcissa[0]) + \
1910 0.5 * (abcissa[1]-abcissa[0] + abcissa[oldsize-1]-abcissa[oldsize-2]) ;
1911 nChan = int( ceil( abs(bw/dnu) ) ) ;
1912 }
1913 // actual regridding
1914 regridChannel( nChan, dnu, irow ) ;
1915
1916 // update FREQUENCIES subtable
1917 if (firstTime[currId]) {
1918 oldincr = abcissa[1]-abcissa[0] ;
1919 factor = dnu/oldincr ;
1920 firstTime[currId] = false ;
1921 freqTable_.getEntry( refpix, refval, increment, currId ) ;
1922
1923 //refval = refval - ( refpix + 0.5 * (1 - factor) ) * increment ;
1924 if (factor > 0 ) {
1925 refpix = (refpix + 0.5)/factor - 0.5;
1926 } else {
1927 refpix = (abcissa.size() - 0.5 - refpix)/abs(factor) - 0.5;
1928 }
1929 freqTable_.setEntry( refpix, refval, increment*factor, currId ) ;
1930 //os << "ID" << currId << ": channel width (Orig) = " << oldincr << " [" << freqTable_.getUnitString() << "], scale factor = " << factor << LogIO::POST ;
1931 //os << " frequency increment (Orig) = " << increment << "-> (New) " << increment*factor << LogIO::POST ;
1932 }
1933 }
1934}
1935
1936void asap::Scantable::regridChannel( int nChan, double dnu )
1937{
1938 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1939 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1940 // assumed that all rows have same nChan
1941 Vector<Float> arr = specCol_( 0 ) ;
1942 int oldsize = arr.nelements() ;
1943
1944 // if oldsize == nChan, nothing to do
1945 if ( oldsize == nChan ) {
1946 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1947 return ;
1948 }
1949
1950 // if oldChan < nChan, unphysical operation
1951 if ( oldsize < nChan ) {
1952 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1953 return ;
1954 }
1955
1956 // change channel number for specCol_, flagCol_, and tsysCol_ (if necessary)
1957 vector<string> coordinfo = getCoordInfo() ;
1958 string oldinfo = coordinfo[0] ;
1959 coordinfo[0] = "Hz" ;
1960 setCoordInfo( coordinfo ) ;
1961 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1962 regridChannel( nChan, dnu, irow ) ;
1963 }
1964 coordinfo[0] = oldinfo ;
1965 setCoordInfo( coordinfo ) ;
1966
1967
1968 // NOTE: this method does not update metadata such as
1969 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1970
1971 return ;
1972}
1973
1974void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1975{
1976 // logging
1977 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1978 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1979
1980 Vector<Float> oldspec = specCol_( irow ) ;
1981 Vector<uChar> oldflag = flagsCol_( irow ) ;
1982 Vector<Float> oldtsys = tsysCol_( irow ) ;
1983 Vector<Float> newspec( nChan, 0 ) ;
1984 Vector<uChar> newflag( nChan, true ) ;
1985 Vector<Float> newtsys ;
1986 bool regridTsys = false ;
1987 if (oldtsys.size() == oldspec.size()) {
1988 regridTsys = true ;
1989 newtsys.resize(nChan,false) ;
1990 newtsys = 0 ;
1991 }
1992
1993 // regrid
1994 vector<double> abcissa = getAbcissa( irow ) ;
1995 int oldsize = abcissa.size() ;
1996 double olddnu = abcissa[1] - abcissa[0] ;
1997 //int ichan = 0 ;
1998 double wsum = 0.0 ;
1999 Vector<double> zi( nChan+1 ) ;
2000 Vector<double> yi( oldsize + 1 ) ;
2001 yi[0] = abcissa[0] - 0.5 * olddnu ;
2002 for ( int ii = 1 ; ii < oldsize ; ii++ )
2003 yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
2004 yi[oldsize] = abcissa[oldsize-1] \
2005 + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
2006 //zi[0] = abcissa[0] - 0.5 * olddnu ;
2007 zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
2008 for ( int ii = 1 ; ii < nChan ; ii++ )
2009 zi[ii] = zi[0] + dnu * ii ;
2010 zi[nChan] = zi[nChan-1] + dnu ;
2011 // Access zi and yi in ascending order
2012 int izs = ((dnu > 0) ? 0 : nChan ) ;
2013 int ize = ((dnu > 0) ? nChan : 0 ) ;
2014 int izincr = ((dnu > 0) ? 1 : -1 ) ;
2015 int ichan = ((olddnu > 0) ? 0 : oldsize ) ;
2016 int iye = ((olddnu > 0) ? oldsize : 0 ) ;
2017 int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
2018 //for ( int ii = izs ; ii != ize ; ii+=izincr ){
2019 int ii = izs ;
2020 while (ii != ize) {
2021 // always zl < zr
2022 double zl = zi[ii] ;
2023 double zr = zi[ii+izincr] ;
2024 // Need to access smaller index for the new spec, flag, and tsys.
2025 // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
2026 int i = min(ii, ii+izincr) ;
2027 //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
2028 int jj = ichan ;
2029 while (jj != iye) {
2030 // always yl < yr
2031 double yl = yi[jj] ;
2032 double yr = yi[jj+iyincr] ;
2033 // Need to access smaller index for the original spec, flag, and tsys.
2034 // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
2035 int j = min(jj, jj+iyincr) ;
2036 if ( yr <= zl ) {
2037 jj += iyincr ;
2038 continue ;
2039 }
2040 else if ( yl <= zl ) {
2041 if ( yr < zr ) {
2042 if (!oldflag[j]) {
2043 newspec[i] += oldspec[j] * ( yr - zl ) ;
2044 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
2045 wsum += ( yr - zl ) ;
2046 }
2047 newflag[i] = newflag[i] && oldflag[j] ;
2048 }
2049 else {
2050 if (!oldflag[j]) {
2051 newspec[i] += oldspec[j] * abs(dnu) ;
2052 if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
2053 wsum += abs(dnu) ;
2054 }
2055 newflag[i] = newflag[i] && oldflag[j] ;
2056 ichan = jj ;
2057 break ;
2058 }
2059 }
2060 else if ( yl < zr ) {
2061 if ( yr <= zr ) {
2062 if (!oldflag[j]) {
2063 newspec[i] += oldspec[j] * ( yr - yl ) ;
2064 if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
2065 wsum += ( yr - yl ) ;
2066 }
2067 newflag[i] = newflag[i] && oldflag[j] ;
2068 }
2069 else {
2070 if (!oldflag[j]) {
2071 newspec[i] += oldspec[j] * ( zr - yl ) ;
2072 if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
2073 wsum += ( zr - yl ) ;
2074 }
2075 newflag[i] = newflag[i] && oldflag[j] ;
2076 ichan = jj ;
2077 break ;
2078 }
2079 }
2080 else {
2081 ichan = jj - iyincr ;
2082 break ;
2083 }
2084 jj += iyincr ;
2085 }
2086 if ( wsum != 0.0 ) {
2087 newspec[i] /= wsum ;
2088 if (regridTsys) newtsys[i] /= wsum ;
2089 }
2090 wsum = 0.0 ;
2091 ii += izincr ;
2092 }
2093// if ( dnu > 0.0 ) {
2094// for ( int ii = 0 ; ii < nChan ; ii++ ) {
2095// double zl = zi[ii] ;
2096// double zr = zi[ii+1] ;
2097// for ( int j = ichan ; j < oldsize ; j++ ) {
2098// double yl = yi[j] ;
2099// double yr = yi[j+1] ;
2100// if ( yl <= zl ) {
2101// if ( yr <= zl ) {
2102// continue ;
2103// }
2104// else if ( yr <= zr ) {
2105// if (!oldflag[j]) {
2106// newspec[ii] += oldspec[j] * ( yr - zl ) ;
2107// if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - zl ) ;
2108// wsum += ( yr - zl ) ;
2109// }
2110// newflag[ii] = newflag[ii] && oldflag[j] ;
2111// }
2112// else {
2113// if (!oldflag[j]) {
2114// newspec[ii] += oldspec[j] * dnu ;
2115// if (regridTsys) newtsys[ii] += oldtsys[j] * dnu ;
2116// wsum += dnu ;
2117// }
2118// newflag[ii] = newflag[ii] && oldflag[j] ;
2119// ichan = j ;
2120// break ;
2121// }
2122// }
2123// else if ( yl < zr ) {
2124// if ( yr <= zr ) {
2125// if (!oldflag[j]) {
2126// newspec[ii] += oldspec[j] * ( yr - yl ) ;
2127// if (regridTsys) newtsys[ii] += oldtsys[j] * ( yr - yl ) ;
2128// wsum += ( yr - yl ) ;
2129// }
2130// newflag[ii] = newflag[ii] && oldflag[j] ;
2131// }
2132// else {
2133// if (!oldflag[j]) {
2134// newspec[ii] += oldspec[j] * ( zr - yl ) ;
2135// if (regridTsys) newtsys[ii] += oldtsys[j] * ( zr - yl ) ;
2136// wsum += ( zr - yl ) ;
2137// }
2138// newflag[ii] = newflag[ii] && oldflag[j] ;
2139// ichan = j ;
2140// break ;
2141// }
2142// }
2143// else {
2144// ichan = j - 1 ;
2145// break ;
2146// }
2147// }
2148// if ( wsum != 0.0 ) {
2149// newspec[ii] /= wsum ;
2150// if (regridTsys) newtsys[ii] /= wsum ;
2151// }
2152// wsum = 0.0 ;
2153// }
2154// }
2155// else if ( dnu < 0.0 ) {
2156// for ( int ii = 0 ; ii < nChan ; ii++ ) {
2157// double zl = zi[ii] ;
2158// double zr = zi[ii+1] ;
2159// for ( int j = ichan ; j < oldsize ; j++ ) {
2160// double yl = yi[j] ;
2161// double yr = yi[j+1] ;
2162// if ( yl >= zl ) {
2163// if ( yr >= zl ) {
2164// continue ;
2165// }
2166// else if ( yr >= zr ) {
2167// if (!oldflag[j]) {
2168// newspec[ii] += oldspec[j] * abs( yr - zl ) ;
2169// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - zl ) ;
2170// wsum += abs( yr - zl ) ;
2171// }
2172// newflag[ii] = newflag[ii] && oldflag[j] ;
2173// }
2174// else {
2175// if (!oldflag[j]) {
2176// newspec[ii] += oldspec[j] * abs( dnu ) ;
2177// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( dnu ) ;
2178// wsum += abs( dnu ) ;
2179// }
2180// newflag[ii] = newflag[ii] && oldflag[j] ;
2181// ichan = j ;
2182// break ;
2183// }
2184// }
2185// else if ( yl > zr ) {
2186// if ( yr >= zr ) {
2187// if (!oldflag[j]) {
2188// newspec[ii] += oldspec[j] * abs( yr - yl ) ;
2189// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( yr - yl ) ;
2190// wsum += abs( yr - yl ) ;
2191// }
2192// newflag[ii] = newflag[ii] && oldflag[j] ;
2193// }
2194// else {
2195// if (!oldflag[j]) {
2196// newspec[ii] += oldspec[j] * abs( zr - yl ) ;
2197// if (regridTsys) newtsys[ii] += oldtsys[j] * abs( zr - yl ) ;
2198// wsum += abs( zr - yl ) ;
2199// }
2200// newflag[ii] = newflag[ii] && oldflag[j] ;
2201// ichan = j ;
2202// break ;
2203// }
2204// }
2205// else {
2206// ichan = j - 1 ;
2207// break ;
2208// }
2209// }
2210// if ( wsum != 0.0 ) {
2211// newspec[ii] /= wsum ;
2212// if (regridTsys) newtsys[ii] /= wsum ;
2213// }
2214// wsum = 0.0 ;
2215// }
2216// }
2217// // //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
2218// // pile += dnu ;
2219// // wedge = olddnu * ( refChan + 1 ) ;
2220// // while ( wedge < pile ) {
2221// // newspec[0] += olddnu * oldspec[refChan] ;
2222// // newflag[0] = newflag[0] || oldflag[refChan] ;
2223// // //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
2224// // refChan++ ;
2225// // wedge += olddnu ;
2226// // wsum += olddnu ;
2227// // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2228// // }
2229// // frac = ( wedge - pile ) / olddnu ;
2230// // wsum += ( 1.0 - frac ) * olddnu ;
2231// // newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2232// // newflag[0] = newflag[0] || oldflag[refChan] ;
2233// // //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
2234// // //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
2235// // newspec[0] /= wsum ;
2236// // //ofs << "newspec[0] = " << newspec[0] << endl ;
2237// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2238
2239// // /***
2240// // * ichan = 1 - nChan-2
2241// // ***/
2242// // for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
2243// // pile += dnu ;
2244// // newspec[ichan] += frac * olddnu * oldspec[refChan] ;
2245// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2246// // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
2247// // refChan++ ;
2248// // wedge += olddnu ;
2249// // wsum = frac * olddnu ;
2250// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2251// // while ( wedge < pile ) {
2252// // newspec[ichan] += olddnu * oldspec[refChan] ;
2253// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2254// // //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
2255// // refChan++ ;
2256// // wedge += olddnu ;
2257// // wsum += olddnu ;
2258// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2259// // }
2260// // frac = ( wedge - pile ) / olddnu ;
2261// // wsum += ( 1.0 - frac ) * olddnu ;
2262// // newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
2263// // newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
2264// // //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
2265// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2266// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
2267// // newspec[ichan] /= wsum ;
2268// // //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
2269// // }
2270
2271// // /***
2272// // * ichan = nChan-1
2273// // ***/
2274// // // NOTE: Assumed that all spectra have the same bandwidth
2275// // pile += dnu ;
2276// // newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
2277// // newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
2278// // //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2279// // refChan++ ;
2280// // wedge += olddnu ;
2281// // wsum = frac * olddnu ;
2282// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2283// // for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
2284// // newspec[nChan-1] += olddnu * oldspec[jchan] ;
2285// // newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
2286// // wsum += olddnu ;
2287// // //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
2288// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2289// // }
2290// // //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
2291// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
2292// // newspec[nChan-1] /= wsum ;
2293// // //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
2294
2295// // // ofs.close() ;
2296
2297 specCol_.put( irow, newspec ) ;
2298 flagsCol_.put( irow, newflag ) ;
2299 if (regridTsys) tsysCol_.put( irow, newtsys );
2300
2301 return ;
2302}
2303
2304std::vector<float> Scantable::getWeather(int whichrow) const
2305{
2306 std::vector<float> out(5);
2307 //Float temperature, pressure, humidity, windspeed, windaz;
2308 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
2309 mweatheridCol_(uInt(whichrow)));
2310
2311
2312 return out;
2313}
2314
2315bool Scantable::getFlagtraFast(uInt whichrow)
2316{
2317 uChar flag;
2318 Vector<uChar> flags;
2319 flagsCol_.get(whichrow, flags);
2320 flag = flags[0];
2321 for (uInt i = 1; i < flags.size(); ++i) {
2322 flag &= flags[i];
2323 }
2324 return ((flag >> 7) == 1);
2325}
2326
2327void Scantable::polyBaseline(const std::vector<bool>& mask, int order, bool getResidual, const std::string& progressInfo, const bool outLogger, const std::string& blfile)
2328{
2329 try {
2330 ofstream ofs;
2331 String coordInfo = "";
2332 bool hasSameNchan = true;
2333 bool outTextFile = false;
2334
2335 if (blfile != "") {
2336 ofs.open(blfile.c_str(), ios::out | ios::app);
2337 if (ofs) outTextFile = true;
2338 }
2339
2340 if (outLogger || outTextFile) {
2341 coordInfo = getCoordInfo()[0];
2342 if (coordInfo == "") coordInfo = "channel";
2343 hasSameNchan = hasSameNchanOverIFs();
2344 }
2345
2346 Fitter fitter = Fitter();
2347 fitter.setExpression("poly", order);
2348 //fitter.setIterClipping(thresClip, nIterClip);
2349
2350 int nRow = nrow();
2351 std::vector<bool> chanMask;
2352 bool showProgress;
2353 int minNRow;
2354 parseProgressInfo(progressInfo, showProgress, minNRow);
2355
2356 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2357 chanMask = getCompositeChanMask(whichrow, mask);
2358 fitBaseline(chanMask, whichrow, fitter);
2359 setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2360 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "polyBaseline()", fitter);
2361 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2362 }
2363
2364 if (outTextFile) ofs.close();
2365
2366 } catch (...) {
2367 throw;
2368 }
2369}
2370
2371void 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)
2372{
2373 try {
2374 ofstream ofs;
2375 String coordInfo = "";
2376 bool hasSameNchan = true;
2377 bool outTextFile = false;
2378
2379 if (blfile != "") {
2380 ofs.open(blfile.c_str(), ios::out | ios::app);
2381 if (ofs) outTextFile = true;
2382 }
2383
2384 if (outLogger || outTextFile) {
2385 coordInfo = getCoordInfo()[0];
2386 if (coordInfo == "") coordInfo = "channel";
2387 hasSameNchan = hasSameNchanOverIFs();
2388 }
2389
2390 Fitter fitter = Fitter();
2391 fitter.setExpression("poly", order);
2392 //fitter.setIterClipping(thresClip, nIterClip);
2393
2394 int nRow = nrow();
2395 std::vector<bool> chanMask;
2396 int minEdgeSize = getIFNos().size()*2;
2397 STLineFinder lineFinder = STLineFinder();
2398 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2399
2400 bool showProgress;
2401 int minNRow;
2402 parseProgressInfo(progressInfo, showProgress, minNRow);
2403
2404 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2405
2406 //-------------------------------------------------------
2407 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2408 //-------------------------------------------------------
2409 int edgeSize = edge.size();
2410 std::vector<int> currentEdge;
2411 if (edgeSize >= 2) {
2412 int idx = 0;
2413 if (edgeSize > 2) {
2414 if (edgeSize < minEdgeSize) {
2415 throw(AipsError("Length of edge element info is less than that of IFs"));
2416 }
2417 idx = 2 * getIF(whichrow);
2418 }
2419 currentEdge.push_back(edge[idx]);
2420 currentEdge.push_back(edge[idx+1]);
2421 } else {
2422 throw(AipsError("Wrong length of edge element"));
2423 }
2424 lineFinder.setData(getSpectrum(whichrow));
2425 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2426 chanMask = lineFinder.getMask();
2427 //-------------------------------------------------------
2428
2429 fitBaseline(chanMask, whichrow, fitter);
2430 setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2431
2432 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoPolyBaseline()", fitter);
2433 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2434 }
2435
2436 if (outTextFile) ofs.close();
2437
2438 } catch (...) {
2439 throw;
2440 }
2441}
2442
2443void 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)
2444{
2445 try {
2446 ofstream ofs;
2447 String coordInfo = "";
2448 bool hasSameNchan = true;
2449 bool outTextFile = false;
2450
2451 if (blfile != "") {
2452 ofs.open(blfile.c_str(), ios::out | ios::app);
2453 if (ofs) outTextFile = true;
2454 }
2455
2456 if (outLogger || outTextFile) {
2457 coordInfo = getCoordInfo()[0];
2458 if (coordInfo == "") coordInfo = "channel";
2459 hasSameNchan = hasSameNchanOverIFs();
2460 }
2461
2462 //Fitter fitter = Fitter();
2463 //fitter.setExpression("cspline", nPiece);
2464 //fitter.setIterClipping(thresClip, nIterClip);
2465
2466 bool showProgress;
2467 int minNRow;
2468 parseProgressInfo(progressInfo, showProgress, minNRow);
2469
2470 int nRow = nrow();
2471 std::vector<bool> chanMask;
2472
2473 //--------------------------------
2474 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2475 chanMask = getCompositeChanMask(whichrow, mask);
2476 //fitBaseline(chanMask, whichrow, fitter);
2477 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2478 std::vector<int> pieceEdges(nPiece+1);
2479 std::vector<float> params(nPiece*4);
2480 int nClipped = 0;
2481 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2482 setSpectrum(res, whichrow);
2483 //
2484
2485 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "cubicSplineBaseline()", pieceEdges, params, nClipped);
2486 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2487 }
2488 //--------------------------------
2489
2490 if (outTextFile) ofs.close();
2491
2492 } catch (...) {
2493 throw;
2494 }
2495}
2496
2497void 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)
2498{
2499 try {
2500 ofstream ofs;
2501 String coordInfo = "";
2502 bool hasSameNchan = true;
2503 bool outTextFile = false;
2504
2505 if (blfile != "") {
2506 ofs.open(blfile.c_str(), ios::out | ios::app);
2507 if (ofs) outTextFile = true;
2508 }
2509
2510 if (outLogger || outTextFile) {
2511 coordInfo = getCoordInfo()[0];
2512 if (coordInfo == "") coordInfo = "channel";
2513 hasSameNchan = hasSameNchanOverIFs();
2514 }
2515
2516 //Fitter fitter = Fitter();
2517 //fitter.setExpression("cspline", nPiece);
2518 //fitter.setIterClipping(thresClip, nIterClip);
2519
2520 int nRow = nrow();
2521 std::vector<bool> chanMask;
2522 int minEdgeSize = getIFNos().size()*2;
2523 STLineFinder lineFinder = STLineFinder();
2524 lineFinder.setOptions(threshold, 3, chanAvgLimit);
2525
2526 bool showProgress;
2527 int minNRow;
2528 parseProgressInfo(progressInfo, showProgress, minNRow);
2529
2530 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
2531
2532 //-------------------------------------------------------
2533 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
2534 //-------------------------------------------------------
2535 int edgeSize = edge.size();
2536 std::vector<int> currentEdge;
2537 if (edgeSize >= 2) {
2538 int idx = 0;
2539 if (edgeSize > 2) {
2540 if (edgeSize < minEdgeSize) {
2541 throw(AipsError("Length of edge element info is less than that of IFs"));
2542 }
2543 idx = 2 * getIF(whichrow);
2544 }
2545 currentEdge.push_back(edge[idx]);
2546 currentEdge.push_back(edge[idx+1]);
2547 } else {
2548 throw(AipsError("Wrong length of edge element"));
2549 }
2550 lineFinder.setData(getSpectrum(whichrow));
2551 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
2552 chanMask = lineFinder.getMask();
2553 //-------------------------------------------------------
2554
2555
2556 //fitBaseline(chanMask, whichrow, fitter);
2557 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
2558 std::vector<int> pieceEdges(nPiece+1);
2559 std::vector<float> params(nPiece*4);
2560 int nClipped = 0;
2561 std::vector<float> res = doCubicSplineFitting(getSpectrum(whichrow), chanMask, nPiece, pieceEdges, params, nClipped, thresClip, nIterClip, getResidual);
2562 setSpectrum(res, whichrow);
2563 //
2564
2565 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoCubicSplineBaseline()", pieceEdges, params, nClipped);
2566 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
2567 }
2568
2569 if (outTextFile) ofs.close();
2570
2571 } catch (...) {
2572 throw;
2573 }
2574}
2575
2576std::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)
2577{
2578 if (data.size() != mask.size()) {
2579 throw(AipsError("data and mask sizes are not identical"));
2580 }
2581 if (nPiece < 1) {
2582 throw(AipsError("number of the sections must be one or more"));
2583 }
2584
2585 int nChan = data.size();
2586 std::vector<int> maskArray(nChan);
2587 std::vector<int> x(nChan);
2588 int j = 0;
2589 for (int i = 0; i < nChan; ++i) {
2590 maskArray[i] = mask[i] ? 1 : 0;
2591 if (mask[i]) {
2592 x[j] = i;
2593 j++;
2594 }
2595 }
2596 int initNData = j;
2597
2598 if (initNData < nPiece) {
2599 throw(AipsError("too few non-flagged channels"));
2600 }
2601
2602 int nElement = (int)(floor(floor((double)(initNData/nPiece))+0.5));
2603 std::vector<double> invEdge(nPiece-1);
2604 idxEdge[0] = x[0];
2605 for (int i = 1; i < nPiece; ++i) {
2606 int valX = x[nElement*i];
2607 idxEdge[i] = valX;
2608 invEdge[i-1] = 1.0/(double)valX;
2609 }
2610 idxEdge[nPiece] = x[initNData-1]+1;
2611
2612 int nData = initNData;
2613 int nDOF = nPiece + 3; //number of parameters to solve, namely, 4+(nPiece-1).
2614
2615 std::vector<double> x1(nChan), x2(nChan), x3(nChan);
2616 std::vector<double> z1(nChan), x1z1(nChan), x2z1(nChan), x3z1(nChan);
2617 std::vector<double> r1(nChan), residual(nChan);
2618 for (int i = 0; i < nChan; ++i) {
2619 double di = (double)i;
2620 double dD = (double)data[i];
2621 x1[i] = di;
2622 x2[i] = di*di;
2623 x3[i] = di*di*di;
2624 z1[i] = dD;
2625 x1z1[i] = dD*di;
2626 x2z1[i] = dD*di*di;
2627 x3z1[i] = dD*di*di*di;
2628 r1[i] = 0.0;
2629 residual[i] = 0.0;
2630 }
2631
2632 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
2633 // xMatrix : horizontal concatenation of
2634 // the least-sq. matrix (left) and an
2635 // identity matrix (right).
2636 // the right part is used to calculate the inverse matrix of the left part.
2637 double xMatrix[nDOF][2*nDOF];
2638 double zMatrix[nDOF];
2639 for (int i = 0; i < nDOF; ++i) {
2640 for (int j = 0; j < 2*nDOF; ++j) {
2641 xMatrix[i][j] = 0.0;
2642 }
2643 xMatrix[i][nDOF+i] = 1.0;
2644 zMatrix[i] = 0.0;
2645 }
2646
2647 for (int n = 0; n < nPiece; ++n) {
2648 int nUseDataInPiece = 0;
2649 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2650
2651 if (maskArray[i] == 0) continue;
2652
2653 xMatrix[0][0] += 1.0;
2654 xMatrix[0][1] += x1[i];
2655 xMatrix[0][2] += x2[i];
2656 xMatrix[0][3] += x3[i];
2657 xMatrix[1][1] += x2[i];
2658 xMatrix[1][2] += x3[i];
2659 xMatrix[1][3] += x2[i]*x2[i];
2660 xMatrix[2][2] += x2[i]*x2[i];
2661 xMatrix[2][3] += x3[i]*x2[i];
2662 xMatrix[3][3] += x3[i]*x3[i];
2663 zMatrix[0] += z1[i];
2664 zMatrix[1] += x1z1[i];
2665 zMatrix[2] += x2z1[i];
2666 zMatrix[3] += x3z1[i];
2667
2668 for (int j = 0; j < n; ++j) {
2669 double q = 1.0 - x1[i]*invEdge[j];
2670 q = q*q*q;
2671 xMatrix[0][j+4] += q;
2672 xMatrix[1][j+4] += q*x1[i];
2673 xMatrix[2][j+4] += q*x2[i];
2674 xMatrix[3][j+4] += q*x3[i];
2675 for (int k = 0; k < j; ++k) {
2676 double r = 1.0 - x1[i]*invEdge[k];
2677 r = r*r*r;
2678 xMatrix[k+4][j+4] += r*q;
2679 }
2680 xMatrix[j+4][j+4] += q*q;
2681 zMatrix[j+4] += q*z1[i];
2682 }
2683
2684 nUseDataInPiece++;
2685 }
2686
2687 if (nUseDataInPiece < 1) {
2688 std::vector<string> suffixOfPieceNumber(4);
2689 suffixOfPieceNumber[0] = "th";
2690 suffixOfPieceNumber[1] = "st";
2691 suffixOfPieceNumber[2] = "nd";
2692 suffixOfPieceNumber[3] = "rd";
2693 int idxNoDataPiece = (n % 10 <= 3) ? n : 0;
2694 ostringstream oss;
2695 oss << "all channels clipped or masked in " << n << suffixOfPieceNumber[idxNoDataPiece];
2696 oss << " piece of the spectrum. can't execute fitting anymore.";
2697 throw(AipsError(String(oss)));
2698 }
2699 }
2700
2701 for (int i = 0; i < nDOF; ++i) {
2702 for (int j = 0; j < i; ++j) {
2703 xMatrix[i][j] = xMatrix[j][i];
2704 }
2705 }
2706
2707 std::vector<double> invDiag(nDOF);
2708 for (int i = 0; i < nDOF; ++i) {
2709 invDiag[i] = 1.0/xMatrix[i][i];
2710 for (int j = 0; j < nDOF; ++j) {
2711 xMatrix[i][j] *= invDiag[i];
2712 }
2713 }
2714
2715 for (int k = 0; k < nDOF; ++k) {
2716 for (int i = 0; i < nDOF; ++i) {
2717 if (i != k) {
2718 double factor1 = xMatrix[k][k];
2719 double factor2 = xMatrix[i][k];
2720 for (int j = k; j < 2*nDOF; ++j) {
2721 xMatrix[i][j] *= factor1;
2722 xMatrix[i][j] -= xMatrix[k][j]*factor2;
2723 xMatrix[i][j] /= factor1;
2724 }
2725 }
2726 }
2727 double xDiag = xMatrix[k][k];
2728 for (int j = k; j < 2*nDOF; ++j) {
2729 xMatrix[k][j] /= xDiag;
2730 }
2731 }
2732
2733 for (int i = 0; i < nDOF; ++i) {
2734 for (int j = 0; j < nDOF; ++j) {
2735 xMatrix[i][nDOF+j] *= invDiag[j];
2736 }
2737 }
2738 //compute a vector y which consists of the coefficients of the best-fit spline curves
2739 //(a0,a1,a2,a3(,b3,c3,...)), namely, the ones for the leftmost piece and the ones of
2740 //cubic terms for the other pieces (in case nPiece>1).
2741 std::vector<double> y(nDOF);
2742 for (int i = 0; i < nDOF; ++i) {
2743 y[i] = 0.0;
2744 for (int j = 0; j < nDOF; ++j) {
2745 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
2746 }
2747 }
2748
2749 double a0 = y[0];
2750 double a1 = y[1];
2751 double a2 = y[2];
2752 double a3 = y[3];
2753
2754 int j = 0;
2755 for (int n = 0; n < nPiece; ++n) {
2756 for (int i = idxEdge[n]; i < idxEdge[n+1]; ++i) {
2757 r1[i] = a0 + a1*x1[i] + a2*x2[i] + a3*x3[i];
2758 }
2759 params[j] = a0;
2760 params[j+1] = a1;
2761 params[j+2] = a2;
2762 params[j+3] = a3;
2763 j += 4;
2764
2765 if (n == nPiece-1) break;
2766
2767 double d = y[4+n];
2768 double iE = invEdge[n];
2769 a0 += d;
2770 a1 -= 3.0*d*iE;
2771 a2 += 3.0*d*iE*iE;
2772 a3 -= d*iE*iE*iE;
2773 }
2774
2775 //subtract constant value for masked regions at the edge of spectrum
2776 if (idxEdge[0] > 0) {
2777 int n = idxEdge[0];
2778 for (int i = 0; i < idxEdge[0]; ++i) {
2779 //--cubic extrapolate--
2780 //r1[i] = params[0] + params[1]*x1[i] + params[2]*x2[i] + params[3]*x3[i];
2781 //--linear extrapolate--
2782 //r1[i] = (r1[n+1] - r1[n])/(x1[n+1] - x1[n])*(x1[i] - x1[n]) + r1[n];
2783 //--constant--
2784 r1[i] = r1[n];
2785 }
2786 }
2787 if (idxEdge[nPiece] < nChan) {
2788 int n = idxEdge[nPiece]-1;
2789 for (int i = idxEdge[nPiece]; i < nChan; ++i) {
2790 //--cubic extrapolate--
2791 //int m = 4*(nPiece-1);
2792 //r1[i] = params[m] + params[m+1]*x1[i] + params[m+2]*x2[i] + params[m+3]*x3[i];
2793 //--linear extrapolate--
2794 //r1[i] = (r1[n-1] - r1[n])/(x1[n-1] - x1[n])*(x1[i] - x1[n]) + r1[n];
2795 //--constant--
2796 r1[i] = r1[n];
2797 }
2798 }
2799
2800 for (int i = 0; i < nChan; ++i) {
2801 residual[i] = z1[i] - r1[i];
2802 }
2803
2804 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
2805 break;
2806 } else {
2807 double stdDev = 0.0;
2808 for (int i = 0; i < nChan; ++i) {
2809 stdDev += residual[i]*residual[i]*(double)maskArray[i];
2810 }
2811 stdDev = sqrt(stdDev/(double)nData);
2812
2813 double thres = stdDev * thresClip;
2814 int newNData = 0;
2815 for (int i = 0; i < nChan; ++i) {
2816 if (abs(residual[i]) >= thres) {
2817 maskArray[i] = 0;
2818 }
2819 if (maskArray[i] > 0) {
2820 newNData++;
2821 }
2822 }
2823 if (newNData == nData) {
2824 break; //no more flag to add. iteration stops.
2825 } else {
2826 nData = newNData;
2827 }
2828 }
2829 }
2830
2831 nClipped = initNData - nData;
2832
2833 std::vector<float> result(nChan);
2834 if (getResidual) {
2835 for (int i = 0; i < nChan; ++i) {
2836 result[i] = (float)residual[i];
2837 }
2838 } else {
2839 for (int i = 0; i < nChan; ++i) {
2840 result[i] = (float)r1[i];
2841 }
2842 }
2843
2844 return result;
2845}
2846
2847void 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)
2848{
2849 nWaves.clear();
2850
2851 if (applyFFT) {
2852 string fftThAttr;
2853 float fftThSigma;
2854 int fftThTop;
2855 parseThresholdExpression(fftThresh, fftThAttr, fftThSigma, fftThTop);
2856 doSelectWaveNumbers(whichrow, chanMask, fftMethod, fftThSigma, fftThTop, fftThAttr, nWaves);
2857 }
2858
2859 addAuxWaveNumbers(whichrow, addNWaves, rejectNWaves, nWaves);
2860}
2861
2862void Scantable::parseThresholdExpression(const std::string& fftThresh, std::string& fftThAttr, float& fftThSigma, int& fftThTop)
2863{
2864 uInt idxSigma = fftThresh.find("sigma");
2865 uInt idxTop = fftThresh.find("top");
2866
2867 if (idxSigma == fftThresh.size() - 5) {
2868 std::istringstream is(fftThresh.substr(0, fftThresh.size() - 5));
2869 is >> fftThSigma;
2870 fftThAttr = "sigma";
2871 } else if (idxTop == 0) {
2872 std::istringstream is(fftThresh.substr(3));
2873 is >> fftThTop;
2874 fftThAttr = "top";
2875 } else {
2876 bool isNumber = true;
2877 for (uInt i = 0; i < fftThresh.size()-1; ++i) {
2878 char ch = (fftThresh.substr(i, 1).c_str())[0];
2879 if (!(isdigit(ch) || (fftThresh.substr(i, 1) == "."))) {
2880 isNumber = false;
2881 break;
2882 }
2883 }
2884 if (isNumber) {
2885 std::istringstream is(fftThresh);
2886 is >> fftThSigma;
2887 fftThAttr = "sigma";
2888 } else {
2889 throw(AipsError("fftthresh has a wrong value"));
2890 }
2891 }
2892}
2893
2894void 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)
2895{
2896 std::vector<float> fspec;
2897 if (fftMethod == "fft") {
2898 fspec = execFFT(whichrow, chanMask, false, true);
2899 //} else if (fftMethod == "lsp") {
2900 // fspec = lombScarglePeriodogram(whichrow);
2901 }
2902
2903 if (fftThAttr == "sigma") {
2904 float mean = 0.0;
2905 float mean2 = 0.0;
2906 for (uInt i = 0; i < fspec.size(); ++i) {
2907 mean += fspec[i];
2908 mean2 += fspec[i]*fspec[i];
2909 }
2910 mean /= float(fspec.size());
2911 mean2 /= float(fspec.size());
2912 float thres = mean + fftThSigma * float(sqrt(mean2 - mean*mean));
2913
2914 for (uInt i = 0; i < fspec.size(); ++i) {
2915 if (fspec[i] >= thres) {
2916 nWaves.push_back(i);
2917 }
2918 }
2919
2920 } else if (fftThAttr == "top") {
2921 for (int i = 0; i < fftThTop; ++i) {
2922 float max = 0.0;
2923 int maxIdx = 0;
2924 for (uInt j = 0; j < fspec.size(); ++j) {
2925 if (fspec[j] > max) {
2926 max = fspec[j];
2927 maxIdx = j;
2928 }
2929 }
2930 nWaves.push_back(maxIdx);
2931 fspec[maxIdx] = 0.0;
2932 }
2933
2934 }
2935
2936 if (nWaves.size() > 1) {
2937 sort(nWaves.begin(), nWaves.end());
2938 }
2939}
2940
2941void Scantable::addAuxWaveNumbers(const int whichrow, const std::vector<int>& addNWaves, const std::vector<int>& rejectNWaves, std::vector<int>& nWaves)
2942{
2943 std::vector<int> tempAddNWaves, tempRejectNWaves;
2944 for (uInt i = 0; i < addNWaves.size(); ++i) {
2945 tempAddNWaves.push_back(addNWaves[i]);
2946 }
2947 if ((tempAddNWaves.size() == 2) && (tempAddNWaves[1] == -999)) {
2948 setWaveNumberListUptoNyquistFreq(whichrow, tempAddNWaves);
2949 }
2950
2951 for (uInt i = 0; i < rejectNWaves.size(); ++i) {
2952 tempRejectNWaves.push_back(rejectNWaves[i]);
2953 }
2954 if ((tempRejectNWaves.size() == 2) && (tempRejectNWaves[1] == -999)) {
2955 setWaveNumberListUptoNyquistFreq(whichrow, tempRejectNWaves);
2956 }
2957
2958 for (uInt i = 0; i < tempAddNWaves.size(); ++i) {
2959 bool found = false;
2960 for (uInt j = 0; j < nWaves.size(); ++j) {
2961 if (nWaves[j] == tempAddNWaves[i]) {
2962 found = true;
2963 break;
2964 }
2965 }
2966 if (!found) nWaves.push_back(tempAddNWaves[i]);
2967 }
2968
2969 for (uInt i = 0; i < tempRejectNWaves.size(); ++i) {
2970 for (std::vector<int>::iterator j = nWaves.begin(); j != nWaves.end(); ) {
2971 if (*j == tempRejectNWaves[i]) {
2972 j = nWaves.erase(j);
2973 } else {
2974 ++j;
2975 }
2976 }
2977 }
2978
2979 if (nWaves.size() > 1) {
2980 sort(nWaves.begin(), nWaves.end());
2981 unique(nWaves.begin(), nWaves.end());
2982 }
2983}
2984
2985void Scantable::setWaveNumberListUptoNyquistFreq(const int whichrow, std::vector<int>& nWaves)
2986{
2987 if ((nWaves.size() == 2)&&(nWaves[1] == -999)) {
2988 int val = nWaves[0];
2989 int nyquistFreq = nchan(getIF(whichrow))/2+1;
2990 nWaves.clear();
2991 if (val > nyquistFreq) { // for safety, at least nWaves contains a constant; CAS-3759
2992 nWaves.push_back(0);
2993 }
2994 while (val <= nyquistFreq) {
2995 nWaves.push_back(val);
2996 val++;
2997 }
2998 }
2999}
3000
3001void 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)
3002{
3003 try {
3004 ofstream ofs;
3005 String coordInfo = "";
3006 bool hasSameNchan = true;
3007 bool outTextFile = false;
3008
3009 if (blfile != "") {
3010 ofs.open(blfile.c_str(), ios::out | ios::app);
3011 if (ofs) outTextFile = true;
3012 }
3013
3014 if (outLogger || outTextFile) {
3015 coordInfo = getCoordInfo()[0];
3016 if (coordInfo == "") coordInfo = "channel";
3017 hasSameNchan = hasSameNchanOverIFs();
3018 }
3019
3020 //Fitter fitter = Fitter();
3021 //fitter.setExpression("sinusoid", nWaves);
3022 //fitter.setIterClipping(thresClip, nIterClip);
3023
3024 int nRow = nrow();
3025 std::vector<bool> chanMask;
3026 std::vector<int> nWaves;
3027
3028 bool showProgress;
3029 int minNRow;
3030 parseProgressInfo(progressInfo, showProgress, minNRow);
3031
3032 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3033 chanMask = getCompositeChanMask(whichrow, mask);
3034 selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
3035
3036 //FOR DEBUGGING------------
3037 /*
3038 if (whichrow < 0) {// == nRow -1) {
3039 cout << "+++ i=" << setw(3) << whichrow << ", IF=" << setw(2) << getIF(whichrow);
3040 if (applyFFT) {
3041 cout << "[ ";
3042 for (uInt j = 0; j < nWaves.size(); ++j) {
3043 cout << nWaves[j] << ", ";
3044 }
3045 cout << " ] " << endl;
3046 }
3047 cout << flush;
3048 }
3049 */
3050 //-------------------------
3051
3052 //fitBaseline(chanMask, whichrow, fitter);
3053 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
3054 std::vector<float> params;
3055 int nClipped = 0;
3056 std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
3057 setSpectrum(res, whichrow);
3058 //
3059
3060 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "sinusoidBaseline()", params, nClipped);
3061 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3062 }
3063
3064 if (outTextFile) ofs.close();
3065
3066 } catch (...) {
3067 throw;
3068 }
3069}
3070
3071void 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)
3072{
3073 try {
3074 ofstream ofs;
3075 String coordInfo = "";
3076 bool hasSameNchan = true;
3077 bool outTextFile = false;
3078
3079 if (blfile != "") {
3080 ofs.open(blfile.c_str(), ios::out | ios::app);
3081 if (ofs) outTextFile = true;
3082 }
3083
3084 if (outLogger || outTextFile) {
3085 coordInfo = getCoordInfo()[0];
3086 if (coordInfo == "") coordInfo = "channel";
3087 hasSameNchan = hasSameNchanOverIFs();
3088 }
3089
3090 //Fitter fitter = Fitter();
3091 //fitter.setExpression("sinusoid", nWaves);
3092 //fitter.setIterClipping(thresClip, nIterClip);
3093
3094 int nRow = nrow();
3095 std::vector<bool> chanMask;
3096 std::vector<int> nWaves;
3097
3098 int minEdgeSize = getIFNos().size()*2;
3099 STLineFinder lineFinder = STLineFinder();
3100 lineFinder.setOptions(threshold, 3, chanAvgLimit);
3101
3102 bool showProgress;
3103 int minNRow;
3104 parseProgressInfo(progressInfo, showProgress, minNRow);
3105
3106 for (int whichrow = 0; whichrow < nRow; ++whichrow) {
3107
3108 //-------------------------------------------------------
3109 //chanMask = getCompositeChanMask(whichrow, mask, edge, minEdgeSize, lineFinder);
3110 //-------------------------------------------------------
3111 int edgeSize = edge.size();
3112 std::vector<int> currentEdge;
3113 if (edgeSize >= 2) {
3114 int idx = 0;
3115 if (edgeSize > 2) {
3116 if (edgeSize < minEdgeSize) {
3117 throw(AipsError("Length of edge element info is less than that of IFs"));
3118 }
3119 idx = 2 * getIF(whichrow);
3120 }
3121 currentEdge.push_back(edge[idx]);
3122 currentEdge.push_back(edge[idx+1]);
3123 } else {
3124 throw(AipsError("Wrong length of edge element"));
3125 }
3126 lineFinder.setData(getSpectrum(whichrow));
3127 lineFinder.findLines(getCompositeChanMask(whichrow, mask), currentEdge, whichrow);
3128 chanMask = lineFinder.getMask();
3129 //-------------------------------------------------------
3130
3131 selectWaveNumbers(whichrow, chanMask, applyFFT, fftMethod, fftThresh, addNWaves, rejectNWaves, nWaves);
3132
3133 //fitBaseline(chanMask, whichrow, fitter);
3134 //setSpectrum((getResidual ? fitter.getResidual() : fitter.getFit()), whichrow);
3135 std::vector<float> params;
3136 int nClipped = 0;
3137 std::vector<float> res = doSinusoidFitting(getSpectrum(whichrow), chanMask, nWaves, params, nClipped, thresClip, nIterClip, getResidual);
3138 setSpectrum(res, whichrow);
3139 //
3140
3141 outputFittingResult(outLogger, outTextFile, chanMask, whichrow, coordInfo, hasSameNchan, ofs, "autoSinusoidBaseline()", params, nClipped);
3142 showProgressOnTerminal(whichrow, nRow, showProgress, minNRow);
3143 }
3144
3145 if (outTextFile) ofs.close();
3146
3147 } catch (...) {
3148 throw;
3149 }
3150}
3151
3152std::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)
3153{
3154 if (data.size() != mask.size()) {
3155 throw(AipsError("data and mask sizes are not identical"));
3156 }
3157 if (data.size() < 2) {
3158 throw(AipsError("data size is too short"));
3159 }
3160 if (waveNumbers.size() == 0) {
3161 throw(AipsError("no wave numbers given"));
3162 }
3163 std::vector<int> nWaves; // sorted and uniqued array of wave numbers
3164 nWaves.reserve(waveNumbers.size());
3165 copy(waveNumbers.begin(), waveNumbers.end(), back_inserter(nWaves));
3166 sort(nWaves.begin(), nWaves.end());
3167 std::vector<int>::iterator end_it = unique(nWaves.begin(), nWaves.end());
3168 nWaves.erase(end_it, nWaves.end());
3169
3170 int minNWaves = nWaves[0];
3171 if (minNWaves < 0) {
3172 throw(AipsError("wave number must be positive or zero (i.e. constant)"));
3173 }
3174 bool hasConstantTerm = (minNWaves == 0);
3175
3176 int nChan = data.size();
3177 std::vector<int> maskArray;
3178 std::vector<int> x;
3179 for (int i = 0; i < nChan; ++i) {
3180 maskArray.push_back(mask[i] ? 1 : 0);
3181 if (mask[i]) {
3182 x.push_back(i);
3183 }
3184 }
3185
3186 int initNData = x.size();
3187
3188 int nData = initNData;
3189 int nDOF = nWaves.size() * 2 - (hasConstantTerm ? 1 : 0); //number of parameters to solve.
3190
3191 const double PI = 6.0 * asin(0.5); // PI (= 3.141592653...)
3192 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)
3193
3194 // xArray : contains elemental values for computing the least-square matrix.
3195 // xArray.size() is nDOF and xArray[*].size() is nChan.
3196 // Each xArray element are as follows:
3197 // xArray[0] = {1.0, 1.0, 1.0, ..., 1.0},
3198 // xArray[2n-1] = {sin(nPI/L*x[0]), sin(nPI/L*x[1]), ..., sin(nPI/L*x[nChan])},
3199 // xArray[2n] = {cos(nPI/L*x[0]), cos(nPI/L*x[1]), ..., cos(nPI/L*x[nChan])},
3200 // where (1 <= n <= nMaxWavesInSW),
3201 // or,
3202 // xArray[2n-1] = {sin(wn[n]PI/L*x[0]), sin(wn[n]PI/L*x[1]), ..., sin(wn[n]PI/L*x[nChan])},
3203 // xArray[2n] = {cos(wn[n]PI/L*x[0]), cos(wn[n]PI/L*x[1]), ..., cos(wn[n]PI/L*x[nChan])},
3204 // where wn[n] denotes waveNumbers[n] (1 <= n <= waveNumbers.size()).
3205 std::vector<std::vector<double> > xArray;
3206 if (hasConstantTerm) {
3207 std::vector<double> xu;
3208 for (int j = 0; j < nChan; ++j) {
3209 xu.push_back(1.0);
3210 }
3211 xArray.push_back(xu);
3212 }
3213 for (uInt i = (hasConstantTerm ? 1 : 0); i < nWaves.size(); ++i) {
3214 double xFactor = baseXFactor*(double)nWaves[i];
3215 std::vector<double> xs, xc;
3216 xs.clear();
3217 xc.clear();
3218 for (int j = 0; j < nChan; ++j) {
3219 xs.push_back(sin(xFactor*(double)j));
3220 xc.push_back(cos(xFactor*(double)j));
3221 }
3222 xArray.push_back(xs);
3223 xArray.push_back(xc);
3224 }
3225
3226 std::vector<double> z1, r1, residual;
3227 for (int i = 0; i < nChan; ++i) {
3228 z1.push_back((double)data[i]);
3229 r1.push_back(0.0);
3230 residual.push_back(0.0);
3231 }
3232
3233 for (int nClip = 0; nClip < nIterClip+1; ++nClip) {
3234 // xMatrix : horizontal concatenation of
3235 // the least-sq. matrix (left) and an
3236 // identity matrix (right).
3237 // the right part is used to calculate the inverse matrix of the left part.
3238 double xMatrix[nDOF][2*nDOF];
3239 double zMatrix[nDOF];
3240 for (int i = 0; i < nDOF; ++i) {
3241 for (int j = 0; j < 2*nDOF; ++j) {
3242 xMatrix[i][j] = 0.0;
3243 }
3244 xMatrix[i][nDOF+i] = 1.0;
3245 zMatrix[i] = 0.0;
3246 }
3247
3248 int nUseData = 0;
3249 for (int k = 0; k < nChan; ++k) {
3250 if (maskArray[k] == 0) continue;
3251
3252 for (int i = 0; i < nDOF; ++i) {
3253 for (int j = i; j < nDOF; ++j) {
3254 xMatrix[i][j] += xArray[i][k] * xArray[j][k];
3255 }
3256 zMatrix[i] += z1[k] * xArray[i][k];
3257 }
3258
3259 nUseData++;
3260 }
3261
3262 if (nUseData < 1) {
3263 throw(AipsError("all channels clipped or masked. can't execute fitting anymore."));
3264 }
3265
3266 for (int i = 0; i < nDOF; ++i) {
3267 for (int j = 0; j < i; ++j) {
3268 xMatrix[i][j] = xMatrix[j][i];
3269 }
3270 }
3271
3272 std::vector<double> invDiag;
3273 for (int i = 0; i < nDOF; ++i) {
3274 invDiag.push_back(1.0/xMatrix[i][i]);
3275 for (int j = 0; j < nDOF; ++j) {
3276 xMatrix[i][j] *= invDiag[i];
3277 }
3278 }
3279
3280 for (int k = 0; k < nDOF; ++k) {
3281 for (int i = 0; i < nDOF; ++i) {
3282 if (i != k) {
3283 double factor1 = xMatrix[k][k];
3284 double factor2 = xMatrix[i][k];
3285 for (int j = k; j < 2*nDOF; ++j) {
3286 xMatrix[i][j] *= factor1;
3287 xMatrix[i][j] -= xMatrix[k][j]*factor2;
3288 xMatrix[i][j] /= factor1;
3289 }
3290 }
3291 }
3292 double xDiag = xMatrix[k][k];
3293 for (int j = k; j < 2*nDOF; ++j) {
3294 xMatrix[k][j] /= xDiag;
3295 }
3296 }
3297
3298 for (int i = 0; i < nDOF; ++i) {
3299 for (int j = 0; j < nDOF; ++j) {
3300 xMatrix[i][nDOF+j] *= invDiag[j];
3301 }
3302 }
3303 //compute a vector y which consists of the coefficients of the sinusoids forming the
3304 //best-fit curves (a0,s1,c1,s2,c2,...), where a0 is constant and s* and c* are of sine
3305 //and cosine functions, respectively.
3306 std::vector<double> y;
3307 params.clear();
3308 for (int i = 0; i < nDOF; ++i) {
3309 y.push_back(0.0);
3310 for (int j = 0; j < nDOF; ++j) {
3311 y[i] += xMatrix[i][nDOF+j]*zMatrix[j];
3312 }
3313 params.push_back(y[i]);
3314 }
3315
3316 for (int i = 0; i < nChan; ++i) {
3317 r1[i] = y[0];
3318 for (int j = 1; j < nDOF; ++j) {
3319 r1[i] += y[j]*xArray[j][i];
3320 }
3321 residual[i] = z1[i] - r1[i];
3322 }
3323
3324 if ((nClip == nIterClip) || (thresClip <= 0.0)) {
3325 break;
3326 } else {
3327 double stdDev = 0.0;
3328 for (int i = 0; i < nChan; ++i) {
3329 stdDev += residual[i]*residual[i]*(double)maskArray[i];
3330 }
3331 stdDev = sqrt(stdDev/(double)nData);
3332
3333 double thres = stdDev * thresClip;
3334 int newNData = 0;
3335 for (int i = 0; i < nChan; ++i) {
3336 if (abs(residual[i]) >= thres) {
3337 maskArray[i] = 0;
3338 }
3339 if (maskArray[i] > 0) {
3340 newNData++;
3341 }
3342 }
3343 if (newNData == nData) {
3344 break; //no more flag to add. iteration stops.
3345 } else {
3346 nData = newNData;
3347 }
3348 }
3349 }
3350
3351 nClipped = initNData - nData;
3352
3353 std::vector<float> result;
3354 if (getResidual) {
3355 for (int i = 0; i < nChan; ++i) {
3356 result.push_back((float)residual[i]);
3357 }
3358 } else {
3359 for (int i = 0; i < nChan; ++i) {
3360 result.push_back((float)r1[i]);
3361 }
3362 }
3363
3364 return result;
3365}
3366
3367void Scantable::fitBaseline(const std::vector<bool>& mask, int whichrow, Fitter& fitter)
3368{
3369 std::vector<double> dAbcissa = getAbcissa(whichrow);
3370 std::vector<float> abcissa;
3371 for (uInt i = 0; i < dAbcissa.size(); ++i) {
3372 abcissa.push_back((float)dAbcissa[i]);
3373 }
3374 std::vector<float> spec = getSpectrum(whichrow);
3375
3376 fitter.setData(abcissa, spec, mask);
3377 fitter.lfit();
3378}
3379
3380std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask)
3381{
3382 std::vector<bool> mask = getMask(whichrow);
3383 uInt maskSize = mask.size();
3384 if (inMask.size() != 0) {
3385 if (maskSize != inMask.size()) {
3386 throw(AipsError("mask sizes are not the same."));
3387 }
3388 for (uInt i = 0; i < maskSize; ++i) {
3389 mask[i] = mask[i] && inMask[i];
3390 }
3391 }
3392
3393 return mask;
3394}
3395
3396/*
3397std::vector<bool> Scantable::getCompositeChanMask(int whichrow, const std::vector<bool>& inMask, const std::vector<int>& edge, const int minEdgeSize, STLineFinder& lineFinder)
3398{
3399 int edgeSize = edge.size();
3400 std::vector<int> currentEdge;
3401 if (edgeSize >= 2) {
3402 int idx = 0;
3403 if (edgeSize > 2) {
3404 if (edgeSize < minEdgeSize) {
3405 throw(AipsError("Length of edge element info is less than that of IFs"));
3406 }
3407 idx = 2 * getIF(whichrow);
3408 }
3409 currentEdge.push_back(edge[idx]);
3410 currentEdge.push_back(edge[idx+1]);
3411 } else {
3412 throw(AipsError("Wrong length of edge element"));
3413 }
3414
3415 lineFinder.setData(getSpectrum(whichrow));
3416 lineFinder.findLines(getCompositeChanMask(whichrow, inMask), currentEdge, whichrow);
3417
3418 return lineFinder.getMask();
3419}
3420*/
3421
3422/* for poly. the variations of outputFittingResult() should be merged into one eventually (2011/3/10 WK) */
3423void 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)
3424{
3425 if (outLogger || outTextFile) {
3426 std::vector<float> params = fitter.getParameters();
3427 std::vector<bool> fixed = fitter.getFixedParameters();
3428 float rms = getRms(chanMask, whichrow);
3429 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3430
3431 if (outLogger) {
3432 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3433 ols << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, false) << LogIO::POST ;
3434 }
3435 if (outTextFile) {
3436 ofs << formatBaselineParams(params, fixed, rms, -1, masklist, whichrow, true) << flush;
3437 }
3438 }
3439}
3440
3441/* for cspline. will be merged once cspline is available in fitter (2011/3/10 WK) */
3442void 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)
3443{
3444 if (outLogger || outTextFile) {
3445 float rms = getRms(chanMask, whichrow);
3446 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3447 std::vector<bool> fixed;
3448 fixed.clear();
3449
3450 if (outLogger) {
3451 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3452 ols << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3453 }
3454 if (outTextFile) {
3455 ofs << formatPiecewiseBaselineParams(edge, params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3456 }
3457 }
3458}
3459
3460/* for sinusoid. will be merged once sinusoid is available in fitter (2011/3/10 WK) */
3461void 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)
3462{
3463 if (outLogger || outTextFile) {
3464 float rms = getRms(chanMask, whichrow);
3465 String masklist = getMaskRangeList(chanMask, whichrow, coordInfo, hasSameNchan);
3466 std::vector<bool> fixed;
3467 fixed.clear();
3468
3469 if (outLogger) {
3470 LogIO ols(LogOrigin("Scantable", funcName, WHERE));
3471 ols << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, false) << LogIO::POST ;
3472 }
3473 if (outTextFile) {
3474 ofs << formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, true) << flush;
3475 }
3476 }
3477}
3478
3479void Scantable::parseProgressInfo(const std::string& progressInfo, bool& showProgress, int& minNRow)
3480{
3481 int idxDelimiter = progressInfo.find(",");
3482 if (idxDelimiter < 0) {
3483 throw(AipsError("wrong value in 'showprogress' parameter")) ;
3484 }
3485 showProgress = (progressInfo.substr(0, idxDelimiter) == "true");
3486 std::istringstream is(progressInfo.substr(idxDelimiter+1));
3487 is >> minNRow;
3488}
3489
3490void Scantable::showProgressOnTerminal(const int nProcessed, const int nTotal, const bool showProgress, const int nTotalThreshold)
3491{
3492 if (showProgress && (nTotal >= nTotalThreshold)) {
3493 int nInterval = int(floor(double(nTotal)/100.0));
3494 if (nInterval == 0) nInterval++;
3495
3496 if (nProcessed % nInterval == 0) {
3497 printf("\r"); //go to the head of line
3498 printf("\x1b[31m\x1b[1m"); //set red color, highlighted
3499 printf("[%3d%%]", (int)(100.0*(double(nProcessed+1))/(double(nTotal))) );
3500 printf("\x1b[39m\x1b[0m"); //set default attributes
3501 fflush(NULL);
3502 }
3503
3504 if (nProcessed == nTotal - 1) {
3505 printf("\r\x1b[K"); //clear
3506 fflush(NULL);
3507 }
3508
3509 }
3510}
3511
3512std::vector<float> Scantable::execFFT(const int whichrow, const std::vector<bool>& inMask, bool getRealImag, bool getAmplitudeOnly)
3513{
3514 std::vector<bool> mask = getMask(whichrow);
3515
3516 if (inMask.size() > 0) {
3517 uInt maskSize = mask.size();
3518 if (maskSize != inMask.size()) {
3519 throw(AipsError("mask sizes are not the same."));
3520 }
3521 for (uInt i = 0; i < maskSize; ++i) {
3522 mask[i] = mask[i] && inMask[i];
3523 }
3524 }
3525
3526 Vector<Float> spec = getSpectrum(whichrow);
3527 mathutil::doZeroOrderInterpolation(spec, mask);
3528
3529 FFTServer<Float,Complex> ffts;
3530 Vector<Complex> fftres;
3531 ffts.fft0(fftres, spec);
3532
3533 std::vector<float> res;
3534 float norm = float(2.0/double(spec.size()));
3535
3536 if (getRealImag) {
3537 for (uInt i = 0; i < fftres.size(); ++i) {
3538 res.push_back(real(fftres[i])*norm);
3539 res.push_back(imag(fftres[i])*norm);
3540 }
3541 } else {
3542 for (uInt i = 0; i < fftres.size(); ++i) {
3543 res.push_back(abs(fftres[i])*norm);
3544 if (!getAmplitudeOnly) res.push_back(arg(fftres[i]));
3545 }
3546 }
3547
3548 return res;
3549}
3550
3551
3552float Scantable::getRms(const std::vector<bool>& mask, int whichrow)
3553{
3554 Vector<Float> spec;
3555 specCol_.get(whichrow, spec);
3556
3557 float mean = 0.0;
3558 float smean = 0.0;
3559 int n = 0;
3560 for (uInt i = 0; i < spec.nelements(); ++i) {
3561 if (mask[i]) {
3562 mean += spec[i];
3563 smean += spec[i]*spec[i];
3564 n++;
3565 }
3566 }
3567
3568 mean /= (float)n;
3569 smean /= (float)n;
3570
3571 return sqrt(smean - mean*mean);
3572}
3573
3574
3575std::string Scantable::formatBaselineParamsHeader(int whichrow, const std::string& masklist, bool verbose) const
3576{
3577 ostringstream oss;
3578
3579 if (verbose) {
3580 oss << " Scan[" << getScan(whichrow) << "]";
3581 oss << " Beam[" << getBeam(whichrow) << "]";
3582 oss << " IF[" << getIF(whichrow) << "]";
3583 oss << " Pol[" << getPol(whichrow) << "]";
3584 oss << " Cycle[" << getCycle(whichrow) << "]: " << endl;
3585 oss << "Fitter range = " << masklist << endl;
3586 oss << "Baseline parameters" << endl;
3587 oss << flush;
3588 }
3589
3590 return String(oss);
3591}
3592
3593std::string Scantable::formatBaselineParamsFooter(float rms, int nClipped, bool verbose) const
3594{
3595 ostringstream oss;
3596
3597 if (verbose) {
3598 oss << "Results of baseline fit" << endl;
3599 oss << " rms = " << setprecision(6) << rms << endl;
3600 if (nClipped >= 0) {
3601 oss << " Number of clipped channels = " << nClipped << endl;
3602 }
3603 for (int i = 0; i < 60; ++i) {
3604 oss << "-";
3605 }
3606 oss << endl;
3607 oss << flush;
3608 }
3609
3610 return String(oss);
3611}
3612
3613std::string Scantable::formatBaselineParams(const std::vector<float>& params,
3614 const std::vector<bool>& fixed,
3615 float rms,
3616 int nClipped,
3617 const std::string& masklist,
3618 int whichrow,
3619 bool verbose,
3620 int start, int count,
3621 bool resetparamid) const
3622{
3623 int nParam = (int)(params.size());
3624
3625 if (nParam < 1) {
3626 return(" Not fitted");
3627 } else {
3628
3629 ostringstream oss;
3630 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3631
3632 if (start < 0) start = 0;
3633 if (count < 0) count = nParam;
3634 int end = start + count;
3635 if (end > nParam) end = nParam;
3636 int paramidoffset = (resetparamid) ? (-start) : 0;
3637
3638 for (int i = start; i < end; ++i) {
3639 if (i > start) {
3640 oss << ",";
3641 }
3642 std::string sFix = ((fixed.size() > 0) && (fixed[i]) && verbose) ? "(fixed)" : "";
3643 oss << " p" << (i+paramidoffset) << sFix << "= " << right << setw(13) << setprecision(6) << params[i];
3644 }
3645
3646 oss << endl;
3647 oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3648
3649 return String(oss);
3650 }
3651
3652}
3653
3654 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
3655{
3656 int nOutParam = (int)(params.size());
3657 int nPiece = (int)(ranges.size()) - 1;
3658
3659 if (nOutParam < 1) {
3660 return(" Not fitted");
3661 } else if (nPiece < 0) {
3662 return formatBaselineParams(params, fixed, rms, nClipped, masklist, whichrow, verbose);
3663 } else if (nPiece < 1) {
3664 return(" Bad count of the piece edge info");
3665 } else if (nOutParam % nPiece != 0) {
3666 return(" Bad count of the output baseline parameters");
3667 } else {
3668
3669 int nParam = nOutParam / nPiece;
3670
3671 ostringstream oss;
3672 oss << formatBaselineParamsHeader(whichrow, masklist, verbose);
3673
3674 stringstream ss;
3675 ss << ranges[nPiece] << flush;
3676 int wRange = ss.str().size() * 2 + 5;
3677
3678 for (int i = 0; i < nPiece; ++i) {
3679 ss.str("");
3680 ss << " [" << ranges[i] << "," << (ranges[i+1]-1) << "]";
3681 oss << left << setw(wRange) << ss.str();
3682 oss << formatBaselineParams(params, fixed, rms, 0, masklist, whichrow, false, i*nParam, nParam, true);
3683 }
3684
3685 oss << formatBaselineParamsFooter(rms, nClipped, verbose);
3686
3687 return String(oss);
3688 }
3689
3690}
3691
3692bool Scantable::hasSameNchanOverIFs()
3693{
3694 int nIF = nif(-1);
3695 int nCh;
3696 int totalPositiveNChan = 0;
3697 int nPositiveNChan = 0;
3698
3699 for (int i = 0; i < nIF; ++i) {
3700 nCh = nchan(i);
3701 if (nCh > 0) {
3702 totalPositiveNChan += nCh;
3703 nPositiveNChan++;
3704 }
3705 }
3706
3707 return (totalPositiveNChan == (nPositiveNChan * nchan(0)));
3708}
3709
3710std::string Scantable::getMaskRangeList(const std::vector<bool>& mask, int whichrow, const casa::String& coordInfo, bool hasSameNchan, bool verbose)
3711{
3712 if (mask.size() <= 0) {
3713 throw(AipsError("The mask elements should be > 0"));
3714 }
3715 int IF = getIF(whichrow);
3716 if (mask.size() != (uInt)nchan(IF)) {
3717 throw(AipsError("Number of channels in scantable != number of mask elements"));
3718 }
3719
3720 if (verbose) {
3721 LogIO logOs(LogOrigin("Scantable", "getMaskRangeList()", WHERE));
3722 logOs << LogIO::WARN << "The current mask window unit is " << coordInfo;
3723 if (!hasSameNchan) {
3724 logOs << endl << "This mask is only valid for IF=" << IF;
3725 }
3726 logOs << LogIO::POST;
3727 }
3728
3729 std::vector<double> abcissa = getAbcissa(whichrow);
3730 std::vector<int> edge = getMaskEdgeIndices(mask);
3731
3732 ostringstream oss;
3733 oss.setf(ios::fixed);
3734 oss << setprecision(1) << "[";
3735 for (uInt i = 0; i < edge.size(); i+=2) {
3736 if (i > 0) oss << ",";
3737 oss << "[" << (float)abcissa[edge[i]] << "," << (float)abcissa[edge[i+1]] << "]";
3738 }
3739 oss << "]" << flush;
3740
3741 return String(oss);
3742}
3743
3744std::vector<int> Scantable::getMaskEdgeIndices(const std::vector<bool>& mask)
3745{
3746 if (mask.size() <= 0) {
3747 throw(AipsError("The mask elements should be > 0"));
3748 }
3749
3750 std::vector<int> out, startIndices, endIndices;
3751 int maskSize = mask.size();
3752
3753 startIndices.clear();
3754 endIndices.clear();
3755
3756 if (mask[0]) {
3757 startIndices.push_back(0);
3758 }
3759 for (int i = 1; i < maskSize; ++i) {
3760 if ((!mask[i-1]) && mask[i]) {
3761 startIndices.push_back(i);
3762 } else if (mask[i-1] && (!mask[i])) {
3763 endIndices.push_back(i-1);
3764 }
3765 }
3766 if (mask[maskSize-1]) {
3767 endIndices.push_back(maskSize-1);
3768 }
3769
3770 if (startIndices.size() != endIndices.size()) {
3771 throw(AipsError("Inconsistent Mask Size: bad data?"));
3772 }
3773 for (uInt i = 0; i < startIndices.size(); ++i) {
3774 if (startIndices[i] > endIndices[i]) {
3775 throw(AipsError("Mask start index > mask end index"));
3776 }
3777 }
3778
3779 out.clear();
3780 for (uInt i = 0; i < startIndices.size(); ++i) {
3781 out.push_back(startIndices[i]);
3782 out.push_back(endIndices[i]);
3783 }
3784
3785 return out;
3786}
3787
3788vector<float> Scantable::getTsysSpectrum( int whichrow ) const
3789{
3790 Vector<Float> tsys( tsysCol_(whichrow) ) ;
3791 vector<float> stlTsys ;
3792 tsys.tovector( stlTsys ) ;
3793 return stlTsys ;
3794}
3795
3796
3797}
3798//namespace asap
Note: See TracBrowser for help on using the repository browser.