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

Last change on this file since 1662 was 1656, checked in by WataruKawasaki, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-1433)

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: functionalities for row-based flagging added

Test Programs:

Put in Release Notes: No

Module(s): sdflag

Description: Added the following functionalities:

sd.scantable._flag_row()

  • to execute row-based flagging

sd.scantable._getflagrow()

  • to get row-based flagging info


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