source: trunk/src/Scantable.cpp@ 938

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

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