source: branches/alma/src/Scantable.cpp@ 1793

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

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.


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