source: trunk/src/Scantable.cpp@ 912

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

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