source: branches/polybatch/src/Scantable.cpp@ 2788

Last change on this file since 2788 was 1924, checked in by Malte Marquarding, 14 years ago

Ticket #206: use STFitEntry as return objetc instead of pointer wrnagling.

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