source: trunk/src/Scantable.cpp@ 2280

Last change on this file since 2280 was 2277, checked in by Malte Marquarding, 13 years ago

revert accidental changes

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