source: trunk/src/Scantable.cpp@ 1958

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:
+ An optional parameter 'prec (unsigned int)' is added to

scantable.get_time, python_Scantable::_gettime, ScantableWrapper::getTime and Scantable::getTime.

+ Also Scantable::fromatTime accepts 'prec' as a parameter.
+ scantable._get_column accepts args which will be passed to callback function.

Test Programs:

Put in Release Notes: No

Module(s): scantable

Description:

Add a parameter prec to scantable.get_time which specifies the precision of time returned.
The default value is prec=0.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.1 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{
[1947]686 return formatTime(me, showdate, 0);
687}
688
689std::string Scantable::formatTime(const MEpoch& me, bool showdate, uInt prec) const
690{
[805]691 MVTime mvt(me.getValue());
692 if (showdate)
[1947]693 //mvt.setFormat(MVTime::YMD);
694 mvt.setFormat(MVTime::YMD, prec);
[805]695 else
[1947]696 //mvt.setFormat(MVTime::TIME);
697 mvt.setFormat(MVTime::TIME, prec);
[805]698 ostringstream oss;
699 oss << mvt;
700 return String(oss);
701}
[488]702
[805]703void Scantable::calculateAZEL()
704{
705 MPosition mp = getAntennaPosition();
706 MEpoch::ROScalarColumn timeCol(table_, "TIME");
707 ostringstream oss;
708 oss << "Computed azimuth/elevation using " << endl
709 << mp << endl;
[996]710 for (Int i=0; i<nrow(); ++i) {
[805]711 MEpoch me = timeCol(i);
[987]712 MDirection md = getDirection(i);
[805]713 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
714 << endl << " => ";
715 MeasFrame frame(mp, me);
716 Vector<Double> azel =
717 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
718 frame)
719 )().getAngle("rad").getValue();
[923]720 azCol_.put(i,Float(azel[0]));
721 elCol_.put(i,Float(azel[1]));
[805]722 oss << "azel: " << azel[0]/C::pi*180.0 << " "
723 << azel[1]/C::pi*180.0 << " (deg)" << endl;
[16]724 }
[805]725 pushLog(String(oss));
726}
[89]727
[1819]728void Scantable::clip(const Float uthres, const Float dthres, bool clipoutside, bool unflag)
729{
730 for (uInt i=0; i<table_.nrow(); ++i) {
731 Vector<uChar> flgs = flagsCol_(i);
732 srchChannelsToClip(i, uthres, dthres, clipoutside, unflag, flgs);
733 flagsCol_.put(i, flgs);
734 }
735}
736
737std::vector<bool> Scantable::getClipMask(int whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag)
738{
739 Vector<uChar> flags;
740 flagsCol_.get(uInt(whichrow), flags);
741 srchChannelsToClip(uInt(whichrow), uthres, dthres, clipoutside, unflag, flags);
742 Vector<Bool> bflag(flags.shape());
743 convertArray(bflag, flags);
744 //bflag = !bflag;
745
746 std::vector<bool> mask;
747 bflag.tovector(mask);
748 return mask;
749}
750
751void Scantable::srchChannelsToClip(uInt whichrow, const Float uthres, const Float dthres, bool clipoutside, bool unflag,
752 Vector<uChar> flgs)
753{
754 Vector<Float> spcs = specCol_(whichrow);
755 uInt nchannel = nchan();
756 if (spcs.nelements() != nchannel) {
757 throw(AipsError("Data has incorrect number of channels"));
758 }
759 uChar userflag = 1 << 7;
760 if (unflag) {
761 userflag = 0 << 7;
762 }
763 if (clipoutside) {
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 } else {
771 for (uInt j = 0; j < nchannel; ++j) {
772 Float spc = spcs(j);
773 if ((spc < uthres) && (spc > dthres)) {
774 flgs(j) = userflag;
775 }
776 }
777 }
778}
779
[1430]780void Scantable::flag(const std::vector<bool>& msk, bool unflag)
[865]781{
[1333]782 std::vector<bool>::const_iterator it;
783 uInt ntrue = 0;
784 for (it = msk.begin(); it != msk.end(); ++it) {
785 if ( *it ) {
786 ntrue++;
787 }
788 }
789 if ( selector_.empty() && (msk.size() == 0 || msk.size() == ntrue) )
[1000]790 throw(AipsError("Trying to flag whole scantable."));
791 if ( msk.size() == 0 ) {
792 uChar userflag = 1 << 7;
[1430]793 if ( unflag ) {
794 userflag = 0 << 7;
795 }
[1000]796 for ( uInt i=0; i<table_.nrow(); ++i) {
797 Vector<uChar> flgs = flagsCol_(i);
798 flgs = userflag;
799 flagsCol_.put(i, flgs);
800 }
801 return;
802 }
803 if ( int(msk.size()) != nchan() ) {
804 throw(AipsError("Mask has incorrect number of channels."));
805 }
806 for ( uInt i=0; i<table_.nrow(); ++i) {
807 Vector<uChar> flgs = flagsCol_(i);
808 if ( flgs.nelements() != msk.size() ) {
809 throw(AipsError("Mask has incorrect number of channels."
810 " Probably varying with IF. Please flag per IF"));
811 }
812 std::vector<bool>::const_iterator it;
813 uInt j = 0;
814 uChar userflag = 1 << 7;
[1430]815 if ( unflag ) {
816 userflag = 0 << 7;
817 }
[1000]818 for (it = msk.begin(); it != msk.end(); ++it) {
819 if ( *it ) {
820 flgs(j) = userflag;
821 }
822 ++j;
823 }
824 flagsCol_.put(i, flgs);
825 }
[865]826}
827
[1819]828void Scantable::flagRow(const std::vector<uInt>& rows, bool unflag)
829{
830 if ( selector_.empty() && (rows.size() == table_.nrow()) )
831 throw(AipsError("Trying to flag whole scantable."));
832
833 uInt rowflag = (unflag ? 0 : 1);
834 std::vector<uInt>::const_iterator it;
835 for (it = rows.begin(); it != rows.end(); ++it)
836 flagrowCol_.put(*it, rowflag);
837}
838
[805]839std::vector<bool> Scantable::getMask(int whichrow) const
840{
841 Vector<uChar> flags;
842 flagsCol_.get(uInt(whichrow), flags);
843 Vector<Bool> bflag(flags.shape());
844 convertArray(bflag, flags);
845 bflag = !bflag;
846 std::vector<bool> mask;
847 bflag.tovector(mask);
848 return mask;
849}
[89]850
[896]851std::vector<float> Scantable::getSpectrum( int whichrow,
[902]852 const std::string& poltype ) const
[805]853{
[905]854 String ptype = poltype;
855 if (poltype == "" ) ptype = getPolType();
[902]856 if ( whichrow < 0 || whichrow >= nrow() )
857 throw(AipsError("Illegal row number."));
[896]858 std::vector<float> out;
[805]859 Vector<Float> arr;
[896]860 uInt requestedpol = polCol_(whichrow);
861 String basetype = getPolType();
[905]862 if ( ptype == basetype ) {
[896]863 specCol_.get(whichrow, arr);
864 } else {
[1598]865 CountedPtr<STPol> stpol(STPol::getPolClass(Scantable::factories_,
[1586]866 basetype));
[1334]867 uInt row = uInt(whichrow);
868 stpol->setSpectra(getPolMatrix(row));
869 Float fang,fhand,parang;
[1586]870 fang = focusTable_.getTotalAngle(mfocusidCol_(row));
[1334]871 fhand = focusTable_.getFeedHand(mfocusidCol_(row));
[1586]872 stpol->setPhaseCorrections(fang, fhand);
[1334]873 arr = stpol->getSpectrum(requestedpol, ptype);
[896]874 }
[902]875 if ( arr.nelements() == 0 )
876 pushLog("Not enough polarisations present to do the conversion.");
[805]877 arr.tovector(out);
878 return out;
[89]879}
[212]880
[1360]881void Scantable::setSpectrum( const std::vector<float>& spec,
[884]882 int whichrow )
883{
884 Vector<Float> spectrum(spec);
885 Vector<Float> arr;
886 specCol_.get(whichrow, arr);
887 if ( spectrum.nelements() != arr.nelements() )
[896]888 throw AipsError("The spectrum has incorrect number of channels.");
[884]889 specCol_.put(whichrow, spectrum);
890}
891
892
[805]893String Scantable::generateName()
[745]894{
[805]895 return (File::newUniqueName("./","temp")).baseName();
[212]896}
897
[805]898const casa::Table& Scantable::table( ) const
[212]899{
[805]900 return table_;
[212]901}
902
[805]903casa::Table& Scantable::table( )
[386]904{
[805]905 return table_;
[386]906}
907
[896]908std::string Scantable::getPolType() const
909{
910 return table_.keywordSet().asString("POLTYPE");
911}
912
[805]913void Scantable::unsetSelection()
[380]914{
[805]915 table_ = originalTable_;
[847]916 attach();
[805]917 selector_.reset();
[380]918}
[386]919
[805]920void Scantable::setSelection( const STSelector& selection )
[430]921{
[805]922 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
923 if ( tab.nrow() == 0 ) {
924 throw(AipsError("Selection contains no data. Not applying it."));
925 }
926 table_ = tab;
[847]927 attach();
[805]928 selector_ = selection;
[430]929}
930
[837]931std::string Scantable::summary( bool verbose )
[447]932{
[805]933 // Format header info
934 ostringstream oss;
935 oss << endl;
936 oss << asap::SEPERATOR << endl;
937 oss << " Scan Table Summary" << endl;
938 oss << asap::SEPERATOR << endl;
939 oss.flags(std::ios_base::left);
940 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
941 << setw(15) << "IFs:" << setw(4) << nif() << endl
[896]942 << setw(15) << "Polarisations:" << setw(4) << npol()
943 << "(" << getPolType() << ")" << endl
[1694]944 << setw(15) << "Channels:" << nchan() << endl;
[805]945 String tmp;
[860]946 oss << setw(15) << "Observer:"
947 << table_.keywordSet().asString("Observer") << endl;
[805]948 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
949 table_.keywordSet().get("Project", tmp);
950 oss << setw(15) << "Project:" << tmp << endl;
951 table_.keywordSet().get("Obstype", tmp);
952 oss << setw(15) << "Obs. Type:" << tmp << endl;
953 table_.keywordSet().get("AntennaName", tmp);
954 oss << setw(15) << "Antenna Name:" << tmp << endl;
955 table_.keywordSet().get("FluxUnit", tmp);
956 oss << setw(15) << "Flux Unit:" << tmp << endl;
[1819]957 //Vector<Double> vec(moleculeTable_.getRestFrequencies());
958 int nid = moleculeTable_.nrow();
959 Bool firstline = True;
[805]960 oss << setw(15) << "Rest Freqs:";
[1819]961 for (int i=0; i<nid; i++) {
962 Table t = table_(table_.col("MOLECULE_ID") == i);
963 if (t.nrow() > 0) {
964 Vector<Double> vec(moleculeTable_.getRestFrequency(i));
965 if (vec.nelements() > 0) {
966 if (firstline) {
967 oss << setprecision(10) << vec << " [Hz]" << endl;
968 firstline=False;
969 }
970 else{
971 oss << setw(15)<<" " << setprecision(10) << vec << " [Hz]" << endl;
972 }
973 } else {
974 oss << "none" << endl;
975 }
976 }
[805]977 }
[941]978
979 oss << setw(15) << "Abcissa:" << getAbcissaLabel(0) << endl;
[805]980 oss << selector_.print() << endl;
981 oss << endl;
982 // main table
983 String dirtype = "Position ("
[987]984 + getDirectionRefString()
[805]985 + ")";
[941]986 oss << setw(5) << "Scan" << setw(15) << "Source"
987 << setw(10) << "Time" << setw(18) << "Integration" << endl;
988 oss << setw(5) << "" << setw(5) << "Beam" << setw(3) << "" << dirtype << endl;
[1694]989 oss << setw(10) << "" << setw(3) << "IF" << setw(3) << ""
[805]990 << setw(8) << "Frame" << setw(16)
[1694]991 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment"
992 << setw(7) << "Channels"
993 << endl;
[805]994 oss << asap::SEPERATOR << endl;
995 TableIterator iter(table_, "SCANNO");
996 while (!iter.pastEnd()) {
997 Table subt = iter.table();
998 ROTableRow row(subt);
999 MEpoch::ROScalarColumn timeCol(subt,"TIME");
1000 const TableRecord& rec = row.get(0);
1001 oss << setw(4) << std::right << rec.asuInt("SCANNO")
1002 << std::left << setw(1) << ""
1003 << setw(15) << rec.asString("SRCNAME")
1004 << setw(10) << formatTime(timeCol(0), false);
1005 // count the cycles in the scan
1006 TableIterator cyciter(subt, "CYCLENO");
1007 int nint = 0;
1008 while (!cyciter.pastEnd()) {
1009 ++nint;
1010 ++cyciter;
1011 }
1012 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
1013 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl;
[447]1014
[805]1015 TableIterator biter(subt, "BEAMNO");
1016 while (!biter.pastEnd()) {
1017 Table bsubt = biter.table();
1018 ROTableRow brow(bsubt);
1019 const TableRecord& brec = brow.get(0);
[1000]1020 uInt row0 = bsubt.rowNumbers(table_)[0];
[941]1021 oss << setw(5) << "" << setw(4) << std::right << brec.asuInt("BEAMNO")<< std::left;
[987]1022 oss << setw(4) << "" << formatDirection(getDirection(row0)) << endl;
[805]1023 TableIterator iiter(bsubt, "IFNO");
1024 while (!iiter.pastEnd()) {
1025 Table isubt = iiter.table();
1026 ROTableRow irow(isubt);
1027 const TableRecord& irec = irow.get(0);
[1694]1028 oss << setw(9) << "";
[941]1029 oss << setw(3) << std::right << irec.asuInt("IFNO") << std::left
[1694]1030 << setw(1) << "" << frequencies().print(irec.asuInt("FREQ_ID"))
1031 << setw(3) << "" << nchan(irec.asuInt("IFNO"))
[1375]1032 << endl;
[447]1033
[805]1034 ++iiter;
1035 }
1036 ++biter;
1037 }
1038 ++iter;
[447]1039 }
[805]1040 /// @todo implement verbose mode
1041 return String(oss);
[447]1042}
1043
[1947]1044// std::string Scantable::getTime(int whichrow, bool showdate) const
1045// {
1046// MEpoch::ROScalarColumn timeCol(table_, "TIME");
1047// MEpoch me;
1048// if (whichrow > -1) {
1049// me = timeCol(uInt(whichrow));
1050// } else {
1051// Double tm;
1052// table_.keywordSet().get("UTC",tm);
1053// me = MEpoch(MVEpoch(tm));
1054// }
1055// return formatTime(me, showdate);
1056// }
1057
1058std::string Scantable::getTime(int whichrow, bool showdate, uInt prec) const
[777]1059{
[805]1060 MEpoch me;
[1947]1061 me = getEpoch(whichrow);
1062 return formatTime(me, showdate, prec);
[777]1063}
[805]1064
[1411]1065MEpoch Scantable::getEpoch(int whichrow) const
1066{
1067 if (whichrow > -1) {
1068 return timeCol_(uInt(whichrow));
1069 } else {
1070 Double tm;
1071 table_.keywordSet().get("UTC",tm);
[1598]1072 return MEpoch(MVEpoch(tm));
[1411]1073 }
1074}
1075
[1068]1076std::string Scantable::getDirectionString(int whichrow) const
1077{
1078 return formatDirection(getDirection(uInt(whichrow)));
1079}
1080
[1598]1081
1082SpectralCoordinate Scantable::getSpectralCoordinate(int whichrow) const {
1083 const MPosition& mp = getAntennaPosition();
1084 const MDirection& md = getDirection(whichrow);
1085 const MEpoch& me = timeCol_(whichrow);
[1819]1086 //Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
1087 Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
[1598]1088 return freqTable_.getSpectralCoordinate(md, mp, me, rf,
1089 mfreqidCol_(whichrow));
1090}
1091
[1360]1092std::vector< double > Scantable::getAbcissa( int whichrow ) const
[865]1093{
[1507]1094 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal row number"));
[865]1095 std::vector<double> stlout;
1096 int nchan = specCol_(whichrow).nelements();
1097 String us = freqTable_.getUnitString();
1098 if ( us == "" || us == "pixel" || us == "channel" ) {
1099 for (int i=0; i<nchan; ++i) {
1100 stlout.push_back(double(i));
1101 }
1102 return stlout;
1103 }
[1598]1104 SpectralCoordinate spc = getSpectralCoordinate(whichrow);
[865]1105 Vector<Double> pixel(nchan);
1106 Vector<Double> world;
1107 indgen(pixel);
1108 if ( Unit(us) == Unit("Hz") ) {
1109 for ( int i=0; i < nchan; ++i) {
1110 Double world;
1111 spc.toWorld(world, pixel[i]);
1112 stlout.push_back(double(world));
1113 }
1114 } else if ( Unit(us) == Unit("km/s") ) {
1115 Vector<Double> world;
1116 spc.pixelToVelocity(world, pixel);
1117 world.tovector(stlout);
1118 }
1119 return stlout;
1120}
[1360]1121void Scantable::setDirectionRefString( const std::string & refstr )
[987]1122{
1123 MDirection::Types mdt;
1124 if (refstr != "" && !MDirection::getType(mdt, refstr)) {
1125 throw(AipsError("Illegal Direction frame."));
1126 }
1127 if ( refstr == "" ) {
1128 String defaultstr = MDirection::showType(dirCol_.getMeasRef().getType());
1129 table_.rwKeywordSet().define("DIRECTIONREF", defaultstr);
1130 } else {
1131 table_.rwKeywordSet().define("DIRECTIONREF", String(refstr));
1132 }
1133}
[865]1134
[1360]1135std::string Scantable::getDirectionRefString( ) const
[987]1136{
1137 return table_.keywordSet().asString("DIRECTIONREF");
1138}
1139
1140MDirection Scantable::getDirection(int whichrow ) const
1141{
1142 String usertype = table_.keywordSet().asString("DIRECTIONREF");
1143 String type = MDirection::showType(dirCol_.getMeasRef().getType());
1144 if ( usertype != type ) {
1145 MDirection::Types mdt;
1146 if (!MDirection::getType(mdt, usertype)) {
1147 throw(AipsError("Illegal Direction frame."));
1148 }
1149 return dirCol_.convert(uInt(whichrow), mdt);
1150 } else {
1151 return dirCol_(uInt(whichrow));
1152 }
1153}
1154
[847]1155std::string Scantable::getAbcissaLabel( int whichrow ) const
1156{
[996]1157 if ( whichrow > int(table_.nrow()) ) throw(AipsError("Illegal ro number"));
[847]1158 const MPosition& mp = getAntennaPosition();
[987]1159 const MDirection& md = getDirection(whichrow);
[847]1160 const MEpoch& me = timeCol_(whichrow);
[1819]1161 //const Double& rf = mmolidCol_(whichrow);
1162 const Vector<Double> rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
[847]1163 SpectralCoordinate spc =
1164 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
1165
1166 String s = "Channel";
1167 Unit u = Unit(freqTable_.getUnitString());
1168 if (u == Unit("km/s")) {
[1170]1169 s = CoordinateUtil::axisLabel(spc, 0, True,True, True);
[847]1170 } else if (u == Unit("Hz")) {
1171 Vector<String> wau(1);wau = u.getName();
1172 spc.setWorldAxisUnits(wau);
[1170]1173 s = CoordinateUtil::axisLabel(spc, 0, True, True, False);
[847]1174 }
1175 return s;
1176
1177}
1178
[1819]1179/**
1180void asap::Scantable::setRestFrequencies( double rf, const std::string& name,
[1170]1181 const std::string& unit )
[1819]1182**/
1183void Scantable::setRestFrequencies( vector<double> rf, const vector<std::string>& name,
1184 const std::string& unit )
1185
[847]1186{
[923]1187 ///@todo lookup in line table to fill in name and formattedname
[847]1188 Unit u(unit);
[1819]1189 //Quantum<Double> urf(rf, u);
1190 Quantum<Vector<Double> >urf(rf, u);
1191 Vector<String> formattedname(0);
1192 //cerr<<"Scantable::setRestFrequnecies="<<urf<<endl;
1193
1194 //uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), name, "");
1195 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), mathutil::toVectorString(name), formattedname);
[847]1196 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
1197 tabvec = id;
1198}
1199
[1819]1200/**
1201void asap::Scantable::setRestFrequencies( const std::string& name )
[847]1202{
1203 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
1204 ///@todo implement
1205}
[1819]1206**/
1207void Scantable::setRestFrequencies( const vector<std::string>& name )
1208{
1209 throw(AipsError("setRestFrequencies( const vector<std::string>& name ) NYI"));
1210 ///@todo implement
1211}
[847]1212
[1360]1213std::vector< unsigned int > Scantable::rownumbers( ) const
[852]1214{
1215 std::vector<unsigned int> stlout;
1216 Vector<uInt> vec = table_.rowNumbers();
1217 vec.tovector(stlout);
1218 return stlout;
1219}
1220
[865]1221
[1360]1222Matrix<Float> Scantable::getPolMatrix( uInt whichrow ) const
[896]1223{
1224 ROTableRow row(table_);
1225 const TableRecord& rec = row.get(whichrow);
1226 Table t =
1227 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
1228 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1229 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
1230 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
1231 ROArrayColumn<Float> speccol(t, "SPECTRA");
1232 return speccol.getColumn();
1233}
[865]1234
[1360]1235std::vector< std::string > Scantable::columnNames( ) const
[902]1236{
1237 Vector<String> vec = table_.tableDesc().columnNames();
1238 return mathutil::tovectorstring(vec);
1239}
[896]1240
[1360]1241MEpoch::Types Scantable::getTimeReference( ) const
[915]1242{
1243 return MEpoch::castType(timeCol_.getMeasRef().getType());
[972]1244}
[915]1245
[1360]1246void Scantable::addFit( const STFitEntry& fit, int row )
[972]1247{
[1819]1248 //cout << mfitidCol_(uInt(row)) << endl;
1249 LogIO os( LogOrigin( "Scantable", "addFit()", WHERE ) ) ;
1250 os << mfitidCol_(uInt(row)) << LogIO::POST ;
[972]1251 uInt id = fitTable_.addEntry(fit, mfitidCol_(uInt(row)));
1252 mfitidCol_.put(uInt(row), id);
1253}
[915]1254
[1360]1255void Scantable::shift(int npix)
1256{
1257 Vector<uInt> fids(mfreqidCol_.getColumn());
1258 genSort( fids, Sort::Ascending,
1259 Sort::QuickSort|Sort::NoDuplicates );
1260 for (uInt i=0; i<fids.nelements(); ++i) {
[1567]1261 frequencies().shiftRefPix(npix, fids[i]);
[1360]1262 }
1263}
[987]1264
[1819]1265String Scantable::getAntennaName() const
[1391]1266{
1267 String out;
1268 table_.keywordSet().get("AntennaName", out);
1269 return out;
[987]1270}
[1391]1271
[1730]1272int Scantable::checkScanInfo(const std::vector<int>& scanlist) const
[1391]1273{
1274 String tbpath;
1275 int ret = 0;
1276 if ( table_.keywordSet().isDefined("GBT_GO") ) {
1277 table_.keywordSet().get("GBT_GO", tbpath);
1278 Table t(tbpath,Table::Old);
1279 // check each scan if other scan of the pair exist
1280 int nscan = scanlist.size();
1281 for (int i = 0; i < nscan; i++) {
1282 Table subt = t( t.col("SCAN") == scanlist[i]+1 );
1283 if (subt.nrow()==0) {
[1819]1284 //cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
1285 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1286 os <<LogIO::WARN<<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<LogIO::POST;
[1391]1287 ret = 1;
1288 break;
1289 }
1290 ROTableRow row(subt);
1291 const TableRecord& rec = row.get(0);
1292 int scan1seqn = rec.asuInt("PROCSEQN");
1293 int laston1 = rec.asuInt("LASTON");
1294 if ( rec.asuInt("PROCSIZE")==2 ) {
1295 if ( i < nscan-1 ) {
1296 Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
1297 if ( subt2.nrow() == 0) {
[1819]1298 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1299
1300 //cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
1301 os<<LogIO::WARN<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<LogIO::POST;
[1391]1302 ret = 1;
1303 break;
1304 }
1305 ROTableRow row2(subt2);
1306 const TableRecord& rec2 = row2.get(0);
1307 int scan2seqn = rec2.asuInt("PROCSEQN");
1308 int laston2 = rec2.asuInt("LASTON");
1309 if (scan1seqn == 1 && scan2seqn == 2) {
1310 if (laston1 == laston2) {
[1819]1311 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1312 //cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1313 os<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
[1391]1314 i +=1;
1315 }
1316 else {
[1819]1317 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1318 //cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
1319 os<<LogIO::WARN<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<LogIO::POST;
[1391]1320 }
1321 }
1322 else if (scan1seqn==2 && scan2seqn == 1) {
1323 if (laston1 == laston2) {
[1819]1324 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1325 //cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
1326 os<<LogIO::WARN<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<LogIO::POST;
[1391]1327 ret = 1;
1328 break;
1329 }
1330 }
1331 else {
[1819]1332 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1333 //cerr<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
1334 os<<LogIO::WARN<<"The other scan for "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<LogIO::POST;
[1391]1335 ret = 1;
1336 break;
1337 }
1338 }
1339 }
1340 else {
[1819]1341 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1342 //cerr<<"The scan does not appear to be standard obsevation."<<endl;
1343 os<<LogIO::WARN<<"The scan does not appear to be standard obsevation."<<LogIO::POST;
[1391]1344 }
1345 //if ( i >= nscan ) break;
1346 }
1347 }
1348 else {
[1819]1349 LogIO os( LogOrigin( "Scantable", "checkScanInfo()", WHERE ) ) ;
1350 //cerr<<"No reference to GBT_GO table."<<endl;
1351 os<<LogIO::WARN<<"No reference to GBT_GO table."<<LogIO::POST;
[1391]1352 ret = 1;
1353 }
1354 return ret;
1355}
1356
[1730]1357std::vector<double> Scantable::getDirectionVector(int whichrow) const
[1391]1358{
1359 Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
1360 std::vector<double> dir;
1361 Dir.tovector(dir);
1362 return dir;
1363}
1364
[1819]1365void asap::Scantable::reshapeSpectrum( int nmin, int nmax )
1366 throw( casa::AipsError )
1367{
1368 // assumed that all rows have same nChan
1369 Vector<Float> arr = specCol_( 0 ) ;
1370 int nChan = arr.nelements() ;
1371
1372 // if nmin < 0 or nmax < 0, nothing to do
1373 if ( nmin < 0 ) {
1374 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1375 }
1376 if ( nmax < 0 ) {
1377 throw( casa::indexError<int>( nmax, "asap::Scantable::reshapeSpectrum: Invalid range. Negative index is specified." ) ) ;
1378 }
1379
1380 // if nmin > nmax, exchange values
1381 if ( nmin > nmax ) {
1382 int tmp = nmax ;
1383 nmax = nmin ;
1384 nmin = tmp ;
1385 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1386 os << "Swap values. Applied range is ["
1387 << nmin << ", " << nmax << "]" << LogIO::POST ;
1388 }
1389
1390 // if nmin exceeds nChan, nothing to do
1391 if ( nmin >= nChan ) {
1392 throw( casa::indexError<int>( nmin, "asap::Scantable::reshapeSpectrum: Invalid range. Specified minimum exceeds nChan." ) ) ;
1393 }
1394
1395 // if nmax exceeds nChan, reset nmax to nChan
1396 if ( nmax >= nChan ) {
1397 if ( nmin == 0 ) {
1398 // nothing to do
1399 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1400 os << "Whole range is selected. Nothing to do." << LogIO::POST ;
1401 return ;
1402 }
1403 else {
1404 LogIO os( LogOrigin( "Scantable", "reshapeSpectrum()", WHERE ) ) ;
1405 os << "Specified maximum exceeds nChan. Applied range is ["
1406 << nmin << ", " << nChan-1 << "]." << LogIO::POST ;
1407 nmax = nChan - 1 ;
1408 }
1409 }
1410
1411 // reshape specCol_ and flagCol_
1412 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1413 reshapeSpectrum( nmin, nmax, irow ) ;
1414 }
1415
1416 // update FREQUENCIES subtable
1417 Double refpix ;
1418 Double refval ;
1419 Double increment ;
1420 int freqnrow = freqTable_.table().nrow() ;
1421 Vector<uInt> oldId( freqnrow ) ;
1422 Vector<uInt> newId( freqnrow ) ;
1423 for ( int irow = 0 ; irow < freqnrow ; irow++ ) {
1424 freqTable_.getEntry( refpix, refval, increment, irow ) ;
1425 /***
1426 * need to shift refpix to nmin
1427 * note that channel nmin in old index will be channel 0 in new one
1428 ***/
1429 refval = refval - ( refpix - nmin ) * increment ;
1430 refpix = 0 ;
1431 freqTable_.setEntry( refpix, refval, increment, irow ) ;
1432 }
1433
1434 // update nchan
1435 int newsize = nmax - nmin + 1 ;
1436 table_.rwKeywordSet().define( "nChan", newsize ) ;
1437
1438 // update bandwidth
1439 // assumed all spectra in the scantable have same bandwidth
1440 table_.rwKeywordSet().define( "Bandwidth", increment * newsize ) ;
1441
1442 return ;
1443}
1444
1445void asap::Scantable::reshapeSpectrum( int nmin, int nmax, int irow )
1446{
1447 // reshape specCol_ and flagCol_
1448 Vector<Float> oldspec = specCol_( irow ) ;
1449 Vector<uChar> oldflag = flagsCol_( irow ) ;
1450 uInt newsize = nmax - nmin + 1 ;
1451 specCol_.put( irow, oldspec( Slice( nmin, newsize, 1 ) ) ) ;
1452 flagsCol_.put( irow, oldflag( Slice( nmin, newsize, 1 ) ) ) ;
1453
1454 return ;
1455}
1456
1457void asap::Scantable::regridChannel( int nChan, double dnu )
1458{
1459 LogIO os( LogOrigin( "Scantable", "regridChannel()", WHERE ) ) ;
1460 os << "Regrid abcissa with channel number " << nChan << " and spectral resoultion " << dnu << "Hz." << LogIO::POST ;
1461 // assumed that all rows have same nChan
1462 Vector<Float> arr = specCol_( 0 ) ;
1463 int oldsize = arr.nelements() ;
1464
1465 // if oldsize == nChan, nothing to do
1466 if ( oldsize == nChan ) {
1467 os << "Specified channel number is same as current one. Nothing to do." << LogIO::POST ;
1468 return ;
1469 }
1470
1471 // if oldChan < nChan, unphysical operation
1472 if ( oldsize < nChan ) {
1473 os << "Unphysical operation. Nothing to do." << LogIO::POST ;
1474 return ;
1475 }
1476
1477 // change channel number for specCol_ and flagCol_
1478 Vector<Float> newspec( nChan, 0 ) ;
1479 Vector<uChar> newflag( nChan, false ) ;
1480 vector<string> coordinfo = getCoordInfo() ;
1481 string oldinfo = coordinfo[0] ;
1482 coordinfo[0] = "Hz" ;
1483 setCoordInfo( coordinfo ) ;
1484 for ( int irow = 0 ; irow < nrow() ; irow++ ) {
1485 regridChannel( nChan, dnu, irow ) ;
1486 }
1487 coordinfo[0] = oldinfo ;
1488 setCoordInfo( coordinfo ) ;
1489
1490
1491 // NOTE: this method does not update metadata such as
1492 // FREQUENCIES subtable, nChan, Bandwidth, etc.
1493
1494 return ;
1495}
1496
1497void asap::Scantable::regridChannel( int nChan, double dnu, int irow )
1498{
1499 // logging
1500 //ofstream ofs( "average.log", std::ios::out | std::ios::app ) ;
1501 //ofs << "IFNO = " << getIF( irow ) << " irow = " << irow << endl ;
1502
1503 Vector<Float> oldspec = specCol_( irow ) ;
1504 Vector<uChar> oldflag = flagsCol_( irow ) ;
1505 Vector<Float> newspec( nChan, 0 ) ;
1506 Vector<uChar> newflag( nChan, false ) ;
1507
1508 // regrid
1509 vector<double> abcissa = getAbcissa( irow ) ;
1510 int oldsize = abcissa.size() ;
1511 double olddnu = abcissa[1] - abcissa[0] ;
1512 //int refChan = 0 ;
1513 //double frac = 0.0 ;
1514 //double wedge = 0.0 ;
1515 //double pile = 0.0 ;
1516 int ichan = 0 ;
1517 double wsum = 0.0 ;
1518 Vector<Float> z( nChan ) ;
1519 z[0] = abcissa[0] - 0.5 * olddnu + 0.5 * dnu ;
1520 for ( int ii = 1 ; ii < nChan ; ii++ )
1521 z[ii] = z[ii-1] + dnu ;
1522 Vector<Float> zi( nChan+1 ) ;
1523 Vector<Float> yi( oldsize + 1 ) ;
1524 zi[0] = z[0] - 0.5 * dnu ;
1525 zi[1] = z[0] + 0.5 * dnu ;
1526 for ( int ii = 2 ; ii < nChan ; ii++ )
1527 zi[ii] = zi[ii-1] + dnu ;
1528 zi[nChan] = z[nChan-1] + 0.5 * dnu ;
1529 yi[0] = abcissa[0] - 0.5 * olddnu ;
1530 yi[1] = abcissa[1] + 0.5 * olddnu ;
1531 for ( int ii = 2 ; ii < oldsize ; ii++ )
1532 yi[ii] = abcissa[ii-1] + olddnu ;
1533 yi[oldsize] = abcissa[oldsize-1] + 0.5 * olddnu ;
1534 if ( dnu > 0.0 ) {
1535 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1536 double zl = zi[ii] ;
1537 double zr = zi[ii+1] ;
1538 for ( int j = ichan ; j < oldsize ; j++ ) {
1539 double yl = yi[j] ;
1540 double yr = yi[j+1] ;
1541 if ( yl <= zl ) {
1542 if ( yr <= zl ) {
1543 continue ;
1544 }
1545 else if ( yr <= zr ) {
1546 newspec[ii] += oldspec[j] * ( yr - zl ) ;
1547 newflag[ii] = newflag[ii] || oldflag[j] ;
1548 wsum += ( yr - zl ) ;
1549 }
1550 else {
1551 newspec[ii] += oldspec[j] * dnu ;
1552 newflag[ii] = newflag[ii] || oldflag[j] ;
1553 wsum += dnu ;
1554 ichan = j ;
1555 break ;
1556 }
1557 }
1558 else if ( yl < zr ) {
1559 if ( yr <= zr ) {
1560 newspec[ii] += oldspec[j] * ( yr - yl ) ;
1561 newflag[ii] = newflag[ii] || oldflag[j] ;
1562 wsum += ( yr - yl ) ;
1563 }
1564 else {
1565 newspec[ii] += oldspec[j] * ( zr - yl ) ;
1566 newflag[ii] = newflag[ii] || oldflag[j] ;
1567 wsum += ( zr - yl ) ;
1568 ichan = j ;
1569 break ;
1570 }
1571 }
1572 else {
1573 ichan = j - 1 ;
1574 break ;
1575 }
1576 }
1577 newspec[ii] /= wsum ;
1578 wsum = 0.0 ;
1579 }
1580 }
1581 else if ( dnu < 0.0 ) {
1582 for ( int ii = 0 ; ii < nChan ; ii++ ) {
1583 double zl = zi[ii] ;
1584 double zr = zi[ii+1] ;
1585 for ( int j = ichan ; j < oldsize ; j++ ) {
1586 double yl = yi[j] ;
1587 double yr = yi[j+1] ;
1588 if ( yl >= zl ) {
1589 if ( yr >= zl ) {
1590 continue ;
1591 }
1592 else if ( yr >= zr ) {
1593 newspec[ii] += oldspec[j] * abs( yr - zl ) ;
1594 newflag[ii] = newflag[ii] || oldflag[j] ;
1595 wsum += abs( yr - zl ) ;
1596 }
1597 else {
1598 newspec[ii] += oldspec[j] * abs( dnu ) ;
1599 newflag[ii] = newflag[ii] || oldflag[j] ;
1600 wsum += abs( dnu ) ;
1601 ichan = j ;
1602 break ;
1603 }
1604 }
1605 else if ( yl > zr ) {
1606 if ( yr >= zr ) {
1607 newspec[ii] += oldspec[j] * abs( yr - yl ) ;
1608 newflag[ii] = newflag[ii] || oldflag[j] ;
1609 wsum += abs( yr - yl ) ;
1610 }
1611 else {
1612 newspec[ii] += oldspec[j] * abs( zr - yl ) ;
1613 newflag[ii] = newflag[ii] || oldflag[j] ;
1614 wsum += abs( zr - yl ) ;
1615 ichan = j ;
1616 break ;
1617 }
1618 }
1619 else {
1620 ichan = j - 1 ;
1621 break ;
1622 }
1623 }
1624 newspec[ii] /= wsum ;
1625 wsum = 0.0 ;
1626 }
1627 }
1628// * ichan = 0
1629// ***/
1630// //ofs << "olddnu = " << olddnu << ", dnu = " << dnu << endl ;
1631// pile += dnu ;
1632// wedge = olddnu * ( refChan + 1 ) ;
1633// while ( wedge < pile ) {
1634// newspec[0] += olddnu * oldspec[refChan] ;
1635// newflag[0] = newflag[0] || oldflag[refChan] ;
1636// //ofs << "channel " << refChan << " is included in new channel 0" << endl ;
1637// refChan++ ;
1638// wedge += olddnu ;
1639// wsum += olddnu ;
1640// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1641// }
1642// frac = ( wedge - pile ) / olddnu ;
1643// wsum += ( 1.0 - frac ) * olddnu ;
1644// newspec[0] += ( 1.0 - frac ) * olddnu * oldspec[refChan] ;
1645// newflag[0] = newflag[0] || oldflag[refChan] ;
1646// //ofs << "channel " << refChan << " is partly included in new channel 0" << " with fraction of " << ( 1.0 - frac ) << endl ;
1647// //ofs << "newspec[0] = " << newspec[0] << " wsum = " << wsum << endl ;
1648// newspec[0] /= wsum ;
1649// //ofs << "newspec[0] = " << newspec[0] << endl ;
1650// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1651
1652// /***
1653// * ichan = 1 - nChan-2
1654// ***/
1655// for ( int ichan = 1 ; ichan < nChan - 1 ; ichan++ ) {
1656// pile += dnu ;
1657// newspec[ichan] += frac * olddnu * oldspec[refChan] ;
1658// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1659// //ofs << "channel " << refChan << " is partly included in new channel " << ichan << " with fraction of " << frac << endl ;
1660// refChan++ ;
1661// wedge += olddnu ;
1662// wsum = frac * olddnu ;
1663// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1664// while ( wedge < pile ) {
1665// newspec[ichan] += olddnu * oldspec[refChan] ;
1666// newflag[ichan] = newflag[ichan] || oldflag[refChan] ;
1667// //ofs << "channel " << refChan << " is included in new channel " << ichan << endl ;
1668// refChan++ ;
1669// wedge += olddnu ;
1670// wsum += olddnu ;
1671// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1672// }
1673// frac = ( wedge - pile ) / olddnu ;
1674// wsum += ( 1.0 - frac ) * olddnu ;
1675// newspec[ichan] += ( 1.0 - 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 " << ( 1.0 - frac ) << endl ;
1678// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1679// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << " wsum = " << wsum << endl ;
1680// newspec[ichan] /= wsum ;
1681// //ofs << "newspec[" << ichan << "] = " << newspec[ichan] << endl ;
1682// }
1683
1684// /***
1685// * ichan = nChan-1
1686// ***/
1687// // NOTE: Assumed that all spectra have the same bandwidth
1688// pile += dnu ;
1689// newspec[nChan-1] += frac * olddnu * oldspec[refChan] ;
1690// newflag[nChan-1] = newflag[nChan-1] || oldflag[refChan] ;
1691// //ofs << "channel " << refChan << " is partly included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1692// refChan++ ;
1693// wedge += olddnu ;
1694// wsum = frac * olddnu ;
1695// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1696// for ( int jchan = refChan ; jchan < oldsize ; jchan++ ) {
1697// newspec[nChan-1] += olddnu * oldspec[jchan] ;
1698// newflag[nChan-1] = newflag[nChan-1] || oldflag[jchan] ;
1699// wsum += olddnu ;
1700// //ofs << "channel " << jchan << " is included in new channel " << nChan-1 << " with fraction of " << frac << endl ;
1701// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1702// }
1703// //ofs << "wedge = " << wedge << ", pile = " << pile << endl ;
1704// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << " wsum = " << wsum << endl ;
1705// newspec[nChan-1] /= wsum ;
1706// //ofs << "newspec[" << nChan - 1 << "] = " << newspec[nChan-1] << endl ;
1707
1708// specCol_.put( irow, newspec ) ;
1709// flagsCol_.put( irow, newflag ) ;
1710
1711// // ofs.close() ;
1712
1713
1714 return ;
1715}
1716
[1730]1717std::vector<float> Scantable::getWeather(int whichrow) const
1718{
1719 std::vector<float> out(5);
1720 //Float temperature, pressure, humidity, windspeed, windaz;
1721 weatherTable_.getEntry(out[0], out[1], out[2], out[3], out[4],
1722 mweatheridCol_(uInt(whichrow)));
1723
1724
1725 return out;
[1391]1726}
[1730]1727
[1907]1728bool Scantable::getFlagtraFast(int whichrow)
1729{
1730 uChar flag;
1731 Vector<uChar> flags;
1732 flagsCol_.get(uInt(whichrow), flags);
1733 for (int i = 0; i < flags.size(); i++) {
1734 if (i==0) {
1735 flag = flags[i];
1736 }
1737 else {
1738 flag &= flags[i];
1739 }
1740 return ((flag >> 7) == 1);
1741 }
[1730]1742}
[1907]1743
1744void Scantable::doPolyBaseline(const std::vector<bool>& mask, int order, int rowno, Fitter& fitter)
1745{
1746 fitter.setExpression("poly", order);
1747
1748 std::vector<double> abcsd = getAbcissa(rowno);
1749 std::vector<float> abcs;
1750 for (int i = 0; i < abcsd.size(); i++) {
1751 abcs.push_back((float)abcsd[i]);
1752 }
1753 std::vector<float> spec = getSpectrum(rowno);
1754 std::vector<bool> fmask = getMask(rowno);
[1908]1755 if (fmask.size() != mask.size()) {
1756 throw(AipsError("different mask sizes"));
1757 }
[1907]1758 for (int i = 0; i < fmask.size(); i++) {
1759 fmask[i] = fmask[i] && mask[i];
1760 }
1761 fitter.setData(abcs, spec, fmask);
1762
1763 fitter.lfit();
1764}
1765
[1931]1766void Scantable::polyBaselineBatch(const std::vector<bool>& mask, int order)
[1907]1767{
1768 Fitter fitter = Fitter();
[1931]1769 for (uInt rowno=0; rowno < nrow(); ++rowno) {
1770 doPolyBaseline(mask, order, rowno, fitter);
1771 setSpectrum(fitter.getResidual(), rowno);
1772 }
[1907]1773}
1774
[1931]1775STFitEntry Scantable::polyBaseline(const std::vector<bool>& mask, int order, int rowno)
[1907]1776{
1777 Fitter fitter = Fitter();
1778 doPolyBaseline(mask, order, rowno, fitter);
1779 setSpectrum(fitter.getResidual(), rowno);
[1931]1780 return fitter.getFitEntry();
[1907]1781}
1782
1783}
[1819]1784//namespace asap
Note: See TracBrowser for help on using the repository browser.