source: trunk/src/SDMemTable.cc@ 402

Last change on this file since 402 was 401, checked in by kil064, 20 years ago

add rest freq specification by line name
add function spectraLines to list known lines
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.5 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMemTable.cc: A MemoryTable container for single dish integrations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31
32#include <map>
33
34#include <casa/aips.h>
35#include <casa/iostream.h>
36#include <casa/iomanip.h>
37#include <casa/Arrays/Array.h>
38#include <casa/Arrays/ArrayMath.h>
39#include <casa/Arrays/MaskArrMath.h>
40#include <casa/Arrays/ArrayLogical.h>
41#include <casa/Arrays/ArrayAccessor.h>
42#include <casa/Arrays/Vector.h>
43#include <casa/Quanta/MVAngle.h>
44
45#include <tables/Tables/TableParse.h>
46#include <tables/Tables/TableDesc.h>
47#include <tables/Tables/SetupNewTab.h>
48#include <tables/Tables/ScaColDesc.h>
49#include <tables/Tables/ArrColDesc.h>
50
51#include <tables/Tables/ExprNode.h>
52#include <tables/Tables/TableRecord.h>
53#include <measures/Measures/MFrequency.h>
54#include <measures/Measures/MeasTable.h>
55#include <coordinates/Coordinates/CoordinateUtil.h>
56#include <casa/Quanta/MVTime.h>
57#include <casa/Quanta/MVAngle.h>
58
59#include "SDDefs.h"
60#include "SDMemTable.h"
61#include "SDContainer.h"
62#include "MathUtils.h"
63
64
65using namespace casa;
66using namespace asap;
67
68SDMemTable::SDMemTable() :
69 IFSel_(0),
70 beamSel_(0),
71 polSel_(0)
72{
73 setup();
74 attach();
75}
76
77SDMemTable::SDMemTable(const std::string& name) :
78 IFSel_(0),
79 beamSel_(0),
80 polSel_(0)
81{
82 Table tab(name);
83 table_ = tab.copyToMemoryTable("dummy");
84 //cerr << "hello from C SDMemTable @ " << this << endl;
85 attach();
86}
87
88SDMemTable::SDMemTable(const SDMemTable& other, Bool clear)
89{
90 IFSel_= other.IFSel_;
91 beamSel_= other.beamSel_;
92 polSel_= other.polSel_;
93 chanMask_ = other.chanMask_;
94 table_ = other.table_.copyToMemoryTable(String("dummy"));
95 // clear all rows()
96 if (clear) {
97 table_.removeRow(this->table_.rowNumbers());
98 } else {
99 IFSel_ = other.IFSel_;
100 beamSel_ = other.beamSel_;
101 polSel_ = other.polSel_;
102 }
103//
104 attach();
105 //cerr << "hello from CC SDMemTable @ " << this << endl;
106}
107
108SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
109 IFSel_(0),
110 beamSel_(0),
111 polSel_(0)
112{
113 cout << exprs << endl;
114 Table t = tableCommand(exprs,tab);
115 if (t.nrow() == 0)
116 throw(AipsError("Query unsuccessful."));
117 table_ = t.copyToMemoryTable("dummy");
118 attach();
119 renumber();
120}
121
122SDMemTable::~SDMemTable()
123{
124 //cerr << "goodbye from SDMemTable @ " << this << endl;
125}
126
127SDMemTable SDMemTable::getScan(Int scanID) const
128{
129 String cond("SELECT * from $1 WHERE SCANID == ");
130 cond += String::toString(scanID);
131 return SDMemTable(table_, cond);
132}
133
134SDMemTable &SDMemTable::operator=(const SDMemTable& other)
135{
136 if (this != &other) {
137 IFSel_= other.IFSel_;
138 beamSel_= other.beamSel_;
139 polSel_= other.polSel_;
140 chanMask_.resize(0);
141 chanMask_ = other.chanMask_;
142 table_ = other.table_.copyToMemoryTable(String("dummy"));
143 attach();
144 }
145 //cerr << "hello from ASS SDMemTable @ " << this << endl;
146 return *this;
147}
148
149SDMemTable SDMemTable::getSource(const std::string& source) const
150{
151 String cond("SELECT * from $1 WHERE SRCNAME == ");
152 cond += source;
153 return SDMemTable(table_, cond);
154}
155
156void SDMemTable::setup()
157{
158 TableDesc td("", "1", TableDesc::Scratch);
159 td.comment() = "A SDMemTable";
160//
161 td.addColumn(ScalarColumnDesc<Double>("TIME"));
162 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
163 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
164 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
165 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
166 td.addColumn(ScalarColumnDesc<Int>("SCANID"));
167 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
168 td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
169 td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID"));
170 td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
171 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
172 td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
173 td.addColumn(ArrayColumnDesc<Float>("TCAL"));
174 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
175 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
176 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
177 td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
178 td.addColumn(ArrayColumnDesc<String>("HISTORY"));
179
180 // Now create a new table from the description.
181
182 SetupNewTable aNewTab("dummy", td, Table::New);
183 table_ = Table(aNewTab, Table::Memory, 0);
184}
185
186void SDMemTable::attach()
187{
188 timeCol_.attach(table_, "TIME");
189 srcnCol_.attach(table_, "SRCNAME");
190 specCol_.attach(table_, "SPECTRA");
191 flagsCol_.attach(table_, "FLAGTRA");
192 tsCol_.attach(table_, "TSYS");
193 scanCol_.attach(table_, "SCANID");
194 integrCol_.attach(table_, "INTERVAL");
195 freqidCol_.attach(table_, "FREQID");
196 restfreqidCol_.attach(table_, "RESTFREQID");
197 dirCol_.attach(table_, "DIRECTION");
198 fldnCol_.attach(table_, "FIELDNAME");
199 tcaltCol_.attach(table_, "TCALTIME");
200 tcalCol_.attach(table_, "TCAL");
201 azCol_.attach(table_, "AZIMUTH");
202 elCol_.attach(table_, "ELEVATION");
203 paraCol_.attach(table_, "PARANGLE");
204 rbeamCol_.attach(table_, "REFBEAM");
205 histCol_.attach(table_, "HISTORY");
206}
207
208
209std::string SDMemTable::getSourceName(Int whichRow) const
210{
211 String name;
212 srcnCol_.get(whichRow, name);
213 return name;
214}
215
216std::string SDMemTable::getTime(Int whichRow, Bool showDate) const
217{
218 Double tm;
219 if (whichRow > -1) {
220 timeCol_.get(whichRow, tm);
221 } else {
222 table_.keywordSet().get("UTC",tm);
223 }
224 MVTime mvt(tm);
225 if (showDate)
226 mvt.setFormat(MVTime::YMD);
227 else
228 mvt.setFormat(MVTime::TIME);
229 ostringstream oss;
230 oss << mvt;
231 return String(oss);
232}
233
234double SDMemTable::getInterval(Int whichRow) const
235{
236 Double intval;
237 integrCol_.get(whichRow, intval);
238 return intval;
239}
240
241bool SDMemTable::setIF(Int whichIF)
242{
243 if ( whichIF >= 0 && whichIF < nIF()) {
244 IFSel_ = whichIF;
245 return true;
246 }
247 return false;
248}
249
250bool SDMemTable::setBeam(Int whichBeam)
251{
252 if ( whichBeam >= 0 && whichBeam < nBeam()) {
253 beamSel_ = whichBeam;
254 return true;
255 }
256 return false;
257}
258
259bool SDMemTable::setPol(Int whichPol)
260{
261 if ( whichPol >= 0 && whichPol < nPol()) {
262 polSel_ = whichPol;
263 return true;
264 }
265 return false;
266}
267
268void SDMemTable::resetCursor()
269{
270 polSel_ = 0;
271 IFSel_ = 0;
272 beamSel_ = 0;
273}
274
275bool SDMemTable::setMask(std::vector<int> whichChans)
276{
277 std::vector<int>::iterator it;
278 uInt n = flagsCol_.shape(0)(3);
279 if (whichChans.empty()) {
280 chanMask_ = std::vector<bool>(n,true);
281 return true;
282 }
283 chanMask_.resize(n,true);
284 for (it = whichChans.begin(); it != whichChans.end(); ++it) {
285 if (*it < n) {
286 chanMask_[*it] = false;
287 }
288 }
289 return true;
290}
291
292std::vector<bool> SDMemTable::getMask(Int whichRow) const {
293 std::vector<bool> mask;
294 Array<uChar> arr;
295 flagsCol_.get(whichRow, arr);
296 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
297 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
298 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
299 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
300 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
301 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
302
303 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
304
305 std::vector<bool> tmp;
306 tmp = chanMask_; // WHY the fxxx do I have to make a copy here
307 std::vector<bool>::iterator miter;
308 miter = tmp.begin();
309
310 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
311 bool out =!static_cast<bool>(*i);
312 if (useUserMask) {
313 out = out && (*miter);
314 miter++;
315 }
316 mask.push_back(out);
317 }
318 return mask;
319}
320std::vector<float> SDMemTable::getSpectrum(Int whichRow) const
321{
322 std::vector<float> spectrum;
323 Array<Float> arr;
324 specCol_.get(whichRow, arr);
325 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
326 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
327 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
328 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
329 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
330 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
331 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
332 spectrum.push_back(*i);
333 }
334 return spectrum;
335}
336std::vector<string> SDMemTable::getCoordInfo() const
337{
338 String un;
339 Table t = table_.keywordSet().asTable("FREQUENCIES");
340 String sunit;
341 t.keywordSet().get("UNIT",sunit);
342 String dpl;
343 t.keywordSet().get("DOPPLER",dpl);
344 if (dpl == "") dpl = "RADIO";
345 String rfrm;
346 t.keywordSet().get("REFFRAME",rfrm);
347 std::vector<string> inf;
348 inf.push_back(sunit);
349 inf.push_back(rfrm);
350 inf.push_back(dpl);
351 t.keywordSet().get("BASEREFFRAME",rfrm);
352 inf.push_back(rfrm);
353 return inf;
354}
355
356void SDMemTable::setCoordInfo(std::vector<string> theinfo)
357{
358 std::vector<string>::iterator it;
359 String un,rfrm, brfrm,dpl;
360 un = theinfo[0]; // Abcissa unit
361 rfrm = theinfo[1]; // Active (or conversion) frame
362 dpl = theinfo[2]; // Doppler
363 brfrm = theinfo[3]; // Base frame
364//
365 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
366//
367 Vector<Double> rstf;
368 t.keywordSet().get("RESTFREQS",rstf);
369//
370 Bool canDo = True;
371 Unit u1("km/s");Unit u2("Hz");
372 if (Unit(un) == u1) {
373 Vector<Double> rstf;
374 t.keywordSet().get("RESTFREQS",rstf);
375 if (rstf.nelements() == 0) {
376 throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
377 }
378 } else if (Unit(un) != u2 && un != "") {
379 throw(AipsError("Unit not conformant with Spectral Coordinates"));
380 }
381 t.rwKeywordSet().define("UNIT", un);
382//
383 MFrequency::Types mdr;
384 if (!MFrequency::getType(mdr, rfrm)) {
385
386 Int a,b;const uInt* c;
387 const String* valid = MFrequency::allMyTypes(a, b, c);
388 String pfix = "Please specify a legal frame type. Types are\n";
389 throw(AipsError(pfix+(*valid)));
390 } else {
391 t.rwKeywordSet().define("REFFRAME",rfrm);
392 }
393//
394 MDoppler::Types dtype;
395 dpl.upcase();
396 if (!MDoppler::getType(dtype, dpl)) {
397 throw(AipsError("Doppler type unknown"));
398 } else {
399 t.rwKeywordSet().define("DOPPLER",dpl);
400 }
401//
402 if (!MFrequency::getType(mdr, brfrm)) {
403 Int a,b;const uInt* c;
404 const String* valid = MFrequency::allMyTypes(a, b, c);
405 String pfix = "Please specify a legal frame type. Types are\n";
406 throw(AipsError(pfix+(*valid)));
407 } else {
408 t.rwKeywordSet().define("BASEREFFRAME",brfrm);
409 }
410}
411
412
413std::vector<double> SDMemTable::getAbcissa(Int whichRow) const
414{
415 std::vector<double> abc(nChan());
416
417// Get header units keyword
418
419 Table t = table_.keywordSet().asTable("FREQUENCIES");
420 String sunit;
421 t.keywordSet().get("UNIT",sunit);
422 if (sunit == "") sunit = "pixel";
423 Unit u(sunit);
424
425// Easy if just wanting pixels
426
427 if (sunit==String("pixel")) {
428 // assume channels/pixels
429 std::vector<double>::iterator it;
430 uInt i=0;
431 for (it = abc.begin(); it != abc.end(); ++it) {
432 (*it) = Double(i++);
433 }
434//
435 return abc;
436 }
437
438// Continue with km/s or Hz. Get FreqIDs
439
440 Vector<uInt> freqIDs;
441 freqidCol_.get(whichRow, freqIDs);
442 uInt freqID = freqIDs(IFSel_);
443 restfreqidCol_.get(whichRow, freqIDs);
444 uInt restFreqID = freqIDs(IFSel_);
445
446// Get SpectralCoordinate, set reference frame conversion,
447// velocity conversion, and rest freq state
448
449 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
450 Vector<Double> pixel(nChan());
451 indgen(pixel);
452//
453 if (u == Unit("km/s")) {
454 Vector<Double> world;
455 spc.pixelToVelocity(world,pixel);
456 std::vector<double>::iterator it;
457 uInt i = 0;
458 for (it = abc.begin(); it != abc.end(); ++it) {
459 (*it) = world[i];
460 i++;
461 }
462 } else if (u == Unit("Hz")) {
463
464// Set world axis units
465
466 Vector<String> wau(1); wau = u.getName();
467 spc.setWorldAxisUnits(wau);
468//
469 std::vector<double>::iterator it;
470 Double tmp;
471 uInt i = 0;
472 for (it = abc.begin(); it != abc.end(); ++it) {
473 spc.toWorld(tmp,pixel[i]);
474 (*it) = tmp;
475 i++;
476 }
477 }
478 return abc;
479}
480
481std::string SDMemTable::getAbcissaString(Int whichRow) const
482{
483 Table t = table_.keywordSet().asTable("FREQUENCIES");
484//
485 String sunit;
486 t.keywordSet().get("UNIT",sunit);
487 if (sunit == "") sunit = "pixel";
488 Unit u(sunit);
489//
490 Vector<uInt> freqIDs;
491 freqidCol_.get(whichRow, freqIDs);
492 uInt freqID = freqIDs(IFSel_);
493 restfreqidCol_.get(whichRow, freqIDs);
494 uInt restFreqID = freqIDs(IFSel_);
495
496// Get SpectralCoordinate, with frame, velocity, rest freq state set
497
498 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
499//
500 String s = "Channel";
501 if (u == Unit("km/s")) {
502 s = CoordinateUtil::axisLabel(spc,0,True,True,True);
503 } else if (u == Unit("Hz")) {
504 Vector<String> wau(1);wau = u.getName();
505 spc.setWorldAxisUnits(wau);
506//
507 s = CoordinateUtil::axisLabel(spc,0,True,True,False);
508 }
509 return s;
510}
511
512void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow)
513{
514 Array<Float> arr;
515 specCol_.get(whichRow, arr);
516 if (spectrum.size() != arr.shape()(3)) {
517 throw(AipsError("Attempting to set spectrum with incorrect length."));
518 }
519
520 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
521 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
522 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
523 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
524 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
525 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
526
527 std::vector<float>::iterator it = spectrum.begin();
528 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
529 (*i) = Float(*it);
530 it++;
531 }
532 specCol_.put(whichRow, arr);
533}
534
535void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const
536{
537 Array<Float> arr;
538 specCol_.get(whichRow, arr);
539 spectrum.resize(arr.shape()(3));
540 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
541 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
542 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
543 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
544 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
545 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
546
547 ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
548 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
549 (*va) = (*i);
550 va++;
551 }
552}
553/*
554void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
555 Array<uChar> arr;
556 flagsCol_.get(whichRow, arr);
557 mask.resize(arr.shape()(3));
558
559 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
560 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
561 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
562 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
563 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
564 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
565
566 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
567
568 ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask);
569 std::vector<bool> tmp;
570 tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
571 // iterator should work on chanMask_??
572 std::vector<bool>::iterator miter;
573 miter = tmp.begin();
574
575 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
576 bool out =!static_cast<bool>(*i);
577 if (useUserMask) {
578 out = out && (*miter);
579 miter++;
580 }
581 (*va) = out;
582 va++;
583 }
584}
585*/
586MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
587 Bool useSelection) const
588{
589 Array<Float> arr;
590 Array<uChar> farr;
591 specCol_.get(whichRow, arr);
592 flagsCol_.get(whichRow, farr);
593 Array<Bool> barr(farr.shape());convertArray(barr, farr);
594 MaskedArray<Float> marr;
595 if (useSelection) {
596 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
597 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
598 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
599 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
600 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
601 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
602
603 ArrayAccessor<Bool, Axis<asap::BeamAxis> > baa0(barr);
604 baa0.reset(baa0.begin(uInt(beamSel_)));//go to beam
605 ArrayAccessor<Bool, Axis<asap::IFAxis> > baa1(baa0);
606 baa1.reset(baa1.begin(uInt(IFSel_)));// go to IF
607 ArrayAccessor<Bool, Axis<asap::PolAxis> > baa2(baa1);
608 baa2.reset(baa2.begin(uInt(polSel_)));// go to pol
609
610 Vector<Float> a(arr.shape()(3));
611 Vector<Bool> b(barr.shape()(3));
612 ArrayAccessor<Float, Axis<asap::BeamAxis> > a0(a);
613 ArrayAccessor<Bool, Axis<asap::BeamAxis> > b0(b);
614
615 ArrayAccessor<Bool, Axis<asap::ChanAxis> > j(baa2);
616 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2);
617 i != i.end(); ++i) {
618 (*a0) = (*i);
619 (*b0) = !(*j);
620 j++;
621 a0++;
622 b0++;
623 }
624 marr.setData(a,b);
625 } else {
626 marr.setData(arr,!barr);
627 }
628 return marr;
629}
630
631Float SDMemTable::getTsys(Int whichRow) const
632{
633 Array<Float> arr;
634 tsCol_.get(whichRow, arr);
635 Float out;
636 IPosition ip(arr.shape());
637 ip(0) = beamSel_;ip(1) = IFSel_;ip(2) = polSel_;ip(3)=0;
638 out = arr(ip);
639 return out;
640}
641
642MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const
643{
644 MDirection::Types mdr = getDirectionReference();
645 Array<Double> posit;
646 dirCol_.get(whichRow,posit);
647 Vector<Double> wpos(2);
648 Int rb;
649 rbeamCol_.get(whichRow,rb);
650 wpos[0] = posit(IPosition(2,beamSel_,0));
651 wpos[1] = posit(IPosition(2,beamSel_,1));
652 if (refBeam && rb > -1) { // use refBeam instead if it exists
653 wpos[0] = posit(IPosition(2,rb,0));
654 wpos[1] = posit(IPosition(2,rb,1));
655 }
656
657 Quantum<Double> lon(wpos[0],Unit(String("rad")));
658 Quantum<Double> lat(wpos[1],Unit(String("rad")));
659 return MDirection(lon, lat, mdr);
660}
661
662MEpoch SDMemTable::getEpoch(Int whichRow) const
663{
664 MEpoch::Types met = getTimeReference();
665//
666 Double obstime;
667 timeCol_.get(whichRow,obstime);
668 MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
669 return MEpoch(tm2, met);
670}
671
672MPosition SDMemTable::getAntennaPosition () const
673{
674 Vector<Double> antpos;
675 table_.keywordSet().get("AntennaPosition", antpos);
676 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
677 return MPosition(mvpos);
678}
679
680
681SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const
682{
683
684 Table t = table_.keywordSet().asTable("FREQUENCIES");
685 if (freqID> t.nrow() ) {
686 throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range"));
687 }
688
689 Double rp,rv,inc;
690 String rf;
691 ROScalarColumn<Double> rpc(t, "REFPIX");
692 ROScalarColumn<Double> rvc(t, "REFVAL");
693 ROScalarColumn<Double> incc(t, "INCREMENT");
694 t.keywordSet().get("BASEREFFRAME",rf);
695
696// Create SpectralCoordinate (units Hz)
697
698 MFrequency::Types mft;
699 if (!MFrequency::getType(mft, rf)) {
700 cerr << "Frequency type unknown assuming TOPO" << endl;
701 mft = MFrequency::TOPO;
702 }
703 rpc.get(freqID, rp);
704 rvc.get(freqID, rv);
705 incc.get(freqID, inc);
706//
707 SpectralCoordinate spec(mft,rv,inc,rp);
708//
709 return spec;
710}
711
712
713SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID, uInt restFreqID,
714 uInt whichRow) const
715{
716
717// Create basic SC
718
719 SpectralCoordinate spec = getSpectralCoordinate (freqID);
720//
721 Table t = table_.keywordSet().asTable("FREQUENCIES");
722
723// Set rest frequency
724
725 Vector<Double> restFreqIDs;
726 t.keywordSet().get("RESTFREQS",restFreqIDs);
727 if (restFreqID < restFreqIDs.nelements()) { // SHould always be true
728 spec.setRestFrequency(restFreqIDs(restFreqID),True);
729 }
730
731// Set up frame conversion layer
732
733 String frm;
734 t.keywordSet().get("REFFRAME",frm);
735 if (frm == "") frm = "TOPO";
736 MFrequency::Types mtype;
737 if (!MFrequency::getType(mtype, frm)) {
738 cerr << "Frequency type unknown assuming TOPO" << endl; // SHould never happen
739 mtype = MFrequency::TOPO;
740 }
741
742// Set reference frame conversion (requires row)
743
744 MDirection dir = getDirection(whichRow);
745 MEpoch epoch = getEpoch(whichRow);
746 MPosition pos = getAntennaPosition();
747//
748 if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) {
749 throw(AipsError("Couldn't convert frequency frame."));
750 }
751
752// Now velocity conversion if appropriate
753
754 String unitStr;
755 t.keywordSet().get("UNIT",unitStr);
756//
757 String dpl;
758 t.keywordSet().get("DOPPLER",dpl);
759 MDoppler::Types dtype;
760 MDoppler::getType(dtype, dpl);
761
762// Only set velocity unit if non-blank and non-Hz
763
764 if (!unitStr.empty()) {
765 Unit unitU(unitStr);
766 if (unitU==Unit("Hz")) {
767 } else {
768 spec.setVelocity(unitStr, dtype);
769 }
770 }
771//
772 return spec;
773}
774
775
776Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
777 uInt freqID) {
778 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
779 if (freqID > t.nrow() ) {
780 throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
781 }
782 ScalarColumn<Double> rpc(t, "REFPIX");
783 ScalarColumn<Double> rvc(t, "REFVAL");
784 ScalarColumn<Double> incc(t, "INCREMENT");
785
786 rpc.put(freqID, speccord.referencePixel()[0]);
787 rvc.put(freqID, speccord.referenceValue()[0]);
788 incc.put(freqID, speccord.increment()[0]);
789
790 return True;
791}
792
793Int SDMemTable::nCoordinates() const
794{
795 return table_.keywordSet().asTable("FREQUENCIES").nrow();
796}
797
798
799std::vector<double> SDMemTable::getRestFreqs() const
800{
801 Table t = table_.keywordSet().asTable("FREQUENCIES");
802 Vector<Double> tvec;
803 t.keywordSet().get("RESTFREQS",tvec);
804 std::vector<double> stlout;
805 tvec.tovector(stlout);
806 return stlout;
807}
808
809bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft)
810{
811 TableDesc td("", "1", TableDesc::Scratch);
812 td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
813 td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
814 td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
815 SetupNewTable aNewTab("freqs", td, Table::New);
816 Table aTable (aNewTab, Table::Memory, sdft.length());
817 ScalarColumn<Double> sc0(aTable, "REFPIX");
818 ScalarColumn<Double> sc1(aTable, "REFVAL");
819 ScalarColumn<Double> sc2(aTable, "INCREMENT");
820 for (uInt i=0; i < sdft.length(); ++i) {
821 sc0.put(i,sdft.referencePixel(i));
822 sc1.put(i,sdft.referenceValue(i));
823 sc2.put(i,sdft.increment(i));
824 }
825 String rf = sdft.refFrame();
826 if (rf.contains("TOPO")) rf = "TOPO";
827
828 aTable.rwKeywordSet().define("BASEREFFRAME", rf);
829 aTable.rwKeywordSet().define("REFFRAME", rf);
830 aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
831 aTable.rwKeywordSet().define("UNIT", String(""));
832 aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
833 Vector<Double> rfvec;
834 String rfunit;
835 sdft.restFrequencies(rfvec,rfunit);
836 Quantum<Vector<Double> > q(rfvec, rfunit);
837 rfvec.resize();
838 rfvec = q.getValue("Hz");
839 aTable.rwKeywordSet().define("RESTFREQS", rfvec);
840 table_.rwKeywordSet().defineTable ("FREQUENCIES", aTable);
841 return True;
842}
843
844SDFrequencyTable SDMemTable::getSDFreqTable() const
845{
846 const Table& t = table_.keywordSet().asTable("FREQUENCIES");
847 SDFrequencyTable sdft;
848
849// Add refpix/refval/incr. What are the units ? Hz I suppose
850// but it's nowhere described...
851
852 Vector<Double> refPix, refVal, incr;
853 ScalarColumn<Double> refPixCol(t, "REFPIX");
854 ScalarColumn<Double> refValCol(t, "REFVAL");
855 ScalarColumn<Double> incrCol(t, "INCREMENT");
856 refPix = refPixCol.getColumn();
857 refVal = refValCol.getColumn();
858 incr = incrCol.getColumn();
859//
860 uInt n = refPix.nelements();
861 for (uInt i=0; i<n; i++) {
862 sdft.addFrequency(refPix[i], refVal[i], incr[i]);
863 }
864
865// Frequency reference frame. I don't know if this
866// is the correct frame. It might be 'REFFRAME'
867// rather than 'BASEREFFRAME' ?
868
869 String baseFrame;
870 t.keywordSet().get("BASEREFFRAME",baseFrame);
871 sdft.setRefFrame(baseFrame);
872
873// Equinox
874
875 Float equinox;
876 t.keywordSet().get("EQUINOX", equinox);
877 sdft.setEquinox(equinox);
878
879// Rest Frequency
880
881 Vector<Double> restFreqs;
882 t.keywordSet().get("RESTFREQS", restFreqs);
883 for (uInt i=0; i<restFreqs.nelements(); i++) {
884 sdft.addRestFrequency(restFreqs[i]);
885 }
886 sdft.setRestFrequencyUnit(String("Hz"));
887//
888 return sdft;
889}
890
891bool SDMemTable::putSDContainer(const SDContainer& sdc)
892{
893 uInt rno = table_.nrow();
894 table_.addRow();
895//
896 timeCol_.put(rno, sdc.timestamp);
897 srcnCol_.put(rno, sdc.sourcename);
898 fldnCol_.put(rno, sdc.fieldname);
899 specCol_.put(rno, sdc.getSpectrum());
900 flagsCol_.put(rno, sdc.getFlags());
901 tsCol_.put(rno, sdc.getTsys());
902 scanCol_.put(rno, sdc.scanid);
903 integrCol_.put(rno, sdc.interval);
904 freqidCol_.put(rno, sdc.getFreqMap());
905 restfreqidCol_.put(rno, sdc.getRestFreqMap());
906 dirCol_.put(rno, sdc.getDirection());
907 rbeamCol_.put(rno, sdc.refbeam);
908 tcalCol_.put(rno, sdc.tcal);
909 tcaltCol_.put(rno, sdc.tcaltime);
910 azCol_.put(rno, sdc.azimuth);
911 elCol_.put(rno, sdc.elevation);
912 paraCol_.put(rno, sdc.parangle);
913 histCol_.put(rno, sdc.getHistory());
914//
915 return true;
916}
917
918SDContainer SDMemTable::getSDContainer(uInt whichRow) const
919{
920 SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
921 timeCol_.get(whichRow, sdc.timestamp);
922 srcnCol_.get(whichRow, sdc.sourcename);
923 integrCol_.get(whichRow, sdc.interval);
924 scanCol_.get(whichRow, sdc.scanid);
925 fldnCol_.get(whichRow, sdc.fieldname);
926 rbeamCol_.get(whichRow, sdc.refbeam);
927 azCol_.get(whichRow, sdc.azimuth);
928 elCol_.get(whichRow, sdc.elevation);
929 paraCol_.get(whichRow, sdc.parangle);
930 Vector<Float> tc;
931 tcalCol_.get(whichRow, tc);
932 sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
933 tcaltCol_.get(whichRow, sdc.tcaltime);
934//
935 Array<Float> spectrum;
936 Array<Float> tsys;
937 Array<uChar> flagtrum;
938 Vector<uInt> fmap;
939 Array<Double> direction;
940 Vector<String> histo;
941//
942 specCol_.get(whichRow, spectrum);
943 sdc.putSpectrum(spectrum);
944 flagsCol_.get(whichRow, flagtrum);
945 sdc.putFlags(flagtrum);
946 tsCol_.get(whichRow, tsys);
947 sdc.putTsys(tsys);
948 freqidCol_.get(whichRow, fmap);
949 sdc.putFreqMap(fmap);
950 restfreqidCol_.get(whichRow, fmap);
951 sdc.putRestFreqMap(fmap);
952 dirCol_.get(whichRow, direction);
953 sdc.putDirection(direction);
954 histCol_.get(whichRow, histo);
955 sdc.putHistory(histo);
956 return sdc;
957}
958
959bool SDMemTable::putSDHeader(const SDHeader& sdh)
960{
961 table_.rwKeywordSet().define("nIF", sdh.nif);
962 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
963 table_.rwKeywordSet().define("nPol", sdh.npol);
964 table_.rwKeywordSet().define("nChan", sdh.nchan);
965 table_.rwKeywordSet().define("Observer", sdh.observer);
966 table_.rwKeywordSet().define("Project", sdh.project);
967 table_.rwKeywordSet().define("Obstype", sdh.obstype);
968 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
969 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
970 table_.rwKeywordSet().define("Equinox", sdh.equinox);
971 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
972 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
973 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
974 table_.rwKeywordSet().define("UTC", sdh.utc);
975 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
976 table_.rwKeywordSet().define("Epoch", sdh.epoch);
977 return true;
978}
979
980SDHeader SDMemTable::getSDHeader() const
981{
982 SDHeader sdh;
983 table_.keywordSet().get("nBeam",sdh.nbeam);
984 table_.keywordSet().get("nIF",sdh.nif);
985 table_.keywordSet().get("nPol",sdh.npol);
986 table_.keywordSet().get("nChan",sdh.nchan);
987 table_.keywordSet().get("Observer", sdh.observer);
988 table_.keywordSet().get("Project", sdh.project);
989 table_.keywordSet().get("Obstype", sdh.obstype);
990 table_.keywordSet().get("AntennaName", sdh.antennaname);
991 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
992 table_.keywordSet().get("Equinox", sdh.equinox);
993 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
994 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
995 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
996 table_.keywordSet().get("UTC", sdh.utc);
997 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
998 table_.keywordSet().get("Epoch", sdh.epoch);
999 return sdh;
1000}
1001void SDMemTable::makePersistent(const std::string& filename)
1002{
1003 table_.deepCopy(filename,Table::New);
1004}
1005
1006Int SDMemTable::nScan() const {
1007 Int n = 0;
1008 Int previous = -1;Int current=0;
1009 for (uInt i=0; i< scanCol_.nrow();i++) {
1010 scanCol_.getScalar(i,current);
1011 if (previous != current) {
1012 previous = current;
1013 n++;
1014 }
1015 }
1016 return n;
1017}
1018
1019String SDMemTable::formatSec(Double x) const
1020{
1021 Double xcop = x;
1022 MVTime mvt(xcop/24./3600.); // make days
1023
1024 if (x < 59.95)
1025 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
1026 else if (x < 3599.95)
1027 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
1028 else {
1029 ostringstream oss;
1030 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
1031 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
1032 return String(oss);
1033 }
1034};
1035
1036String SDMemTable::formatDirection(const MDirection& md) const
1037{
1038 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
1039 Int prec = 7;
1040
1041 MVAngle mvLon(t[0]);
1042 String sLon = mvLon.string(MVAngle::TIME,prec);
1043 MVAngle mvLat(t[1]);
1044 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
1045 return sLon + String(" ") + sLat;
1046}
1047
1048
1049std::string SDMemTable::getFluxUnit() const
1050{
1051 String tmp;
1052 table_.keywordSet().get("FluxUnit", tmp);
1053 return tmp;
1054}
1055
1056void SDMemTable::setFluxUnit(const std::string& unit)
1057{
1058 String tmp(unit);
1059 Unit tU(tmp);
1060 if (tU==Unit("K") || tU==Unit("Jy")) {
1061 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
1062 } else {
1063 throw AipsError("Illegal unit - must be compatible with Jy or K");
1064 }
1065}
1066
1067
1068void SDMemTable::setInstrument(const std::string& name)
1069{
1070 Bool throwIt = True;
1071 Instrument ins = convertInstrument (name, throwIt);
1072 String nameU(name);
1073 nameU.upcase();
1074 table_.rwKeywordSet().define(String("AntennaName"), nameU);
1075}
1076
1077std::string SDMemTable::summary(bool verbose) const {
1078
1079// Format header info
1080
1081 ostringstream oss;
1082 oss << endl;
1083 oss << "--------------------------------------------------------------------------------" << endl;
1084 oss << " Scan Table Summary" << endl;
1085 oss << "--------------------------------------------------------------------------------" << endl;
1086 oss.flags(std::ios_base::left);
1087 oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
1088 << setw(15) << "IFs:" << setw(4) << nIF() << endl
1089 << setw(15) << "Polarisations:" << setw(4) << nPol() << endl
1090 << setw(15) << "Channels:" << setw(4) << nChan() << endl;
1091 oss << endl;
1092 String tmp;
1093 table_.keywordSet().get("Observer", tmp);
1094 oss << setw(15) << "Observer:" << tmp << endl;
1095 oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl;
1096 table_.keywordSet().get("Project", tmp);
1097 oss << setw(15) << "Project:" << tmp << endl;
1098 table_.keywordSet().get("Obstype", tmp);
1099 oss << setw(15) << "Obs. Type:" << tmp << endl;
1100 table_.keywordSet().get("AntennaName", tmp);
1101 oss << setw(15) << "Antenna Name:" << tmp << endl;
1102 table_.keywordSet().get("FluxUnit", tmp);
1103 oss << setw(15) << "Flux Unit:" << tmp << endl;
1104 Table t = table_.keywordSet().asTable("FREQUENCIES");
1105 Vector<Double> vec;
1106 t.keywordSet().get("RESTFREQS",vec);
1107 oss << setw(15) << "Rest Freqs:";
1108 if (vec.nelements() > 0) {
1109 oss << setprecision(0) << vec << " [Hz]" << endl;
1110 } else {
1111 oss << "None set" << endl;
1112 }
1113 oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
1114 oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
1115 << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
1116 oss << endl;
1117//
1118 String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")";
1119 oss << setw(5) << "Scan"
1120 << setw(15) << "Source"
1121 << setw(24) << dirtype
1122 << setw(10) << "Time"
1123 << setw(18) << "Integration"
1124 << setw(7) << "FreqIDs" << endl;
1125 oss << "--------------------------------------------------------------------------------" << endl;
1126
1127// Generate list of scan start and end integrations
1128
1129 Vector<Int> scanIDs = scanCol_.getColumn();
1130 Vector<uInt> startInt, endInt;
1131 mathutil::scanBoundaries(startInt, endInt, scanIDs);
1132//
1133 const uInt nScans = startInt.nelements();
1134 String name;
1135 Vector<uInt> freqIDs, listFQ;
1136 uInt nInt;
1137//
1138 for (uInt i=0; i<nScans; i++) {
1139
1140// Get things from first integration of scan
1141
1142 String time = getTime(startInt(i),False);
1143 String tInt = formatSec(Double(getInterval(startInt(i))));
1144 String posit = formatDirection(getDirection(startInt(i),True));
1145 srcnCol_.getScalar(startInt(i),name);
1146
1147// Find all the FreqIDs in this scan
1148
1149 listFQ.resize(0);
1150 for (uInt j=startInt(i); j<endInt(i)+1; j++) {
1151 freqidCol_.get(j, freqIDs);
1152 for (uInt k=0; k<freqIDs.nelements(); k++) {
1153 mathutil::addEntry(listFQ, freqIDs(k));
1154 }
1155 }
1156//
1157 nInt = endInt(i) - startInt(i) + 1;
1158 oss << setw(3) << std::right << i << std::left << setw(2) << " "
1159 << setw(15) << name
1160 << setw(24) << posit
1161 << setw(10) << time
1162 << setw(3) << std::right << nInt << setw(3) << " x " << std::left
1163 << setw(6) << tInt
1164 << " " << listFQ << endl;
1165 }
1166 oss << endl;
1167 oss << "Table contains " << table_.nrow() << " integration(s) in " << nScans << " scan(s)." << endl;
1168
1169// Frequency Table
1170
1171 if (verbose) {
1172 std::vector<string> info = getCoordInfo();
1173 SDFrequencyTable sdft = getSDFreqTable();
1174 oss << endl << endl;
1175 oss << "FreqID Frame RefFreq(Hz) RefPix Increment(Hz)" << endl;
1176 oss << "--------------------------------------------------------------------------------" << endl;
1177 for (uInt i=0; i<sdft.length(); i++) {
1178 oss << setw(8) << i << setw(8)
1179 << info[3] << setw(16) << setprecision(8)
1180 << sdft.referenceValue(i) << setw(10)
1181 << sdft.referencePixel(i) << setw(12)
1182 << sdft.increment(i) << endl;
1183 }
1184 oss << "--------------------------------------------------------------------------------" << endl;
1185 }
1186 return String(oss);
1187}
1188
1189Int SDMemTable::nBeam() const
1190{
1191 Int n;
1192 table_.keywordSet().get("nBeam",n);
1193 return n;
1194}
1195Int SDMemTable::nIF() const {
1196 Int n;
1197 table_.keywordSet().get("nIF",n);
1198 return n;
1199}
1200Int SDMemTable::nPol() const {
1201 Int n;
1202 table_.keywordSet().get("nPol",n);
1203 return n;
1204}
1205Int SDMemTable::nChan() const {
1206 Int n;
1207 table_.keywordSet().get("nChan",n);
1208 return n;
1209}
1210bool SDMemTable::appendHistory(const std::string& hist, int whichRow)
1211{
1212 Vector<String> history;
1213 histCol_.get(whichRow, history);
1214 history.resize(history.nelements()+1,True);
1215 history[history.nelements()-1] = hist;
1216 histCol_.put(whichRow, history);
1217}
1218
1219std::vector<std::string> SDMemTable::history(int whichRow) const
1220{
1221 Vector<String> history;
1222 histCol_.get(whichRow, history);
1223 std::vector<std::string> stlout;
1224 // there is no Array<String>.tovector(std::vector<std::string>), so
1225 // do it by hand
1226 for (uInt i=0; i<history.nelements(); ++i) {
1227 stlout.push_back(history[i]);
1228 }
1229 return stlout;
1230}
1231/*
1232void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
1233
1234 std::vector<int>::iterator it;
1235 ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_);
1236 for (it = whichChans.begin(); it != whichChans.end(); it++) {
1237 j.reset(j.begin(uInt(*it)));
1238 for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) {
1239 for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) {
1240 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii);
1241 iii != iii.end(); ++iii) {
1242 (*iii) =
1243 }
1244 }
1245 }
1246 }
1247
1248}
1249*/
1250void SDMemTable::flag(int whichRow)
1251{
1252 Array<uChar> arr;
1253 flagsCol_.get(whichRow, arr);
1254
1255 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
1256 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
1257 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
1258 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
1259 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
1260 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
1261
1262 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1263 (*i) = uChar(True);
1264 }
1265
1266 flagsCol_.put(whichRow, arr);
1267}
1268
1269MDirection::Types SDMemTable::getDirectionReference() const
1270{
1271 Float eq;
1272 table_.keywordSet().get("Equinox",eq);
1273 std::map<float,string> mp;
1274 mp[2000.0] = "J2000";
1275 mp[1950.0] = "B1950";
1276 MDirection::Types mdr;
1277 if (!MDirection::getType(mdr, mp[eq])) {
1278 mdr = MDirection::J2000;
1279 cerr << "Unknown equinox using J2000" << endl;
1280 }
1281//
1282 return mdr;
1283}
1284
1285MEpoch::Types SDMemTable::getTimeReference() const
1286{
1287 MEpoch::Types met;
1288 String ep;
1289 table_.keywordSet().get("Epoch",ep);
1290 if (!MEpoch::getType(met, ep)) {
1291 cerr << "Epoch type unknown - using UTC" << endl;
1292 met = MEpoch::UTC;
1293 }
1294//
1295 return met;
1296}
1297
1298
1299Instrument SDMemTable::convertInstrument(const String& instrument,
1300 Bool throwIt)
1301{
1302 String t(instrument);
1303 t.upcase();
1304
1305// The strings are what SDReader returns, after cunning interrogation
1306// of station names... :-(
1307
1308 Instrument inst = asap::UNKNOWN;
1309 if (t==String("DSS-43")) {
1310 inst = TIDBINBILLA;
1311 } else if (t==String("ATPKSMB")) {
1312 inst = ATPKSMB;
1313 } else if (t==String("ATPKSHOH")) {
1314 inst = ATPKSHOH;
1315 } else if (t==String("ATMOPRA")) {
1316 inst = ATMOPRA;
1317 } else if (t==String("CEDUNA")) {
1318 inst = CEDUNA;
1319 } else if (t==String("HOBART")) {
1320 inst = HOBART;
1321 } else {
1322 if (throwIt) {
1323 throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
1324 }
1325 }
1326 return inst;
1327}
1328
1329Bool SDMemTable::setRestFreqs (const Vector<Double>& restFreqsIn,
1330 const String& sUnit,
1331 const vector<string>& lines,
1332 const String& source,
1333 Int whichIF)
1334{
1335 const Int nIFs = nIF();
1336 if (whichIF>=nIFs) {
1337 throw(AipsError("Illegal IF index"));
1338 }
1339
1340// FInd vector of restfrequencies
1341// Double takes precedence over String
1342
1343 Unit unit;
1344 Vector<Double> restFreqs;
1345 if (restFreqsIn.nelements()>0) {
1346 restFreqs.resize(restFreqsIn.nelements());
1347 restFreqs = restFreqsIn;
1348 unit = Unit(sUnit);
1349 } else if (lines.size()>0) {
1350 const uInt nLines = lines.size();
1351 unit = Unit("Hz");
1352 restFreqs.resize(nLines);
1353 MFrequency lineFreq;
1354 for (uInt i=0; i<nLines; i++) {
1355 String tS(lines[i]);
1356 tS.upcase();
1357 if (MeasTable::Line(lineFreq, tS)) {
1358 restFreqs[i] = lineFreq.getValue().getValue(); // Hz
1359 } else {
1360 String s = String(lines[i]) + String(" is an unrecognized spectral line");
1361 throw(AipsError(s));
1362 }
1363 }
1364 } else {
1365 throw(AipsError("You have not specified any rest frequencies or lines"));
1366 }
1367
1368// If multiple restfreqs, must be length nIF. In this
1369// case we will just replace the rest frequencies
1370//
1371
1372 const uInt nRestFreqs = restFreqs.nelements();
1373 Int idx = -1;
1374 SDFrequencyTable sdft = getSDFreqTable();
1375
1376 if (nRestFreqs>1) {
1377
1378// Replace restFreqs, one per IF
1379
1380 if (nRestFreqs != nIFs) {
1381 throw (AipsError("Number of rest frequencies must be equal to the number of IFs"));
1382 }
1383 cerr << "Replacing rest frequencies with given list, one per IF" << endl;
1384//
1385 sdft.deleteRestFrequencies();
1386 for (uInt i=0; i<nRestFreqs; i++) {
1387 Quantum<Double> rf(restFreqs[i], unit);
1388 sdft.addRestFrequency(rf.getValue("Hz"));
1389 }
1390 } else {
1391
1392// Add new rest freq
1393
1394 Quantum<Double> rf(restFreqs[0], unit);
1395 idx = sdft.addRestFrequency(rf.getValue("Hz"));
1396 cerr << "Selecting given rest frequency" << endl;
1397 }
1398
1399// Replace
1400
1401 Bool empty = source.empty();
1402 Bool ok = False;
1403 if (putSDFreqTable(sdft)) {
1404 const uInt nRow = table_.nrow();
1405 String srcName;
1406 Vector<uInt> restFreqIDs;
1407 for (uInt i=0; i<nRow; i++) {
1408 srcnCol_.get(i, srcName);
1409 restfreqidCol_.get(i,restFreqIDs);
1410//
1411 if (idx==-1) {
1412
1413// Replace vector of restFreqs; one per IF.
1414// No selection possible
1415
1416 for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i;
1417 } else {
1418
1419// Set RestFreqID for selected data
1420
1421 if (empty || source==srcName) {
1422 if (whichIF<0) {
1423 restFreqIDs = idx;
1424 } else {
1425 restFreqIDs[whichIF] = idx;
1426 }
1427 }
1428 }
1429//
1430 restfreqidCol_.put(i,restFreqIDs);
1431 }
1432 ok = True;
1433 } else {
1434 ok = False;
1435 }
1436//
1437 return ok;
1438}
1439
1440void SDMemTable::spectralLines () const
1441{
1442 Vector<String> lines = MeasTable::Lines();
1443 MFrequency lineFreq;
1444 Double freq;
1445//
1446 cerr.flags(std::ios_base::left);
1447 cerr << "Line Frequency (Hz)" << endl;
1448 cerr << "-----------------------" << endl;
1449 for (uInt i=0; i<lines.nelements(); i++) {
1450 MeasTable::Line(lineFreq, lines[i]);
1451 freq = lineFreq.getValue().getValue(); // Hz
1452//
1453 cerr << setw(11) << lines[i] << setprecision(10) << freq << endl;
1454 }
1455}
1456
1457void SDMemTable::renumber()
1458{
1459 uInt nRow = scanCol_.nrow();
1460 Int newscanid = 0;
1461 Int cIdx;// the current scanid
1462 // get the first scanid
1463 scanCol_.getScalar(0,cIdx);
1464 Int pIdx = cIdx;// the scanid of the previous row
1465 for (uInt i=0; i<nRow;++i) {
1466 scanCol_.getScalar(i,cIdx);
1467 if (pIdx == cIdx) {
1468 // renumber
1469 scanCol_.put(i,newscanid);
1470 } else {
1471 ++newscanid;
1472 pIdx = cIdx; // store scanid
1473 --i; // don't increment next loop
1474 }
1475 }
1476}
1477
Note: See TracBrowser for help on using the repository browser.