source: branches/parallel/src/Scantable.cpp@ 2041

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

New Development: No

JIRA Issue: Yes (CAS-2817/ASAP-232)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: added a new function SrcType::getName(int)

Test Programs: sdlist (CASA) or scantable.summary()

Put in Release Notes: No

Module(s): scantable, sdlist, sdplot

Description:

Added a new function SrcType::getName(int srctype) which returns
a source type in string.
The source type of each scan is retuned by Scantable::summary.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.8 KB
Line 
1//
2// C++ Implementation: Scantable
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <map>
13#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 "MathUtils.h"
70
71using namespace casa;
72
73namespace asap {
74
75std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
76
77void Scantable::initFactories() {
78 if ( factories_.empty() ) {
79 Scantable::factories_["linear"] = &STPolLinear::myFactory;
80 Scantable::factories_["circular"] = &STPolCircular::myFactory;
81 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
82 }
83}
84
85Scantable::Scantable(Table::TableType ttype) :
86 type_(ttype)
87{
88 initFactories();
89 setupMainTable();
90 freqTable_ = STFrequencies(*this);
91 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
92 weatherTable_ = STWeather(*this);
93 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
94 focusTable_ = STFocus(*this);
95 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
96 tcalTable_ = STTcal(*this);
97 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
98 moleculeTable_ = STMolecules(*this);
99 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
100 historyTable_ = STHistory(*this);
101 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
102 fitTable_ = STFit(*this);
103 table_.rwKeywordSet().defineTable("FIT", fitTable_.table());
104 table_.tableInfo().setType( "Scantable" ) ;
105 originalTable_ = table_;
106 attach();
107}
108
109Scantable::Scantable(const std::string& name, Table::TableType ttype) :
110 type_(ttype)
111{
112 initFactories();
113
114 Table tab(name, Table::Update);
115 uInt version = tab.keywordSet().asuInt("VERSION");
116 if (version != version_) {
117 throw(AipsError("Unsupported version of ASAP file."));
118 }
119 if ( type_ == Table::Memory ) {
120 table_ = tab.copyToMemoryTable(generateName());
121 } else {
122 table_ = tab;
123 }
124 table_.tableInfo().setType( "Scantable" ) ;
125
126 attachSubtables();
127 originalTable_ = table_;
128 attach();
129}
130/*
131Scantable::Scantable(const std::string& name, Table::TableType ttype) :
132 type_(ttype)
133{
134 initFactories();
135 Table tab(name, Table::Update);
136 uInt version = tab.keywordSet().asuInt("VERSION");
137 if (version != version_) {
138 throw(AipsError("Unsupported version of ASAP file."));
139 }
140 if ( type_ == Table::Memory ) {
141 table_ = tab.copyToMemoryTable(generateName());
142 } else {
143 table_ = tab;
144 }
145
146 attachSubtables();
147 originalTable_ = table_;
148 attach();
149}
150*/
151
152Scantable::Scantable( const Scantable& other, bool clear )
153{
154 // with or without data
155 String newname = String(generateName());
156 type_ = other.table_.tableType();
157 if ( other.table_.tableType() == Table::Memory ) {
158 if ( clear ) {
159 table_ = TableCopy::makeEmptyMemoryTable(newname,
160 other.table_, True);
161 } else
162 table_ = other.table_.copyToMemoryTable(newname);
163 } else {
164 other.table_.deepCopy(newname, Table::New, False,
165 other.table_.endianFormat(),
166 Bool(clear));
167 table_ = Table(newname, Table::Update);
168 table_.markForDelete();
169 }
170 table_.tableInfo().setType( "Scantable" ) ;
171 /// @todo reindex SCANNO, recompute nbeam, nif, npol
172 if ( clear ) copySubtables(other);
173 attachSubtables();
174 originalTable_ = table_;
175 attach();
176}
177
178void Scantable::copySubtables(const Scantable& other) {
179 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
180 TableCopy::copyRows(t, other.freqTable_.table());
181 t = table_.rwKeywordSet().asTable("FOCUS");
182 TableCopy::copyRows(t, other.focusTable_.table());
183 t = table_.rwKeywordSet().asTable("WEATHER");
184 TableCopy::copyRows(t, other.weatherTable_.table());
185 t = table_.rwKeywordSet().asTable("TCAL");
186 TableCopy::copyRows(t, other.tcalTable_.table());
187 t = table_.rwKeywordSet().asTable("MOLECULES");
188 TableCopy::copyRows(t, other.moleculeTable_.table());
189 t = table_.rwKeywordSet().asTable("HISTORY");
190 TableCopy::copyRows(t, other.historyTable_.table());
191 t = table_.rwKeywordSet().asTable("FIT");
192 TableCopy::copyRows(t, other.fitTable_.table());
193}
194
195void Scantable::attachSubtables()
196{
197 freqTable_ = STFrequencies(table_);
198 focusTable_ = STFocus(table_);
199 weatherTable_ = STWeather(table_);
200 tcalTable_ = STTcal(table_);
201 moleculeTable_ = STMolecules(table_);
202 historyTable_ = STHistory(table_);
203 fitTable_ = STFit(table_);
204}
205
206Scantable::~Scantable()
207{
208 //cout << "~Scantable() " << this << endl;
209}
210
211void Scantable::setupMainTable()
212{
213 TableDesc td("", "1", TableDesc::Scratch);
214 td.comment() = "An ASAP Scantable";
215 td.rwKeywordSet().define("VERSION", uInt(version_));
216
217 // n Cycles
218 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
219 // new index every nBeam x nIF x nPol
220 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
221
222 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
223 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
224 // linear, circular, stokes
225 td.rwKeywordSet().define("POLTYPE", String("linear"));
226 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
227
228 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
229 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
230
231 ScalarColumnDesc<Int> refbeamnoColumn("REFBEAMNO");
232 refbeamnoColumn.setDefault(Int(-1));
233 td.addColumn(refbeamnoColumn);
234
235 ScalarColumnDesc<uInt> flagrowColumn("FLAGROW");
236 flagrowColumn.setDefault(uInt(0));
237 td.addColumn(flagrowColumn);
238
239 td.addColumn(ScalarColumnDesc<Double>("TIME"));
240 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
241 TableMeasValueDesc measVal(td, "TIME");
242 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
243 mepochCol.write(td);
244
245 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
246
247 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
248 // Type of source (on=0, off=1, other=-1)
249 ScalarColumnDesc<Int> stypeColumn("SRCTYPE");
250 stypeColumn.setDefault(Int(-1));
251 td.addColumn(stypeColumn);
252 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
253
254 //The actual Data Vectors
255 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
256 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
257 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
258
259 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
260 IPosition(1,2),
261 ColumnDesc::Direct));
262 TableMeasRefDesc mdirRef(MDirection::J2000); // default
263 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
264 // the TableMeasDesc gives the column a type
265 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
266 // a uder set table type e.g. GALCTIC, B1950 ...
267 td.rwKeywordSet().define("DIRECTIONREF", String("J2000"));
268 // writing create the measure column
269 mdirCol.write(td);
270 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
271 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
272 td.addColumn(ScalarColumnDesc<Float>("OPACITY"));
273
274 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
275 ScalarColumnDesc<Int> fitColumn("FIT_ID");
276 fitColumn.setDefault(Int(-1));
277 td.addColumn(fitColumn);
278
279 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
280 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
281
282 // columns which just get dragged along, as they aren't used in asap
283 td.addColumn(ScalarColumnDesc<Double>("SRCVELOCITY"));
284 td.addColumn(ArrayColumnDesc<Double>("SRCPROPERMOTION"));
285 td.addColumn(ArrayColumnDesc<Double>("SRCDIRECTION"));
286 td.addColumn(ArrayColumnDesc<Double>("SCANRATE"));
287
288 td.rwKeywordSet().define("OBSMODE", String(""));
289
290 // Now create Table SetUp from the description.
291 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
292 table_ = Table(aNewTab, type_, 0);
293 originalTable_ = table_;
294}
295
296void Scantable::attach()
297{
298 timeCol_.attach(table_, "TIME");
299 srcnCol_.attach(table_, "SRCNAME");
300 srctCol_.attach(table_, "SRCTYPE");
301 specCol_.attach(table_, "SPECTRA");
302 flagsCol_.attach(table_, "FLAGTRA");
303 tsysCol_.attach(table_, "TSYS");
304 cycleCol_.attach(table_,"CYCLENO");
305 scanCol_.attach(table_, "SCANNO");
306 beamCol_.attach(table_, "BEAMNO");
307 ifCol_.attach(table_, "IFNO");
308 polCol_.attach(table_, "POLNO");
309 integrCol_.attach(table_, "INTERVAL");
310 azCol_.attach(table_, "AZIMUTH");
311 elCol_.attach(table_, "ELEVATION");
312 dirCol_.attach(table_, "DIRECTION");
313 fldnCol_.attach(table_, "FIELDNAME");
314 rbeamCol_.attach(table_, "REFBEAMNO");
315
316 mweatheridCol_.attach(table_,"WEATHER_ID");
317 mfitidCol_.attach(table_,"FIT_ID");
318 mfreqidCol_.attach(table_, "FREQ_ID");
319 mtcalidCol_.attach(table_, "TCAL_ID");
320 mfocusidCol_.attach(table_, "FOCUS_ID");
321 mmolidCol_.attach(table_, "MOLECULE_ID");
322
323 //Add auxiliary column for row-based flagging (CAS-1433 Wataru Kawasaki)
324 attachAuxColumnDef(flagrowCol_, "FLAGROW", 0);
325
326}
327
328template<class T, class T2>
329void Scantable::attachAuxColumnDef(ScalarColumn<T>& col,
330 const String& colName,
331 const T2& defValue)
332{
333 try {
334 col.attach(table_, colName);
335 } catch (TableError& err) {
336 String errMesg = err.getMesg();
337 if (errMesg == "Table column " + colName + " is unknown") {
338 table_.addColumn(ScalarColumnDesc<T>(colName));
339 col.attach(table_, colName);
340 col.fillColumn(static_cast<T>(defValue));
341 } else {
342 throw;
343 }
344 } catch (...) {
345 throw;
346 }
347}
348
349template<class T, class T2>
350void Scantable::attachAuxColumnDef(ArrayColumn<T>& col,
351 const String& colName,
352 const Array<T2>& defValue)
353{
354 try {
355 col.attach(table_, colName);
356 } catch (TableError& err) {
357 String errMesg = err.getMesg();
358 if (errMesg == "Table column " + colName + " is unknown") {
359 table_.addColumn(ArrayColumnDesc<T>(colName));
360 col.attach(table_, colName);
361
362 int size = 0;
363 ArrayIterator<T2>& it = defValue.begin();
364 while (it != defValue.end()) {
365 ++size;
366 ++it;
367 }
368 IPosition ip(1, size);
369 Array<T>& arr(ip);
370 for (int i = 0; i < size; ++i)
371 arr[i] = static_cast<T>(defValue[i]);
372
373 col.fillColumn(arr);
374 } else {
375 throw;
376 }
377 } catch (...) {
378 throw;
379 }
380}
381
382void Scantable::setHeader(const STHeader& sdh)
383{
384 table_.rwKeywordSet().define("nIF", sdh.nif);
385 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
386 table_.rwKeywordSet().define("nPol", sdh.npol);
387 table_.rwKeywordSet().define("nChan", sdh.nchan);
388 table_.rwKeywordSet().define("Observer", sdh.observer);
389 table_.rwKeywordSet().define("Project", sdh.project);
390 table_.rwKeywordSet().define("Obstype", sdh.obstype);
391 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
392 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
393 table_.rwKeywordSet().define("Equinox", sdh.equinox);
394 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
395 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
396 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
397 table_.rwKeywordSet().define("UTC", sdh.utc);
398 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
399 table_.rwKeywordSet().define("Epoch", sdh.epoch);
400 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
401}
402
403STHeader Scantable::getHeader() const
404{
405 STHeader sdh;
406 table_.keywordSet().get("nBeam",sdh.nbeam);
407 table_.keywordSet().get("nIF",sdh.nif);
408 table_.keywordSet().get("nPol",sdh.npol);
409 table_.keywordSet().get("nChan",sdh.nchan);
410 table_.keywordSet().get("Observer", sdh.observer);
411 table_.keywordSet().get("Project", sdh.project);
412 table_.keywordSet().get("Obstype", sdh.obstype);
413 table_.keywordSet().get("AntennaName", sdh.antennaname);
414 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
415 table_.keywordSet().get("Equinox", sdh.equinox);
416 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
417 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
418 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
419 table_.keywordSet().get("UTC", sdh.utc);
420 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
421 table_.keywordSet().get("Epoch", sdh.epoch);
422 table_.keywordSet().get("POLTYPE", sdh.poltype);
423 return sdh;
424}
425
426void Scantable::setSourceType( int stype )
427{
428 if ( stype < 0 || stype > 1 )
429 throw(AipsError("Illegal sourcetype."));
430 TableVector<Int> tabvec(table_, "SRCTYPE");
431 tabvec = Int(stype);
432}
433
434bool Scantable::conformant( const Scantable& other )
435{
436 return this->getHeader().conformant(other.getHeader());
437}
438
439
440
441std::string Scantable::formatSec(Double x) const
442{
443 Double xcop = x;
444 MVTime mvt(xcop/24./3600.); // make days
445
446 if (x < 59.95)
447 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
448 else if (x < 3599.95)
449 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
450 else {
451 ostringstream oss;
452 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
453 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
454 return String(oss);
455 }
456};
457
458std::string Scantable::formatDirection(const MDirection& md) const
459{
460 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
461 Int prec = 7;
462
463 MVAngle mvLon(t[0]);
464 String sLon = mvLon.string(MVAngle::TIME,prec);
465 uInt tp = md.getRef().getType();
466 if (tp == MDirection::GALACTIC ||
467 tp == MDirection::SUPERGAL ) {
468 sLon = mvLon(0.0).string(MVAngle::ANGLE_CLEAN,prec);
469 }
470 MVAngle mvLat(t[1]);
471 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
472 return sLon + String(" ") + sLat;
473}
474
475
476std::string Scantable::getFluxUnit() const
477{
478 return table_.keywordSet().asString("FluxUnit");
479}
480
481void Scantable::setFluxUnit(const std::string& unit)
482{
483 String tmp(unit);
484 Unit tU(tmp);
485 if (tU==Unit("K") || tU==Unit("Jy")) {
486 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
487 } else {
488 throw AipsError("Illegal unit - must be compatible with Jy or K");
489 }
490}
491
492void Scantable::setInstrument(const std::string& name)
493{
494 bool throwIt = true;
495 // create an Instrument to see if this is valid
496 STAttr::convertInstrument(name, throwIt);
497 String nameU(name);
498 nameU.upcase();
499 table_.rwKeywordSet().define(String("AntennaName"), nameU);
500}
501
502void Scantable::setFeedType(const std::string& feedtype)
503{
504 if ( Scantable::factories_.find(feedtype) == Scantable::factories_.end() ) {
505 std::string msg = "Illegal feed type "+ feedtype;
506 throw(casa::AipsError(msg));
507 }
508 table_.rwKeywordSet().define(String("POLTYPE"), feedtype);
509}
510
511MPosition Scantable::getAntennaPosition() const
512{
513 Vector<Double> antpos;
514 table_.keywordSet().get("AntennaPosition", antpos);
515 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
516 return MPosition(mvpos);
517}
518
519void Scantable::makePersistent(const std::string& filename)
520{
521 String inname(filename);
522 Path path(inname);
523 /// @todo reindex SCANNO, recompute nbeam, nif, npol
524 inname = path.expandedName();
525 // WORKAROUND !!! for Table bug
526 // Remove when fixed in casacore
527 if ( table_.tableType() == Table::Memory && !selector_.empty() ) {
528 Table tab = table_.copyToMemoryTable(generateName());
529 tab.deepCopy(inname, Table::New);
530 tab.markForDelete();
531
532 } else {
533 table_.deepCopy(inname, Table::New);
534 }
535}
536
537int Scantable::nbeam( int scanno ) const
538{
539 if ( scanno < 0 ) {
540 Int n;
541 table_.keywordSet().get("nBeam",n);
542 return int(n);
543 } else {
544 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
545 Table t = table_(table_.col("SCANNO") == scanno);
546 ROTableRow row(t);
547 const TableRecord& rec = row.get(0);
548 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
549 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
550 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
551 ROTableVector<uInt> v(subt, "BEAMNO");
552 return int(v.nelements());
553 }
554 return 0;
555}
556
557int Scantable::nif( int scanno ) const
558{
559 if ( scanno < 0 ) {
560 Int n;
561 table_.keywordSet().get("nIF",n);
562 return int(n);
563 } else {
564 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
565 Table t = table_(table_.col("SCANNO") == scanno);
566 ROTableRow row(t);
567 const TableRecord& rec = row.get(0);
568 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
569 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
570 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
571 if ( subt.nrow() == 0 ) return 0;
572 ROTableVector<uInt> v(subt, "IFNO");
573 return int(v.nelements());
574 }
575 return 0;
576}
577
578int Scantable::npol( int scanno ) const
579{
580 if ( scanno < 0 ) {
581 Int n;
582 table_.keywordSet().get("nPol",n);
583 return n;
584 } else {
585 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
586 Table t = table_(table_.col("SCANNO") == scanno);
587 ROTableRow row(t);
588 const TableRecord& rec = row.get(0);
589 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
590 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
591 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
592 if ( subt.nrow() == 0 ) return 0;
593 ROTableVector<uInt> v(subt, "POLNO");
594 return int(v.nelements());
595 }
596 return 0;
597}
598
599int Scantable::ncycle( int scanno ) const
600{
601 if ( scanno < 0 ) {
602 Block<String> cols(2);
603 cols[0] = "SCANNO";
604 cols[1] = "CYCLENO";
605 TableIterator it(table_, cols);
606 int n = 0;
607 while ( !it.pastEnd() ) {
608 ++n;
609 ++it;
610 }
611 return n;
612 } else {
613 Table t = table_(table_.col("SCANNO") == scanno);
614 ROTableRow row(t);
615 const TableRecord& rec = row.get(0);
616 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
617 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
618 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
619 if ( subt.nrow() == 0 ) return 0;
620 return int(subt.nrow());
621 }
622 return 0;
623}
624
625
626int Scantable::nrow( int scanno ) const
627{
628 return int(table_.nrow());
629}
630
631int Scantable::nchan( int ifno ) const
632{
633 if ( ifno < 0 ) {
634 Int n;
635 table_.keywordSet().get("nChan",n);
636 return int(n);
637 } else {
638 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't
639 // vary with these
640 Table t = table_(table_.col("IFNO") == ifno);
641 if ( t.nrow() == 0 ) return 0;
642 ROArrayColumn<Float> v(t, "SPECTRA");
643 return v.shape(0)(0);
644 }
645 return 0;
646}
647
648int Scantable::nscan() const {
649 Vector<uInt> scannos(scanCol_.getColumn());
650 uInt nout = genSort( scannos, Sort::Ascending,
651 Sort::QuickSort|Sort::NoDuplicates );
652 return int(nout);
653}
654
655int Scantable::getChannels(int whichrow) const
656{
657 return specCol_.shape(whichrow)(0);
658}
659
660int Scantable::getBeam(int whichrow) const
661{
662 return beamCol_(whichrow);
663}
664
665std::vector<uint> Scantable::getNumbers(const ScalarColumn<uInt>& col) const
666{
667 Vector<uInt> nos(col.getColumn());
668 uInt n = genSort( nos, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates );
669 nos.resize(n, True);
670 std::vector<uint> stlout;
671 nos.tovector(stlout);
672 return stlout;
673}
674
675int Scantable::getIF(int whichrow) const
676{
677 return ifCol_(whichrow);
678}
679
680int Scantable::getPol(int whichrow) const
681{
682 return polCol_(whichrow);
683}
684
685std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
686{
687 return formatTime(me, showdate, 0);
688}
689
690std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
691{
692 MVTime mvt(me.getValue());
693 if (showdate)
694 //mvt.setFormat(MVTime::YMD);
695 mvt.setFormat(MVTime::YMD, prec);
696 else
697 //mvt.setFormat(MVTime::TIME);
698 mvt.setFormat(MVTime::TIME, prec);
699 ostringstream oss;
700 oss << mvt;
701 return String(oss);
702}
703
704void Scantable::calculateAZEL()
705{
706 MPosition mp = getAntennaPosition();
707 MEpoch::ROScalarColumn timeCol(table_, "TIME");
708 ostringstream oss;
709 oss << "Computed azimuth/elevation using " << endl
710 << mp << endl;
711 for (Int i=0; i<nrow(); ++i) {
712 MEpoch me = timeCol(i);
713 MDirection md = getDirection(i);
714 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
715 << endl << " => ";
716 MeasFrame frame(mp, me);
717 Vector<Double> azel =
718 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
719 frame)
720 )().getAngle("rad").getValue();
721 azCol_.put(i,Float(azel[0]));
722 elCol_.put(i,Float(azel[1]));
723 oss << "azel: " << azel[0]/C::pi*180.0 << " "
724 << azel[1]/C::pi*180.0 << " (deg)" << endl;
725 }
726 pushLog(String(oss));
727}
728
729void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
730{
731 for (uInt i=0; i<table_.nrow(); ++i) {
732 Vector<uChar> flgs = flagsCol_(i);
733 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
734 flagsCol_.put(i, flgs);
735 }
736}
737
738std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
739{
740 Vector<uChar> flags;
741 flagsCol_.get(uInt(whichrow), flags);
742 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
743 Vector<Bool> bflag(flags.shape());
744 convertArray(bflag, flags);
745 //bflag = !bflag;
746
747 std::vector<bool> mask;
748 bflag.tovector(mask);
749 return mask;
750}
751
752void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
753 Vector<uChar> flgs)
754{
755 Vector<Float> spcs = specCol_(whichrow);
756 uInt nchannel = nchan();
757 if (spcs.nelements() != nchannel) {
758 throw(AipsError("Data has incorrect number of channels"));
759 }
760 uChar userflag = 1 << 7;
761 if (unflag) {
762 userflag = 0 << 7;
763 }
764 if (clipoutside) {
765 for (uInt j = 0; j < nchannel; ++j) {
766 Float spc = spcs(j);
767 if ((spc >= uthres) || (spc <= dthres)) {
768 flgs(j) = userflag;
769 }
770 }
771 } else {
772 for (uInt j = 0; j < nchannel; ++j) {
773 Float spc = spcs(j);
774 if ((spc < uthres) && (spc > dthres)) {
775 flgs(j) = userflag;
776 }
777 }
778 }
779}
780
781
782void Scantable::flag( int whichrow, const std::vector<bool>& msk, bool unflag ) {
783 std::vector<bool>::const_iterator it;
784 uInt ntrue = 0;
785 if (whichrow >= int(table_.nrow()) ) {
786 throw(AipsError("Invalid row number"));
787 }
788 for (it = msk.begin(); it != msk.end(); ++it) {
789 if ( *it ) {
790 ntrue++;
791 }
792 }
793 //if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
794 if ( whichrow == -1 && !unflag && selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
795 throw(AipsError("Trying to flag whole scantable."));
796 uChar userflag = 1 << 7;
797 if ( unflag ) {
798 userflag = 0 << 7;
799 }
800 if (whichrow > -1 ) {
801 applyChanFlag(uInt(whichrow), msk, userflag);
802 } else {
803 for ( uInt i=0; i<table_.nrow(); ++i) {
804 applyChanFlag(i, msk, userflag);
805 }
806 }
807}
808
809void Scantable::applyChanFlag( uInt whichrow, const std::vector<bool>& msk, uChar flagval )
810{
811 if (whichrow >= table_.nrow() ) {
812 throw( casa::indexError<int>( whichrow, "asap::Scantable::applyChanFlag: Invalid row number" ) );
813 }
814 Vector<uChar> flgs = flagsCol_(whichrow);
815 if ( msk.size() == 0 ) {
816 flgs = flagval;
817 flagsCol_.put(whichrow, flgs);
818 return;
819 }
820 if ( int(msk.size()) != nchan() ) {
821 throw(AipsError("Mask has incorrect number of channels."));
822 }
823 if ( flgs.nelements() != msk.size() ) {
824 throw(AipsError("Mask has incorrect number of channels."
825 " Probably varying with IF. Please flag per IF"));
826 }
827 std::vector<bool>::const_iterator it;
828 uInt j = 0;
829 for (it = msk.begin(); it != msk.end(); ++it) {
830 if ( *it ) {
831 flgs(j) = flagval;
832 }
833 ++j;
834 }
835 flagsCol_.put(whichrow, flgs);
836}
837
838void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
839{
840 if ( selector_.empty() && (rows.size() == table_.nrow()) )
841 throw(AipsError("Trying to flag whole scantable."));
842
843 uInt rowflag = (unflag ? 0 : 1);
844 std::vector<uInt>::const_iterator it;
845 for (it = rows.begin(); it != rows.end(); ++it)
846 flagrowCol_.put(*it, rowflag);
847}
848
849std::vector<bool> Scantable::getMask(int whichrow) const
850{
851 Vector<uChar> flags;
852 flagsCol_.get(uInt(whichrow), flags);
853 Vector<Bool> bflag(flags.shape());
854 convertArray(bflag, flags);
855 bflag = !bflag;
856 std::vector<bool> mask;
857 bflag.tovector(mask);
858 return mask;
859}
860
861std::vector<float> Scantable::getSpectrum( int whichrow,
862 const std::string& poltype ) const
863{
864 String ptype = poltype;
865 if (poltype == "" ) ptype = getPolType();
866 if ( whichrow < 0 || whichrow >= nrow() )
867 throw(AipsError("Illegal row number."));
868 std::vector<float> out;
869 Vector<Float> arr;
870 uInt requestedpol = polCol_(whichrow);
871 String basetype = getPolType();
872 if ( ptype == basetype ) {
873 specCol_.get(whichrow, arr);
874 } else {
875 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
876 basetype));
877 uInt row = uInt(whichrow);
878 stpol->setSpectra(getPolMatrix(row));
879 Float fang,fhand,parang;
880 fang = focusTable_.getTotalAngle(mfocusidCol_(row));
881 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
882 stpol->setPhaseCorrections(fang, fhand);
883 arr = stpol->getSpectrum(requestedpol, ptype);
884 }
885 if ( arr.nelements() == 0 )
886 pushLog("Not enough polarisations present to do the conversion.");
887 arr.tovector(out);
888 return out;
889}
890
891void Scantable::setSpectrum( const std::vector<float>& spec,
892 int whichrow )
893{
894 Vector<Float> spectrum(spec);
895 Vector<Float> arr;
896 specCol_.get(whichrow, arr);
897 if ( spectrum.nelements() != arr.nelements() )
898 throw AipsError("The spectrum has incorrect number of channels.");
899 specCol_.put(whichrow, spectrum);
900}
901
902
903String Scantable::generateName()
904{
905 return (File::newUniqueName("./","temp")).baseName();
906}
907
908const casa::Table& Scantable::table( ) const
909{
910 return table_;
911}
912
913casa::Table& Scantable::table( )
914{
915 return table_;
916}
917
918std::string Scantable::getPolType() const
919{
920 return table_.keywordSet().asString("POLTYPE");
921}
922
923void Scantable::unsetSelection()
924{
925 table_ = originalTable_;
926 attach();
927 selector_.reset();
928}
929
930void Scantable::setSelection( const STSelector& selection )
931{
932 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
933 if ( tab.nrow() == 0 ) {
934 throw(AipsError("Selection contains no data. Not applying it."));
935 }
936 table_ = tab;
937 attach();
938 selector_ = selection;
939}
940
941std::string Scantable::summary( bool verbose )
942{
943 // Format header info
944 ostringstream oss;
945 oss << endl;
946 oss << asap::SEPERATOR << endl;
947 oss << " Scan Table Summary" << endl;
948 oss << asap::SEPERATOR << endl;
949 oss.flags(std::ios_base::left);
950 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
951 << setw(15) << "IFs:" << setw(4) << nif() << endl
952 << setw(15) << "Polarisations:" << setw(4) << npol()
953 << "(" << getPolType() << ")" << endl
954 << setw(15) << "Channels:" << nchan() << endl;
955 String tmp;
956 oss << setw(15) << "Observer:"
957 << table_.keywordSet().asString("Observer") << endl;
958 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
959 table_.keywordSet().get("Project", tmp);
960 oss << setw(15) << "Project:" << tmp << endl;
961 table_.keywordSet().get("Obstype", tmp);
962 oss << setw(15) << "Obs. Type:" << tmp << endl;
963 table_.keywordSet().get("AntennaName", tmp);
964 oss << setw(15) << "Antenna Name:" << tmp << endl;
965 table_.keywordSet().get("FluxUnit", tmp);
966 oss << setw(15) << "Flux Unit:" << tmp << endl;
967 //Vector<Double> vec(moleculeTable_.getRestFrequencies());
968 int nid = moleculeTable_.nrow();
969 Bool firstline = True;
970 oss << setw(15) << "Rest Freqs:";
971 for (int i=0; i<nid; i++) {
972 Table t = table_(table_.col("MOLECULE_ID") == i);
973 if (t.nrow() > 0) {
974 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
975 if (vec.nelements() > 0) {
976 if (firstline) {
977 oss << setprecision(10) << vec << " [Hz]" << endl;
978 firstline=False;
979 }
980 else{
981 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
982 }
983 } else {
984 oss << "none" << endl;
985 }
986 }
987 }
988
989 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
990 oss << selector_.print() << endl;
991 oss << endl;
992 // main table
993 String dirtype = "Position ("
994 + getDirectionRefString()
995 + ")";
996 oss << setw(5) << "Scan" << setw(15) << "Source"
997 << setw(10) << "Time" << setw(18) << "Integration"
998 << setw(15) << "Source Type" << endl;
999 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
1000 oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
1001 << setw(8) << "Frame" << setw(16)
1002 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
1003 << setw(7) << "Channels"
1004 << endl;
1005 oss << asap::SEPERATOR << endl;
1006 TableIterator iter(table_, "SCANNO");
1007 while (!iter.pastEnd()) {
1008 Table subt = iter.table();
1009 ROTableRow row(subt);
1010 MEpoch::ROScalarColumn timeCol(subt,"TIME");
1011 const TableRecord& rec = row.get(0);
1012 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1013 << std::left << setw(1) << ""
1014 << setw(15) << rec.asString("SRCNAME")
1015 << setw(10) << formatTime(timeCol(0), false);
1016 // count the cycles in the scan
1017 TableIterator cyciter(subt, "CYCLENO");
1018 int nint = 0;
1019 while (!cyciter.pastEnd()) {
1020 ++nint;
1021 ++cyciter;
1022 }
1023 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
1024 << setw(11) << formatSec(rec.asFloat("INTERVAL")) << setw(1) << ""
1025 << setw(15) << SrcType::getName(rec.asInt("SRCTYPE")) << endl;
1026
1027 TableIterator biter(subt, "BEAMNO");
1028 while (!biter.pastEnd()) {
1029 Table bsubt = biter.table();
1030 ROTableRow brow(bsubt);
1031 const TableRecord& brec = brow.get(0);
1032 uInt row0 = bsubt.rowNumbers(table_)[0];
1033 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
1034 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
1035 TableIterator iiter(bsubt, "IFNO");
1036 while (!iiter.pastEnd()) {
1037 Table isubt = iiter.table();
1038 ROTableRow irow(isubt);
1039 const TableRecord& irec = irow.get(0);
1040 oss << setw(9) << "";
1041 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
1042 << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1043 << setw(3) << "" << nchan(irec.asuInt("IFNO"))
1044 << endl;
1045
1046 ++iiter;
1047 }
1048 ++biter;
1049 }
1050 ++iter;
1051 }
1052 /// @todo implement verbose mode
1053 return String(oss);
1054}
1055
1056// std::string Scantable::getTime(int whichrow, bool showdate) const
1057// {
1058// MEpoch::ROScalarColumn timeCol(table_, "TIME");
1059// MEpoch me;
1060// if (whichrow > -1) {
1061// me = timeCol(uInt(whichrow));
1062// } else {
1063// Double tm;
1064// table_.keywordSet().get("UTC",tm);
1065// me = MEpoch(MVEpoch(tm));
1066// }
1067// return formatTime(me, showdate);
1068// }
1069
1070std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
1071{
1072 MEpoch me;
1073 me = getEpoch(whichrow);
1074 return formatTime(me, showdate, prec);
1075}
1076
1077MEpoch Scantable::getEpoch(int whichrow) const
1078{
1079 if (whichrow > -1) {
1080 return timeCol_(uInt(whichrow));
1081 } else {
1082 Double tm;
1083 table_.keywordSet().get("UTC",tm);
1084 return MEpoch(MVEpoch(tm));
1085 }
1086}
1087
1088std::string Scantable::getDirectionString(int whichrow) const
1089{
1090 return formatDirection(getDirection(uInt(whichrow)));
1091}
1092
1093
1094SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1095 const MPosition& mp = getAntennaPosition();
1096 const MDirection& md = getDirection(whichrow);
1097 const MEpoch& me = timeCol_(whichrow);
1098 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1099 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1100 return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1101 mfreqidCol_(whichrow));
1102}
1103
1104std::vector< double > Scantable::getAbcissa( int whichrow ) const
1105{
1106 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
1107 std::vector<double> stlout;
1108 int nchan = specCol_(whichrow).nelements();
1109 String us = freqTable_.getUnitString();
1110 if ( us == "" || us == "pixel" || us == "channel" ) {
1111 for (int i=0; i<nchan; ++i) {
1112 stlout.push_back(double(i));
1113 }
1114 return stlout;
1115 }
1116 SpectralCoordinate spc = getSpectralCoordinate(whichrow);
1117 Vector<Double> pixel(nchan);
1118 Vector<Double> world;
1119 indgen(pixel);
1120 if ( Unit(us) == Unit("Hz") ) {
1121 for ( int i=0; i < nchan; ++i) {
1122 Double world;
1123 spc.toWorld(world, pixel[i]);
1124 stlout.push_back(double(world));
1125 }
1126 } else if ( Unit(us) == Unit("km/s") ) {
1127 Vector<Double> world;
1128 spc.pixelToVelocity(world, pixel);
1129 world.tovector(stlout);
1130 }
1131 return stlout;
1132}
1133void Scantable::setDirectionRefString( const std::string & refstr )
1134{
1135 MDirection::Types mdt;
1136 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1137 throw(AipsError("Illegal Direction frame."));
1138 }
1139 if ( refstr == "" ) {
1140 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1141 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1142 } else {
1143 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1144 }
1145}
1146
1147std::string Scantable::getDirectionRefString( ) const
1148{
1149 return table_.keywordSet().asString("DIRECTIONREF");
1150}
1151
1152MDirection Scantable::getDirection(int whichrow ) const
1153{
1154 String usertype = table_.keywordSet().asString("DIRECTIONREF");
1155 String type = MDirection::showType(dirCol_.getMeasRef().getType());
1156 if ( usertype != type ) {
1157 MDirection::Types mdt;
1158 if (!MDirection::getType(mdt, usertype)) {
1159 throw(AipsError("Illegal Direction frame."));
1160 }
1161 return dirCol_.convert(uInt(whichrow), mdt);
1162 } else {
1163 return dirCol_(uInt(whichrow));
1164 }
1165}
1166
1167std::string Scantable::getAbcissaLabel( int whichrow ) const
1168{
1169 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
1170 const MPosition& mp = getAntennaPosition();
1171 const MDirection& md = getDirection(whichrow);
1172 const MEpoch& me = timeCol_(whichrow);
1173 //const Double& rf = mmolidCol_(whichrow);
1174 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1175 SpectralCoordinate spc =
1176 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1177
1178 String s = "Channel";
1179 Unit u = Unit(freqTable_.getUnitString());
1180 if (u == Unit("km/s")) {
1181 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
1182 } else if (u == Unit("Hz")) {
1183 Vector<String> wau(1);wau = u.getName();
1184 spc.setWorldAxisUnits(wau);
1185 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
1186 }
1187 return s;
1188
1189}
1190
1191/**
1192void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
1193 const std::string& unit )
1194**/
1195void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1196 const std::string& unit )
1197
1198{
1199 ///@todo lookup in line table to fill in name and formattedname
1200 Unit u(unit);
1201 //Quantum<Double> urf(rf, u);
1202 Quantum<Vector<Double> >urf(rf, u);
1203 Vector<String> formattedname(0);
1204 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1205
1206 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1207 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
1208 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1209 tabvec = id;
1210}
1211
1212/**
1213void asap::Scantable::setRestFrequencies( const std::string& name )
1214{
1215 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1216 ///@todo implement
1217}
1218**/
1219void Scantable::setRestFrequencies( const vector<std::string>& name )
1220{
1221 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1222 ///@todo implement
1223}
1224
1225std::vector< unsigned int > Scantable::rownumbers( ) const
1226{
1227 std::vector<unsigned int> stlout;
1228 Vector<uInt> vec = table_.rowNumbers();
1229 vec.tovector(stlout);
1230 return stlout;
1231}
1232
1233
1234Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
1235{
1236 ROTableRow row(table_);
1237 const TableRecord& rec = row.get(whichrow);
1238 Table t =
1239 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1240 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1241 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1242 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1243 ROArrayColumn<Float> speccol(t, "SPECTRA");
1244 return speccol.getColumn();
1245}
1246
1247std::vector< std::string > Scantable::columnNames( ) const
1248{
1249 Vector<String> vec = table_.tableDesc().columnNames();
1250 return mathutil::tovectorstring(vec);
1251}
1252
1253MEpoch::Types Scantable::getTimeReference( ) const
1254{
1255 return MEpoch::castType(timeCol_.getMeasRef().getType());
1256}
1257
1258void Scantable::addFit( const STFitEntry& fit, int row )
1259{
1260 //cout << mfitidCol_(uInt(row)) << endl;
1261 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1262 os << mfitidCol_(uInt(row)) << LogIO::POST ;
1263 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1264 mfitidCol_.put(uInt(row), id);
1265}
1266
1267void Scantable::shift(int npix)
1268{
1269 Vector<uInt> fids(mfreqidCol_.getColumn());
1270 genSort( fids, Sort::Ascending,
1271 Sort::QuickSort|Sort::NoDuplicates );
1272 for (uInt i=0; i<fids.nelements(); ++i) {
1273 frequencies().shiftRefPix(npix, fids[i]);
1274 }
1275}
1276
1277String Scantable::getAntennaName() const
1278{
1279 String out;
1280 table_.keywordSet().get("AntennaName", out);
1281 String::size_type pos1 = out.find("@") ;
1282 String::size_type pos2 = out.find("//") ;
1283 if ( pos2 != String::npos )
1284 out = out.substr(pos2+2,pos1) ;
1285 else if ( pos1 != String::npos )
1286 out = out.substr(0,pos1) ;
1287 return out;
1288}
1289
1290int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
1291{
1292 String tbpath;
1293 int ret = 0;
1294 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1295 table_.keywordSet().get("GBT_GO", tbpath);
1296 Table t(tbpath,Table::Old);
1297 // check each scan if other scan of the pair exist
1298 int nscan = scanlist.size();
1299 for (int i = 0; i < nscan; i++) {
1300 Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1301 if (subt.nrow()==0) {
1302 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1303 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1304 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
1305 ret = 1;
1306 break;
1307 }
1308 ROTableRow row(subt);
1309 const TableRecord& rec = row.get(0);
1310 int scan1seqn = rec.asuInt("PROCSEQN");
1311 int laston1 = rec.asuInt("LASTON");
1312 if ( rec.asuInt("PROCSIZE")==2 ) {
1313 if ( i < nscan-1 ) {
1314 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1315 if ( subt2.nrow() == 0) {
1316 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1317
1318 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1319 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
1320 ret = 1;
1321 break;
1322 }
1323 ROTableRow row2(subt2);
1324 const TableRecord& rec2 = row2.get(0);
1325 int scan2seqn = rec2.asuInt("PROCSEQN");
1326 int laston2 = rec2.asuInt("LASTON");
1327 if (scan1seqn == 1 && scan2seqn == 2) {
1328 if (laston1 == laston2) {
1329 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1330 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1331 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1332 i +=1;
1333 }
1334 else {
1335 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1336 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1337 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
1338 }
1339 }
1340 else if (scan1seqn==2 && scan2seqn == 1) {
1341 if (laston1 == laston2) {
1342 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1343 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1344 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
1345 ret = 1;
1346 break;
1347 }
1348 }
1349 else {
1350 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1351 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1352 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
1353 ret = 1;
1354 break;
1355 }
1356 }
1357 }
1358 else {
1359 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1360 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1361 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
1362 }
1363 //if ( i >= nscan ) break;
1364 }
1365 }
1366 else {
1367 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1368 //cerr<<"No reference to GBT_GO table."<<endl;
1369 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
1370 ret = 1;
1371 }
1372 return ret;
1373}
1374
1375std::vector<double> Scantable::getDirectionVector(int whichrow) const
1376{
1377 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1378 std::vector<double> dir;
1379 Dir.tovector(dir);
1380 return dir;
1381}
1382
1383void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1384 throw( casa::AipsError )
1385{
1386 // assumed that all rows have same nChan
1387 Vector<Float> arr = specCol_( 0 ) ;
1388 int nChan = arr.nelements() ;
1389
1390 // if nmin < 0 or nmax < 0, nothing to do
1391 if ( nmin < 0 ) {
1392 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1393 }
1394 if ( nmax < 0 ) {
1395 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1396 }
1397
1398 // if nmin > nmax, exchange values
1399 if ( nmin > nmax ) {
1400 int tmp = nmax ;
1401 nmax = nmin ;
1402 nmin = tmp ;
1403 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1404 os << "Swap values. Applied range is ["
1405 << nmin << ", " << nmax << "]" << LogIO::POST ;
1406 }
1407
1408 // if nmin exceeds nChan, nothing to do
1409 if ( nmin >= nChan ) {
1410 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1411 }
1412
1413 // if nmax exceeds nChan, reset nmax to nChan
1414 if ( nmax >= nChan ) {
1415 if ( nmin == 0 ) {
1416 // nothing to do
1417 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1418 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1419 return ;
1420 }
1421 else {
1422 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1423 os << "Specified maximum exceeds nChan. Applied range is ["
1424 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1425 nmax = nChan - 1 ;
1426 }
1427 }
1428
1429 // reshape specCol_ and flagCol_
1430 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1431 reshapeSpectrum( nmin, nmax, irow ) ;
1432 }
1433
1434 // update FREQUENCIES subtable
1435 Double refpix ;
1436 Double refval ;
1437 Double increment ;
1438 int freqnrow = freqTable_.table().nrow() ;
1439 Vector<uInt> oldId( freqnrow ) ;
1440 Vector<uInt> newId( freqnrow ) ;
1441 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1442 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1443 /***
1444 * need to shift refpix to nmin
1445 * note that channel nmin in old index will be channel 0 in new one
1446 ***/
1447 refval = refval - ( refpix - nmin ) * increment ;
1448 refpix = 0 ;
1449 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1450 }
1451
1452 // update nchan
1453 int newsize = nmax - nmin + 1 ;
1454 table_.rwKeywordSet().define( "nChan", newsize ) ;
1455
1456 // update bandwidth
1457 // assumed all spectra in the scantable have same bandwidth
1458 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1459
1460 return ;
1461}
1462
1463void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1464{
1465 // reshape specCol_ and flagCol_
1466 Vector<Float> oldspec = specCol_( irow ) ;
1467 Vector<uChar> oldflag = flagsCol_( irow ) ;
1468 uInt newsize = nmax - nmin + 1 ;
1469 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1470 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1471
1472 return ;
1473}
1474
1475void asap::Scantable::regridChannel( int nChan, double dnu )
1476{
1477 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1478 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1479 // assumed that all rows have same nChan
1480 Vector<Float> arr = specCol_( 0 ) ;
1481 int oldsize = arr.nelements() ;
1482
1483 // if oldsize == nChan, nothing to do
1484 if ( oldsize == nChan ) {
1485 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1486 return ;
1487 }
1488
1489 // if oldChan < nChan, unphysical operation
1490 if ( oldsize < nChan ) {
1491 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1492 return ;
1493 }
1494
1495 // change channel number for specCol_ and flagCol_
1496 Vector<Float> newspec( nChan, 0 ) ;
1497 Vector<uChar> newflag( nChan, false ) ;
1498 vector<string> coordinfo = getCoordInfo() ;
1499 string oldinfo = coordinfo[0] ;
1500 coordinfo[0] = "Hz" ;
1501 setCoordInfo( coordinfo ) ;
1502 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1503 regridChannel( nChan, dnu, irow ) ;
1504 }
1505 coordinfo[0] = oldinfo ;
1506 setCoordInfo( coordinfo ) ;
1507
1508
1509 // NOTE: this method does not update metadata such as
1510 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1511
1512 return ;
1513}
1514
1515void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1516{
1517 // logging
1518 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1519 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1520
1521 Vector<Float> oldspec = specCol_( irow ) ;
1522 Vector<uChar> oldflag = flagsCol_( irow ) ;
1523 Vector<Float> newspec( nChan, 0 ) ;
1524 Vector<uChar> newflag( nChan, false ) ;
1525
1526 // regrid
1527 vector<double> abcissa = getAbcissa( irow ) ;
1528 int oldsize = abcissa.size() ;
1529 double olddnu = abcissa[1] - abcissa[0] ;
1530 //int refChan = 0 ;
1531 //double frac = 0.0 ;
1532 //double wedge = 0.0 ;
1533 //double pile = 0.0 ;
1534 int ichan = 0 ;
1535 double wsum = 0.0 ;
1536 Vector<Float> z( nChan ) ;
1537 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
1538 for ( int ii = 1 ; ii < nChan ; ii++ )
1539 z[ii] = z[ii-1] + dnu ;
1540 Vector<Float> zi( nChan+1 ) ;
1541 Vector<Float> yi( oldsize + 1 ) ;
1542 zi[0] = z[0] - 0.5 * dnu ;
1543 zi[1] = z[0] + 0.5 * dnu ;
1544 for ( int ii = 2 ; ii < nChan ; ii++ )
1545 zi[ii] = zi[ii-1] + dnu ;
1546 zi[nChan] = z[nChan-1] + 0.5 * dnu ;
1547 yi[0] = abcissa[0] - 0.5 * olddnu ;
1548 yi[1] = abcissa[1] + 0.5 * olddnu ;
1549 for ( int ii = 2 ; ii < oldsize ; ii++ )
1550 yi[ii] = abcissa[ii-1] + olddnu ;
1551 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
1552 if ( dnu > 0.0 ) {
1553 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1554 double zl = zi[ii] ;
1555 double zr = zi[ii+1] ;
1556 for ( int j = ichan ; j < oldsize ; j++ ) {
1557 double yl = yi[j] ;
1558 double yr = yi[j+1] ;
1559 if ( yl <= zl ) {
1560 if ( yr <= zl ) {
1561 continue ;
1562 }
1563 else if ( yr <= zr ) {
1564 newspec[ii] += oldspec[j] * ( yr - zl ) ;
1565 newflag[ii] = newflag[ii] || oldflag[j] ;
1566 wsum += ( yr - zl ) ;
1567 }
1568 else {
1569 newspec[ii] += oldspec[j] * dnu ;
1570 newflag[ii] = newflag[ii] || oldflag[j] ;
1571 wsum += dnu ;
1572 ichan = j ;
1573 break ;
1574 }
1575 }
1576 else if ( yl < zr ) {
1577 if ( yr <= zr ) {
1578 newspec[ii] += oldspec[j] * ( yr - yl ) ;
1579 newflag[ii] = newflag[ii] || oldflag[j] ;
1580 wsum += ( yr - yl ) ;
1581 }
1582 else {
1583 newspec[ii] += oldspec[j] * ( zr - yl ) ;
1584 newflag[ii] = newflag[ii] || oldflag[j] ;
1585 wsum += ( zr - yl ) ;
1586 ichan = j ;
1587 break ;
1588 }
1589 }
1590 else {
1591 ichan = j - 1 ;
1592 break ;
1593 }
1594 }
1595 newspec[ii] /= wsum ;
1596 wsum = 0.0 ;
1597 }
1598 }
1599 else if ( dnu < 0.0 ) {
1600 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1601 double zl = zi[ii] ;
1602 double zr = zi[ii+1] ;
1603 for ( int j = ichan ; j < oldsize ; j++ ) {
1604 double yl = yi[j] ;
1605 double yr = yi[j+1] ;
1606 if ( yl >= zl ) {
1607 if ( yr >= zl ) {
1608 continue ;
1609 }
1610 else if ( yr >= zr ) {
1611 newspec[ii] += oldspec[j] * abs( yr - zl ) ;
1612 newflag[ii] = newflag[ii] || oldflag[j] ;
1613 wsum += abs( yr - zl ) ;
1614 }
1615 else {
1616 newspec[ii] += oldspec[j] * abs( dnu ) ;
1617 newflag[ii] = newflag[ii] || oldflag[j] ;
1618 wsum += abs( dnu ) ;
1619 ichan = j ;
1620 break ;
1621 }
1622 }
1623 else if ( yl > zr ) {
1624 if ( yr >= zr ) {
1625 newspec[ii] += oldspec[j] * abs( yr - yl ) ;
1626 newflag[ii] = newflag[ii] || oldflag[j] ;
1627 wsum += abs( yr - yl ) ;
1628 }
1629 else {
1630 newspec[ii] += oldspec[j] * abs( zr - yl ) ;
1631 newflag[ii] = newflag[ii] || oldflag[j] ;
1632 wsum += abs( zr - yl ) ;
1633 ichan = j ;
1634 break ;
1635 }
1636 }
1637 else {
1638 ichan = j - 1 ;
1639 break ;
1640 }
1641 }
1642 newspec[ii] /= wsum ;
1643 wsum = 0.0 ;
1644 }
1645 }
1646// * ichan = 0
1647// ***/
1648// //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
1649// pile += dnu ;
1650// wedge = olddnu * ( refChan + 1 ) ;
1651// while ( wedge < pile ) {
1652// newspec[0] += olddnu * oldspec[refChan] ;
1653// newflag[0] = newflag[0] || oldflag[refChan] ;
1654// //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
1655// refChan++ ;
1656// wedge += olddnu ;
1657// wsum += olddnu ;
1658// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1659// }
1660// frac = ( wedge - pile ) / olddnu ;
1661// wsum += ( 1.0 - frac ) * olddnu ;
1662// newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1663// newflag[0] = newflag[0] || oldflag[refChan] ;
1664// //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
1665// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1666// newspec[0] /= wsum ;
1667// //ofs << "newspec[0] = " << newspec[0] << endl ;
1668// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1669
1670// /***
1671// * ichan = 1 - nChan-2
1672// ***/
1673// for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
1674// pile += dnu ;
1675// newspec[ichan] += frac * olddnu * oldspec[refChan] ;
1676// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1677// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
1678// refChan++ ;
1679// wedge += olddnu ;
1680// wsum = frac * olddnu ;
1681// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1682// while ( wedge < pile ) {
1683// newspec[ichan] += olddnu * oldspec[refChan] ;
1684// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1685// //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
1686// refChan++ ;
1687// wedge += olddnu ;
1688// wsum += olddnu ;
1689// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1690// }
1691// frac = ( wedge - pile ) / olddnu ;
1692// wsum += ( 1.0 - frac ) * olddnu ;
1693// newspec[ichan] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1694// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1695// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << ( 1.0 - frac ) << endl ;
1696// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1697// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1698// newspec[ichan] /= wsum ;
1699// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
1700// }
1701
1702// /***
1703// * ichan = nChan-1
1704// ***/
1705// // NOTE: Assumed that all spectra have the same bandwidth
1706// pile += dnu ;
1707// newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
1708// newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
1709// //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1710// refChan++ ;
1711// wedge += olddnu ;
1712// wsum = frac * olddnu ;
1713// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1714// for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
1715// newspec[nChan-1] += olddnu * oldspec[jchan] ;
1716// newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
1717// wsum += olddnu ;
1718// //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1719// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1720// }
1721// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1722// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1723// newspec[nChan-1] /= wsum ;
1724// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
1725
1726// specCol_.put( irow, newspec ) ;
1727// flagsCol_.put( irow, newflag ) ;
1728
1729// // ofs.close() ;
1730
1731
1732 return ;
1733}
1734
1735std::vector<float> Scantable::getWeather(int whichrow) const
1736{
1737 std::vector<float> out(5);
1738 //Float temperature, pressure, humidity, windspeed, windaz;
1739 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
1740 mweatheridCol_(uInt(whichrow)));
1741
1742
1743 return out;
1744}
1745
1746bool Scantable::getFlagtraFast(int whichrow)
1747{
1748 uChar flag;
1749 Vector<uChar> flags;
1750 flagsCol_.get(uInt(whichrow), flags);
1751 for (int i = 0; i < flags.size(); i++) {
1752 if (i==0) {
1753 flag = flags[i];
1754 }
1755 else {
1756 flag &= flags[i];
1757 }
1758 return ((flag >> 7) == 1);
1759 }
1760}
1761
1762void Scantable::doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter)
1763{
1764 fitter.setExpression("poly", order);
1765
1766 std::vector<double> abcsd = getAbcissa(rowno);
1767 std::vector<float> abcs;
1768 for (int i = 0; i < abcsd.size(); i++) {
1769 abcs.push_back((float)abcsd[i]);
1770 }
1771 std::vector<float> spec = getSpectrum(rowno);
1772 std::vector<bool> fmask = getMask(rowno);
1773 if (fmask.size() != mask.size()) {
1774 throw(AipsError("different mask sizes"));
1775 }
1776 for (int i = 0; i < fmask.size(); i++) {
1777 fmask[i] = fmask[i] && mask[i];
1778 }
1779 fitter.setData(abcs, spec, fmask);
1780
1781 fitter.lfit();
1782}
1783
1784void Scantable::polyBaselineBatch(const std::vector<bool>& mask, int order)
1785{
1786 Fitter fitter = Fitter();
1787 for (uInt rowno=0; rowno < nrow(); ++rowno) {
1788 doPolyBaseline(mask, order, rowno, fitter);
1789 setSpectrum(fitter.getResidual(), rowno);
1790 }
1791}
1792
1793STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
1794{
1795 Fitter fitter = Fitter();
1796 doPolyBaseline(mask, order, rowno, fitter);
1797 setSpectrum(fitter.getResidual(), rowno);
1798 return fitter.getFitEntry();
1799}
1800
1801}
1802//namespace asap
Note: See TracBrowser for help on using the repository browser.