source: trunk/src/Scantable.cpp@ 2123

Last change on this file since 2123 was 2111, checked in by Kana Sugimoto, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs: sdlist unit tests

Put in Release Notes: No

Module(s): scantable

Description:

Put header summarizing part of code Scantable::summary to a new function Scantable::headerSummary.
The function is accessible from python level using scantable._list_header()


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