source: branches/mergetest/src/Scantable.cpp@ 1859

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

New Development: No

JIRA Issue: No (merge alma branch to trunk)

Ready for Test: Yes

Interface Changes: No

Test Programs: regressions may work

Module(s): all single dish modules

Description:

Merged all changes in alma (r1386:1818) and newfiller (r1774:1818) branch.


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