source: trunk/src/Scantable.cpp@ 913

Last change on this file since 913 was 913, checked in by mar637, 19 years ago

Stokes polarimetry handling and conversion
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.9 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
14#include <casa/aips.h>
15#include <casa/iostream.h>
16#include <casa/iomanip.h>
17#include <casa/OS/Path.h>
18#include <casa/OS/File.h>
19#include <casa/Arrays/Array.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/MaskArrMath.h>
22#include <casa/Arrays/ArrayLogical.h>
23#include <casa/Arrays/ArrayAccessor.h>
24#include <casa/Arrays/VectorSTLIterator.h>
25#include <casa/Arrays/Vector.h>
26#include <casa/BasicMath/Math.h>
27#include <casa/BasicSL/Constants.h>
28#include <casa/Quanta/MVAngle.h>
29#include <casa/Containers/RecordField.h>
30#include <casa/Utilities/GenSort.h>
31
32#include <tables/Tables/TableParse.h>
33#include <tables/Tables/TableDesc.h>
34#include <tables/Tables/TableCopy.h>
35#include <tables/Tables/SetupNewTab.h>
36#include <tables/Tables/ScaColDesc.h>
37#include <tables/Tables/ArrColDesc.h>
38#include <tables/Tables/TableRow.h>
39#include <tables/Tables/TableVector.h>
40#include <tables/Tables/TableIter.h>
41
42#include <tables/Tables/ExprNode.h>
43#include <tables/Tables/TableRecord.h>
44#include <measures/Measures/MFrequency.h>
45#include <measures/Measures/MEpoch.h>
46#include <measures/Measures/MeasTable.h>
47#include <measures/Measures/MeasRef.h>
48#include <measures/TableMeasures/TableMeasRefDesc.h>
49#include <measures/TableMeasures/TableMeasValueDesc.h>
50#include <measures/TableMeasures/TableMeasDesc.h>
51#include <measures/TableMeasures/ScalarMeasColumn.h>
52#include <coordinates/Coordinates/CoordinateUtil.h>
53#include <casa/Quanta/MVTime.h>
54#include <casa/Quanta/MVAngle.h>
55
56#include "Scantable.h"
57#include "STPolLinear.h"
58#include "STPolStokes.h"
59#include "STAttr.h"
60#include "MathUtils.h"
61
62using namespace casa;
63
64namespace asap {
65
66std::map<std::string, STPol::STPolFactory *> Scantable::factories_;
67
68void Scantable::initFactories() {
69 if ( factories_.empty() ) {
70 Scantable::factories_["linear"] = &STPolLinear::myFactory;
71 Scantable::factories_["stokes"] = &STPolStokes::myFactory;
72 }
73}
74
75Scantable::Scantable(Table::TableType ttype) :
76 type_(ttype)
77{
78 initFactories();
79 setupMainTable();
80 freqTable_ = STFrequencies(*this);
81 table_.rwKeywordSet().defineTable("FREQUENCIES", freqTable_.table());
82 weatherTable_ = STWeather(*this);
83 table_.rwKeywordSet().defineTable("WEATHER", weatherTable_.table());
84 focusTable_ = STFocus(*this);
85 table_.rwKeywordSet().defineTable("FOCUS", focusTable_.table());
86 tcalTable_ = STTcal(*this);
87 table_.rwKeywordSet().defineTable("TCAL", tcalTable_.table());
88 moleculeTable_ = STMolecules(*this);
89 table_.rwKeywordSet().defineTable("MOLECULES", moleculeTable_.table());
90 historyTable_ = STHistory(*this);
91 table_.rwKeywordSet().defineTable("HISTORY", historyTable_.table());
92 setupFitTable();
93 fitTable_ = table_.keywordSet().asTable("FITS");
94 originalTable_ = table_;
95 attach();
96}
97
98Scantable::Scantable(const std::string& name, Table::TableType ttype) :
99 type_(ttype)
100{
101 initFactories();
102 Table tab(name, Table::Update);
103 Int version;
104 tab.keywordSet().get("VERSION", version);
105 if (version != version_) {
106 throw(AipsError("Unsupported version of ASAP file."));
107 }
108 if ( type_ == Table::Memory )
109 table_ = tab.copyToMemoryTable(generateName());
110 else
111 table_ = tab;
112 attachSubtables();
113 originalTable_ = table_;
114 attach();
115}
116
117Scantable::Scantable( const Scantable& other, bool clear )
118{
119 // with or without data
120 String newname = String(generateName());
121 type_ = other.table_.tableType();
122 if ( other.table_.tableType() == Table::Memory ) {
123 if ( clear ) {
124 table_ = TableCopy::makeEmptyMemoryTable(newname,
125 other.table_, True);
126 } else
127 table_ = other.table_.copyToMemoryTable(newname);
128 } else {
129 other.table_.deepCopy(newname, Table::New, False, Table::AipsrcEndian,
130 Bool(clear));
131 table_ = Table(newname, Table::Update);
132 if ( clear ) copySubtables(other);
133 table_.markForDelete();
134 }
135
136 attachSubtables();
137 originalTable_ = table_;
138 attach();
139}
140
141void Scantable::copySubtables(const Scantable& other) {
142 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
143 TableCopy::copyRows(t, other.freqTable_.table());
144 t = table_.rwKeywordSet().asTable("FOCUS");
145 TableCopy::copyRows(t, other.focusTable_.table());
146 t = table_.rwKeywordSet().asTable("WEATHER");
147 TableCopy::copyRows(t, other.weatherTable_.table());
148 t = table_.rwKeywordSet().asTable("TCAL");
149 TableCopy::copyRows(t, other.tcalTable_.table());
150 t = table_.rwKeywordSet().asTable("MOLECULES");
151 TableCopy::copyRows(t, other.moleculeTable_.table());
152 t = table_.rwKeywordSet().asTable("HISTORY");
153 TableCopy::copyRows(t, other.historyTable_.table());
154}
155
156void Scantable::attachSubtables()
157{
158 freqTable_ = STFrequencies(table_);
159 focusTable_ = STFocus(table_);
160 weatherTable_ = STWeather(table_);
161 tcalTable_ = STTcal(table_);
162 moleculeTable_ = STMolecules(table_);
163 historyTable_ = STHistory(table_);
164}
165
166Scantable::~Scantable()
167{
168 cout << "~Scantable() " << this << endl;
169}
170
171void Scantable::setupMainTable()
172{
173 TableDesc td("", "1", TableDesc::Scratch);
174 td.comment() = "An ASAP Scantable";
175 td.rwKeywordSet().define("VERSION", Int(version_));
176
177 // n Cycles
178 td.addColumn(ScalarColumnDesc<uInt>("SCANNO"));
179 // new index every nBeam x nIF x nPol
180 td.addColumn(ScalarColumnDesc<uInt>("CYCLENO"));
181
182 td.addColumn(ScalarColumnDesc<uInt>("BEAMNO"));
183 td.addColumn(ScalarColumnDesc<uInt>("IFNO"));
184 td.rwKeywordSet().define("POLTYPE", String("linear"));
185 td.addColumn(ScalarColumnDesc<uInt>("POLNO"));
186
187 td.addColumn(ScalarColumnDesc<uInt>("FREQ_ID"));
188 td.addColumn(ScalarColumnDesc<uInt>("MOLECULE_ID"));
189 // linear, circular, stokes [I Q U V], stokes1 [I Plinear Pangle V]
190 td.addColumn(ScalarColumnDesc<Int>("REFBEAMNO"));
191
192 td.addColumn(ScalarColumnDesc<Double>("TIME"));
193 TableMeasRefDesc measRef(MEpoch::UTC); // UTC as default
194 TableMeasValueDesc measVal(td, "TIME");
195 TableMeasDesc<MEpoch> mepochCol(measVal, measRef);
196 mepochCol.write(td);
197
198 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
199
200 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
201 // Type of source (on=0, off=1, other=-1)
202 td.addColumn(ScalarColumnDesc<Int>("SRCTYPE", Int(-1)));
203 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
204
205 //The actual Data Vectors
206 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
207 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
208 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
209
210 td.addColumn(ArrayColumnDesc<Double>("DIRECTION",
211 IPosition(1,2),
212 ColumnDesc::Direct));
213 TableMeasRefDesc mdirRef(MDirection::J2000); // default
214 TableMeasValueDesc tmvdMDir(td, "DIRECTION");
215 // the TableMeasDesc gives the column a type
216 TableMeasDesc<MDirection> mdirCol(tmvdMDir, mdirRef);
217 // writing create the measure column
218 mdirCol.write(td);
219 td.addColumn(ScalarColumnDesc<Double>("AZIMUTH"));
220 td.addColumn(ScalarColumnDesc<Double>("ELEVATION"));
221 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
222
223 td.addColumn(ScalarColumnDesc<uInt>("TCAL_ID"));
224 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID"));
225
226 td.addColumn(ScalarColumnDesc<uInt>("FOCUS_ID"));
227 td.addColumn(ScalarColumnDesc<uInt>("WEATHER_ID"));
228
229 td.rwKeywordSet().define("OBSMODE", String(""));
230
231 // Now create Table SetUp from the description.
232 SetupNewTable aNewTab(generateName(), td, Table::Scratch);
233 table_ = Table(aNewTab, type_, 0);
234 originalTable_ = table_;
235
236}
237
238void Scantable::setupFitTable()
239{
240 TableDesc td("", "1", TableDesc::Scratch);
241 td.addColumn(ScalarColumnDesc<uInt>("FIT_ID"));
242 td.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
243 td.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
244 td.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
245 td.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
246 td.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
247 SetupNewTable aNewTab("fits", td, Table::Scratch);
248 Table aTable(aNewTab, Table::Memory);
249 table_.rwKeywordSet().defineTable("FITS", aTable);
250}
251
252void Scantable::attach()
253{
254 timeCol_.attach(table_, "TIME");
255 srcnCol_.attach(table_, "SRCNAME");
256 specCol_.attach(table_, "SPECTRA");
257 flagsCol_.attach(table_, "FLAGTRA");
258 tsysCol_.attach(table_, "TSYS");
259 cycleCol_.attach(table_,"CYCLENO");
260 scanCol_.attach(table_, "SCANNO");
261 beamCol_.attach(table_, "BEAMNO");
262 ifCol_.attach(table_, "IFNO");
263 polCol_.attach(table_, "POLNO");
264 integrCol_.attach(table_, "INTERVAL");
265 azCol_.attach(table_, "AZIMUTH");
266 elCol_.attach(table_, "ELEVATION");
267 dirCol_.attach(table_, "DIRECTION");
268 paraCol_.attach(table_, "PARANGLE");
269 fldnCol_.attach(table_, "FIELDNAME");
270 rbeamCol_.attach(table_, "REFBEAMNO");
271
272 mfitidCol_.attach(table_,"FIT_ID");
273 //fitidCol_.attach(fitTable_,"FIT_ID");
274
275 mfreqidCol_.attach(table_, "FREQ_ID");
276
277 mtcalidCol_.attach(table_, "TCAL_ID");
278
279 mfocusidCol_.attach(table_, "FOCUS_ID");
280
281 mmolidCol_.attach(table_, "MOLECULE_ID");
282}
283
284void Scantable::setHeader(const STHeader& sdh)
285{
286 table_.rwKeywordSet().define("nIF", sdh.nif);
287 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
288 table_.rwKeywordSet().define("nPol", sdh.npol);
289 table_.rwKeywordSet().define("nChan", sdh.nchan);
290 table_.rwKeywordSet().define("Observer", sdh.observer);
291 table_.rwKeywordSet().define("Project", sdh.project);
292 table_.rwKeywordSet().define("Obstype", sdh.obstype);
293 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
294 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
295 table_.rwKeywordSet().define("Equinox", sdh.equinox);
296 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
297 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
298 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
299 table_.rwKeywordSet().define("UTC", sdh.utc);
300 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
301 table_.rwKeywordSet().define("Epoch", sdh.epoch);
302 table_.rwKeywordSet().define("POLTYPE", sdh.poltype);
303}
304
305STHeader Scantable::getHeader() const
306{
307 STHeader sdh;
308 table_.keywordSet().get("nBeam",sdh.nbeam);
309 table_.keywordSet().get("nIF",sdh.nif);
310 table_.keywordSet().get("nPol",sdh.npol);
311 table_.keywordSet().get("nChan",sdh.nchan);
312 table_.keywordSet().get("Observer", sdh.observer);
313 table_.keywordSet().get("Project", sdh.project);
314 table_.keywordSet().get("Obstype", sdh.obstype);
315 table_.keywordSet().get("AntennaName", sdh.antennaname);
316 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
317 table_.keywordSet().get("Equinox", sdh.equinox);
318 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
319 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
320 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
321 table_.keywordSet().get("UTC", sdh.utc);
322 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
323 table_.keywordSet().get("Epoch", sdh.epoch);
324 table_.keywordSet().get("POLTYPE", sdh.poltype);
325 return sdh;
326}
327
328bool Scantable::conformant( const Scantable& other )
329{
330 return this->getHeader().conformant(other.getHeader());
331}
332
333
334int Scantable::nscan() const {
335 int n = 0;
336 Int previous = -1; Int current = 0;
337 Vector<uInt> scannos(scanCol_.getColumn());
338 uInt nout = GenSort<uInt>::sort( scannos, Sort::Ascending,
339 Sort::QuickSort|Sort::NoDuplicates );
340 return int(nout);
341}
342
343std::string Scantable::formatSec(Double x) const
344{
345 Double xcop = x;
346 MVTime mvt(xcop/24./3600.); // make days
347
348 if (x < 59.95)
349 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
350 else if (x < 3599.95)
351 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
352 else {
353 ostringstream oss;
354 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
355 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
356 return String(oss);
357 }
358};
359
360std::string Scantable::formatDirection(const MDirection& md) const
361{
362 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
363 Int prec = 7;
364
365 MVAngle mvLon(t[0]);
366 String sLon = mvLon.string(MVAngle::TIME,prec);
367 MVAngle mvLat(t[1]);
368 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
369 return sLon + String(" ") + sLat;
370}
371
372
373std::string Scantable::getFluxUnit() const
374{
375 return table_.keywordSet().asString("FluxUnit");
376}
377
378void Scantable::setFluxUnit(const std::string& unit)
379{
380 String tmp(unit);
381 Unit tU(tmp);
382 if (tU==Unit("K") || tU==Unit("Jy")) {
383 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
384 } else {
385 throw AipsError("Illegal unit - must be compatible with Jy or K");
386 }
387}
388
389void Scantable::setInstrument(const std::string& name)
390{
391 bool throwIt = true;
392 Instrument ins = STAttr::convertInstrument(name, throwIt);
393 String nameU(name);
394 nameU.upcase();
395 table_.rwKeywordSet().define(String("AntennaName"), nameU);
396}
397
398MPosition Scantable::getAntennaPosition () const
399{
400 Vector<Double> antpos;
401 table_.keywordSet().get("AntennaPosition", antpos);
402 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
403 return MPosition(mvpos);
404}
405
406void Scantable::makePersistent(const std::string& filename)
407{
408 String inname(filename);
409 Path path(inname);
410 inname = path.expandedName();
411 table_.deepCopy(inname, Table::New);
412}
413
414int Scantable::nbeam( int scanno ) const
415{
416 if ( scanno < 0 ) {
417 Int n;
418 table_.keywordSet().get("nBeam",n);
419 return int(n);
420 } else {
421 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
422 Table t = table_(table_.col("SCANNO") == scanno);
423 ROTableRow row(t);
424 const TableRecord& rec = row.get(0);
425 Table subt = t( t.col("IFNO") == Int(rec.asuInt("IFNO"))
426 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
427 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
428 ROTableVector<uInt> v(subt, "BEAMNO");
429 return int(v.nelements());
430 }
431 return 0;
432}
433
434int Scantable::nif( int scanno ) const
435{
436 if ( scanno < 0 ) {
437 Int n;
438 table_.keywordSet().get("nIF",n);
439 return int(n);
440 } else {
441 // take the first POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
442 Table t = table_(table_.col("SCANNO") == scanno);
443 ROTableRow row(t);
444 const TableRecord& rec = row.get(0);
445 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
446 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
447 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
448 if ( subt.nrow() == 0 ) return 0;
449 ROTableVector<uInt> v(subt, "IFNO");
450 return int(v.nelements());
451 }
452 return 0;
453}
454
455int Scantable::npol( int scanno ) const
456{
457 if ( scanno < 0 ) {
458 Int n;
459 table_.keywordSet().get("nPol",n);
460 return n;
461 } else {
462 // take the first POLNO,IFNO,CYCLENO as nbeam shouldn't vary with these
463 Table t = table_(table_.col("SCANNO") == scanno);
464 ROTableRow row(t);
465 const TableRecord& rec = row.get(0);
466 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
467 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
468 && t.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
469 if ( subt.nrow() == 0 ) return 0;
470 ROTableVector<uInt> v(subt, "POLNO");
471 return int(v.nelements());
472 }
473 return 0;
474}
475
476int Scantable::ncycle( int scanno ) const
477{
478 if ( scanno < 0 ) {
479 Block<String> cols(2);
480 cols[0] = "SCANNO";
481 cols[1] = "CYCLENO";
482 TableIterator it(table_, cols);
483 int n = 0;
484 while ( !it.pastEnd() ) {
485 ++n;
486 ++it;
487 }
488 return n;
489 } else {
490 Table t = table_(table_.col("SCANNO") == scanno);
491 ROTableRow row(t);
492 const TableRecord& rec = row.get(0);
493 Table subt = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
494 && t.col("POLNO") == Int(rec.asuInt("POLNO"))
495 && t.col("IFNO") == Int(rec.asuInt("IFNO")) );
496 if ( subt.nrow() == 0 ) return 0;
497 return int(subt.nrow());
498 }
499 return 0;
500}
501
502
503int Scantable::nrow( int scanno ) const
504{
505 return int(table_.nrow());
506}
507
508int Scantable::nchan( int ifno ) const
509{
510 if ( ifno < 0 ) {
511 Int n;
512 table_.keywordSet().get("nChan",n);
513 return int(n);
514 } else {
515 // take the first SCANNO,POLNO,BEAMNO,CYCLENO as nbeam shouldn't vary with these
516 Table t = table_(table_.col("IFNO") == ifno);
517 if ( t.nrow() == 0 ) return 0;
518 ROArrayColumn<Float> v(t, "SPECTRA");
519 return v(0).nelements();
520 }
521 return 0;
522}
523
524
525int Scantable::getBeam(int whichrow) const
526{
527 return beamCol_(whichrow);
528}
529
530int Scantable::getIF(int whichrow) const
531{
532 return ifCol_(whichrow);
533}
534
535int Scantable::getPol(int whichrow) const
536{
537 return polCol_(whichrow);
538}
539
540std::string Scantable::formatTime(const MEpoch& me, bool showdate) const
541{
542 MVTime mvt(me.getValue());
543 if (showdate)
544 mvt.setFormat(MVTime::YMD);
545 else
546 mvt.setFormat(MVTime::TIME);
547 ostringstream oss;
548 oss << mvt;
549 return String(oss);
550}
551
552void Scantable::calculateAZEL()
553{
554 MPosition mp = getAntennaPosition();
555 MEpoch::ROScalarColumn timeCol(table_, "TIME");
556 ostringstream oss;
557 oss << "Computed azimuth/elevation using " << endl
558 << mp << endl;
559 for (uInt i=0; i<nrow(); ++i) {
560 MEpoch me = timeCol(i);
561 MDirection md = dirCol_(i);
562 dirCol_.get(i,md);
563 oss << " Time: " << formatTime(me,False) << " Direction: " << formatDirection(md)
564 << endl << " => ";
565 MeasFrame frame(mp, me);
566 Vector<Double> azel =
567 MDirection::Convert(md, MDirection::Ref(MDirection::AZEL,
568 frame)
569 )().getAngle("rad").getValue();
570 azCol_.put(i,azel[0]);
571 elCol_.put(i,azel[1]);
572 oss << "azel: " << azel[0]/C::pi*180.0 << " "
573 << azel[1]/C::pi*180.0 << " (deg)" << endl;
574 }
575 pushLog(String(oss));
576}
577
578void Scantable::flag()
579{
580 if ( selector_.empty() )
581 throw(AipsError("Trying to flag whole scantable. Aborted."));
582 TableVector<uChar> tvec(table_, "FLAGTRA");
583 uChar userflag = 1 << 7;
584 tvec = userflag;
585}
586
587std::vector<bool> Scantable::getMask(int whichrow) const
588{
589 Vector<uChar> flags;
590 flagsCol_.get(uInt(whichrow), flags);
591 Vector<Bool> bflag(flags.shape());
592 convertArray(bflag, flags);
593 bflag = !bflag;
594 std::vector<bool> mask;
595 bflag.tovector(mask);
596 return mask;
597}
598
599std::vector<float> Scantable::getSpectrum( int whichrow,
600 const std::string& poltype ) const
601{
602 String ptype = poltype;
603 if (poltype == "" ) ptype = getPolType();
604 if ( whichrow < 0 || whichrow >= nrow() )
605 throw(AipsError("Illegal row number."));
606 std::vector<float> out;
607 Vector<Float> arr;
608 uInt requestedpol = polCol_(whichrow);
609 String basetype = getPolType();
610 if ( ptype == basetype ) {
611 specCol_.get(whichrow, arr);
612 } else {
613 STPol* stpol = 0;
614 stpol =STPol::getPolClass(Scantable::factories_, basetype);
615 try {
616 uInt row = uInt(whichrow);
617 stpol->setSpectra(getPolMatrix(row));
618 Float frot,fang,ftan;
619 focusTable_.getEntry(frot, fang, ftan, mfocusidCol_(row));
620 stpol->setPhaseCorrections(frot, fang, ftan);
621 arr = stpol->getSpectrum(requestedpol, ptype);
622 delete stpol;
623 } catch (AipsError& e) {
624 delete stpol;
625 throw(e);
626 }
627 }
628 if ( arr.nelements() == 0 )
629 pushLog("Not enough polarisations present to do the conversion.");
630 arr.tovector(out);
631 return out;
632}
633
634void asap::Scantable::setSpectrum( const std::vector<float>& spec,
635 int whichrow )
636{
637 Vector<Float> spectrum(spec);
638 Vector<Float> arr;
639 specCol_.get(whichrow, arr);
640 if ( spectrum.nelements() != arr.nelements() )
641 throw AipsError("The spectrum has incorrect number of channels.");
642 specCol_.put(whichrow, spectrum);
643}
644
645
646String Scantable::generateName()
647{
648 return (File::newUniqueName("./","temp")).baseName();
649}
650
651const casa::Table& Scantable::table( ) const
652{
653 return table_;
654}
655
656casa::Table& Scantable::table( )
657{
658 return table_;
659}
660
661std::string Scantable::getPolType() const
662{
663 return table_.keywordSet().asString("POLTYPE");
664}
665
666void Scantable::unsetSelection()
667{
668 table_ = originalTable_;
669 attach();
670 selector_.reset();
671}
672
673void Scantable::setSelection( const STSelector& selection )
674{
675 Table tab = const_cast<STSelector&>(selection).apply(originalTable_);
676 if ( tab.nrow() == 0 ) {
677 throw(AipsError("Selection contains no data. Not applying it."));
678 }
679 table_ = tab;
680 attach();
681 selector_ = selection;
682}
683
684std::string Scantable::summary( bool verbose )
685{
686 // Format header info
687 ostringstream oss;
688 oss << endl;
689 oss << asap::SEPERATOR << endl;
690 oss << " Scan Table Summary" << endl;
691 oss << asap::SEPERATOR << endl;
692 oss.flags(std::ios_base::left);
693 oss << setw(15) << "Beams:" << setw(4) << nbeam() << endl
694 << setw(15) << "IFs:" << setw(4) << nif() << endl
695 << setw(15) << "Polarisations:" << setw(4) << npol()
696 << "(" << getPolType() << ")" << endl
697 << setw(15) << "Channels:" << setw(4) << nchan() << endl;
698 oss << endl;
699 String tmp;
700 oss << setw(15) << "Observer:"
701 << table_.keywordSet().asString("Observer") << endl;
702 oss << setw(15) << "Obs Date:" << getTime(-1,true) << endl;
703 table_.keywordSet().get("Project", tmp);
704 oss << setw(15) << "Project:" << tmp << endl;
705 table_.keywordSet().get("Obstype", tmp);
706 oss << setw(15) << "Obs. Type:" << tmp << endl;
707 table_.keywordSet().get("AntennaName", tmp);
708 oss << setw(15) << "Antenna Name:" << tmp << endl;
709 table_.keywordSet().get("FluxUnit", tmp);
710 oss << setw(15) << "Flux Unit:" << tmp << endl;
711 Vector<Float> vec;
712 oss << setw(15) << "Rest Freqs:";
713 if (vec.nelements() > 0) {
714 oss << setprecision(10) << vec << " [Hz]" << endl;
715 } else {
716 oss << "none" << endl;
717 }
718 oss << setw(15) << "Abcissa:" << "channel" << endl;
719 oss << selector_.print() << endl;
720 oss << endl;
721 // main table
722 String dirtype = "Position ("
723 + MDirection::showType(dirCol_.getMeasRef().getType())
724 + ")";
725 oss << setw(5) << "Scan"
726 << setw(15) << "Source"
727// << setw(24) << dirtype
728 << setw(10) << "Time"
729 << setw(18) << "Integration" << endl
730 << setw(5) << "" << setw(10) << "Beam" << dirtype << endl
731 << setw(15) << "" << setw(5) << "IF"
732 << setw(8) << "Frame" << setw(16)
733 << "RefVal" << setw(10) << "RefPix" << setw(12) << "Increment" <<endl;
734 oss << asap::SEPERATOR << endl;
735 TableIterator iter(table_, "SCANNO");
736 while (!iter.pastEnd()) {
737 Table subt = iter.table();
738 ROTableRow row(subt);
739 MEpoch::ROScalarColumn timeCol(subt,"TIME");
740 const TableRecord& rec = row.get(0);
741 oss << setw(4) << std::right << rec.asuInt("SCANNO")
742 << std::left << setw(1) << ""
743 << setw(15) << rec.asString("SRCNAME")
744 << setw(10) << formatTime(timeCol(0), false);
745 // count the cycles in the scan
746 TableIterator cyciter(subt, "CYCLENO");
747 int nint = 0;
748 while (!cyciter.pastEnd()) {
749 ++nint;
750 ++cyciter;
751 }
752 oss << setw(3) << std::right << nint << setw(3) << " x " << std::left
753 << setw(6) << formatSec(rec.asFloat("INTERVAL")) << endl;
754
755 TableIterator biter(subt, "BEAMNO");
756 while (!biter.pastEnd()) {
757 Table bsubt = biter.table();
758 ROTableRow brow(bsubt);
759 MDirection::ROScalarColumn bdirCol(bsubt,"DIRECTION");
760 const TableRecord& brec = brow.get(0);
761 oss << setw(6) << "" << setw(10) << brec.asuInt("BEAMNO");
762 oss << setw(24) << formatDirection(bdirCol(0)) << endl;
763 TableIterator iiter(bsubt, "IFNO");
764 while (!iiter.pastEnd()) {
765 Table isubt = iiter.table();
766 ROTableRow irow(isubt);
767 const TableRecord& irec = irow.get(0);
768 oss << std::right <<setw(8) << "" << std::left << irec.asuInt("IFNO");
769 oss << frequencies().print(irec.asuInt("FREQ_ID"));
770
771 ++iiter;
772 }
773 ++biter;
774 }
775 ++iter;
776 }
777 /// @todo implement verbose mode
778 return String(oss);
779}
780
781std::string Scantable::getTime(int whichrow, bool showdate) const
782{
783 MEpoch::ROScalarColumn timeCol(table_, "TIME");
784 MEpoch me;
785 if (whichrow > -1) {
786 me = timeCol(uInt(whichrow));
787 } else {
788 Double tm;
789 table_.keywordSet().get("UTC",tm);
790 me = MEpoch(MVEpoch(tm));
791 }
792 return formatTime(me, showdate);
793}
794
795std::vector< double > asap::Scantable::getAbcissa( int whichrow ) const
796{
797 if ( whichrow > table_.nrow() ) throw(AipsError("Illegal ro number"));
798 std::vector<double> stlout;
799 int nchan = specCol_(whichrow).nelements();
800 String us = freqTable_.getUnitString();
801 if ( us == "" || us == "pixel" || us == "channel" ) {
802 for (int i=0; i<nchan; ++i) {
803 stlout.push_back(double(i));
804 }
805 return stlout;
806 }
807
808 const MPosition& mp = getAntennaPosition();
809 const MDirection& md = dirCol_(whichrow);
810 const MEpoch& me = timeCol_(whichrow);
811 Double rf = moleculeTable_.getRestFrequency(mmolidCol_(whichrow));
812 SpectralCoordinate spc =
813 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
814 Vector<Double> pixel(nchan);
815 Vector<Double> world;
816 indgen(pixel);
817 if ( Unit(us) == Unit("Hz") ) {
818 for ( int i=0; i < nchan; ++i) {
819 Double world;
820 spc.toWorld(world, pixel[i]);
821 stlout.push_back(double(world));
822 }
823 } else if ( Unit(us) == Unit("km/s") ) {
824 Vector<Double> world;
825 spc.pixelToVelocity(world, pixel);
826 world.tovector(stlout);
827 }
828 return stlout;
829}
830
831std::string Scantable::getAbcissaLabel( int whichrow ) const
832{
833 if ( whichrow > table_.nrow() ) throw(AipsError("Illegal ro number"));
834 const MPosition& mp = getAntennaPosition();
835 const MDirection& md = dirCol_(whichrow);
836 const MEpoch& me = timeCol_(whichrow);
837 const Double& rf = mmolidCol_(whichrow);
838 SpectralCoordinate spc =
839 freqTable_.getSpectralCoordinate(md, mp, me, rf, mfreqidCol_(whichrow));
840
841 String s = "Channel";
842 Unit u = Unit(freqTable_.getUnitString());
843 if (u == Unit("km/s")) {
844 s = CoordinateUtil::axisLabel(spc,0,True,True,True);
845 } else if (u == Unit("Hz")) {
846 Vector<String> wau(1);wau = u.getName();
847 spc.setWorldAxisUnits(wau);
848
849 s = CoordinateUtil::axisLabel(spc,0,True,True,False);
850 }
851 return s;
852
853}
854
855void asap::Scantable::setRestFrequencies( double rf, const std::string& unit )
856{
857 ///@todo lookup in line table
858 Unit u(unit);
859 Quantum<Double> urf(rf, u);
860 uInt id = moleculeTable_.addEntry(urf.getValue("Hz"), "", "");
861 TableVector<uInt> tabvec(table_, "MOLECULE_ID");
862 tabvec = id;
863}
864
865void asap::Scantable::setRestFrequencies( const std::string& name )
866{
867 throw(AipsError("setRestFrequencies( const std::string& name ) NYI"));
868 ///@todo implement
869}
870
871std::vector< unsigned int > asap::Scantable::rownumbers( ) const
872{
873 std::vector<unsigned int> stlout;
874 Vector<uInt> vec = table_.rowNumbers();
875 vec.tovector(stlout);
876 return stlout;
877}
878
879
880Matrix<Float> asap::Scantable::getPolMatrix( uInt whichrow ) const
881{
882 ROTableRow row(table_);
883 const TableRecord& rec = row.get(whichrow);
884 Table t =
885 originalTable_( originalTable_.col("SCANNO") == Int(rec.asuInt("SCANNO"))
886 && originalTable_.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
887 && originalTable_.col("IFNO") == Int(rec.asuInt("IFNO"))
888 && originalTable_.col("CYCLENO") == Int(rec.asuInt("CYCLENO")) );
889 ROArrayColumn<Float> speccol(t, "SPECTRA");
890 return speccol.getColumn();
891}
892
893std::vector< std::string > asap::Scantable::columnNames( ) const
894{
895 Vector<String> vec = table_.tableDesc().columnNames();
896 return mathutil::tovectorstring(vec);
897}
898
899} //namespace asap
Note: See TracBrowser for help on using the repository browser.