source: trunk/src/SDMemTable.cc@ 471

Last change on this file since 471 was 465, checked in by mar637, 20 years ago

Added SDFitTable to handle fits and expose them to python vi the sdfit class.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 49.9 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/VectorSTLIterator.h>
43#include <casa/Arrays/Vector.h>
44#include <casa/BasicMath/Math.h>
45#include <casa/Quanta/MVAngle.h>
46
47#include <tables/Tables/TableParse.h>
48#include <tables/Tables/TableDesc.h>
49#include <tables/Tables/SetupNewTab.h>
50#include <tables/Tables/ScaColDesc.h>
51#include <tables/Tables/ArrColDesc.h>
52
53#include <tables/Tables/ExprNode.h>
54#include <tables/Tables/TableRecord.h>
55#include <measures/Measures/MFrequency.h>
56#include <measures/Measures/MeasTable.h>
57#include <coordinates/Coordinates/CoordinateUtil.h>
58#include <casa/Quanta/MVTime.h>
59#include <casa/Quanta/MVAngle.h>
60
61#include "SDDefs.h"
62#include "SDContainer.h"
63#include "MathUtils.h"
64#include "SDPol.h"
65
66#include "SDMemTable.h"
67
68using namespace casa;
69using namespace asap;
70
71SDMemTable::SDMemTable() :
72 IFSel_(0),
73 beamSel_(0),
74 polSel_(0)
75{
76 setup();
77 attach();
78}
79
80SDMemTable::SDMemTable(const std::string& name) :
81 IFSel_(0),
82 beamSel_(0),
83 polSel_(0)
84{
85 Table tab(name);
86 table_ = tab.copyToMemoryTable("dummy");
87 //cerr << "hello from C SDMemTable @ " << this << endl;
88 attach();
89}
90
91SDMemTable::SDMemTable(const SDMemTable& other, Bool clear)
92{
93 IFSel_= other.IFSel_;
94 beamSel_= other.beamSel_;
95 polSel_= other.polSel_;
96 chanMask_ = other.chanMask_;
97 table_ = other.table_.copyToMemoryTable(String("dummy"));
98 // clear all rows()
99 if (clear) {
100 table_.removeRow(this->table_.rowNumbers());
101 } else {
102 IFSel_ = other.IFSel_;
103 beamSel_ = other.beamSel_;
104 polSel_ = other.polSel_;
105 }
106//
107 attach();
108 //cerr << "hello from CC SDMemTable @ " << this << endl;
109}
110
111SDMemTable::SDMemTable(const Table& tab, const std::string& exprs) :
112 IFSel_(0),
113 beamSel_(0),
114 polSel_(0)
115{
116 Table t = tableCommand(exprs,tab);
117 if (t.nrow() == 0)
118 throw(AipsError("Query unsuccessful."));
119 table_ = t.copyToMemoryTable("dummy");
120 attach();
121 renumber();
122}
123
124SDMemTable::~SDMemTable()
125{
126 //cerr << "goodbye from SDMemTable @ " << this << endl;
127}
128
129SDMemTable SDMemTable::getScan(Int scanID) const
130{
131 String cond("SELECT * from $1 WHERE SCANID == ");
132 cond += String::toString(scanID);
133 return SDMemTable(table_, cond);
134}
135
136SDMemTable &SDMemTable::operator=(const SDMemTable& other)
137{
138 if (this != &other) {
139 IFSel_= other.IFSel_;
140 beamSel_= other.beamSel_;
141 polSel_= other.polSel_;
142 chanMask_.resize(0);
143 chanMask_ = other.chanMask_;
144 table_ = other.table_.copyToMemoryTable(String("dummy"));
145 attach();
146 }
147 //cerr << "hello from ASS SDMemTable @ " << this << endl;
148 return *this;
149}
150
151SDMemTable SDMemTable::getSource(const std::string& source) const
152{
153 String cond("SELECT * from $1 WHERE SRCNAME == ");
154 cond += source;
155 return SDMemTable(table_, cond);
156}
157
158void SDMemTable::setup()
159{
160 TableDesc td("", "1", TableDesc::Scratch);
161 td.comment() = "A SDMemTable";
162
163 td.addColumn(ScalarColumnDesc<Double>("TIME"));
164 td.addColumn(ScalarColumnDesc<String>("SRCNAME"));
165 td.addColumn(ArrayColumnDesc<Float>("SPECTRA"));
166 td.addColumn(ArrayColumnDesc<uChar>("FLAGTRA"));
167 td.addColumn(ArrayColumnDesc<Float>("TSYS"));
168 td.addColumn(ArrayColumnDesc<Float>("STOKES"));
169 td.addColumn(ScalarColumnDesc<Int>("SCANID"));
170 td.addColumn(ScalarColumnDesc<Double>("INTERVAL"));
171 td.addColumn(ArrayColumnDesc<uInt>("FREQID"));
172 td.addColumn(ArrayColumnDesc<uInt>("RESTFREQID"));
173 td.addColumn(ArrayColumnDesc<Double>("DIRECTION"));
174 td.addColumn(ScalarColumnDesc<String>("FIELDNAME"));
175 td.addColumn(ScalarColumnDesc<String>("TCALTIME"));
176 td.addColumn(ArrayColumnDesc<Float>("TCAL"));
177 td.addColumn(ScalarColumnDesc<Float>("AZIMUTH"));
178 td.addColumn(ScalarColumnDesc<Float>("ELEVATION"));
179 td.addColumn(ScalarColumnDesc<Float>("PARANGLE"));
180 td.addColumn(ScalarColumnDesc<Int>("REFBEAM"));
181 td.addColumn(ArrayColumnDesc<String>("HISTORY"));
182 td.addColumn(ArrayColumnDesc<Int>("FITID"));
183
184
185 // Now create Table SetUp from the description.
186
187 SetupNewTable aNewTab("dummy", td, Table::New);
188
189 // Bind the Stokes Virtual machine to the STOKES column Because we
190 // don't know how many polarizations will be in the data at this
191 // point, we must bind the Virtual Engine regardless. The STOKES
192 // column won't be accessed if not appropriate (nPol=4)
193
194
195 SDStokesEngine::registerClass();
196 SDStokesEngine stokesEngine(String("STOKES"), String("SPECTRA"));
197 aNewTab.bindColumn ("STOKES", stokesEngine);
198
199 // Create Table
200 table_ = Table(aNewTab, Table::Memory, 0);
201 // add subtable
202 TableDesc tdf("", "1", TableDesc::Scratch);
203 tdf.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
204 tdf.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
205 tdf.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
206 tdf.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
207 tdf.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
208 SetupNewTable fittab("fits", tdf, Table::New);
209 Table fitTable(fittab, Table::Memory);
210 table_.rwKeywordSet().defineTable("FITS", fitTable);
211
212
213}
214
215void SDMemTable::attach()
216{
217 timeCol_.attach(table_, "TIME");
218 srcnCol_.attach(table_, "SRCNAME");
219 specCol_.attach(table_, "SPECTRA");
220 flagsCol_.attach(table_, "FLAGTRA");
221 tsCol_.attach(table_, "TSYS");
222 stokesCol_.attach(table_, "STOKES");
223 scanCol_.attach(table_, "SCANID");
224 integrCol_.attach(table_, "INTERVAL");
225 freqidCol_.attach(table_, "FREQID");
226 restfreqidCol_.attach(table_, "RESTFREQID");
227 dirCol_.attach(table_, "DIRECTION");
228 fldnCol_.attach(table_, "FIELDNAME");
229 tcaltCol_.attach(table_, "TCALTIME");
230 tcalCol_.attach(table_, "TCAL");
231 azCol_.attach(table_, "AZIMUTH");
232 elCol_.attach(table_, "ELEVATION");
233 paraCol_.attach(table_, "PARANGLE");
234 rbeamCol_.attach(table_, "REFBEAM");
235 histCol_.attach(table_, "HISTORY");
236 fitCol_.attach(table_,"FITID");
237}
238
239
240std::string SDMemTable::getSourceName(Int whichRow) const
241{
242 String name;
243 srcnCol_.get(whichRow, name);
244 return name;
245}
246
247std::string SDMemTable::getTime(Int whichRow, Bool showDate) const
248{
249 Double tm;
250 if (whichRow > -1) {
251 timeCol_.get(whichRow, tm);
252 } else {
253 table_.keywordSet().get("UTC",tm);
254 }
255 MVTime mvt(tm);
256 if (showDate)
257 mvt.setFormat(MVTime::YMD);
258 else
259 mvt.setFormat(MVTime::TIME);
260 ostringstream oss;
261 oss << mvt;
262 return String(oss);
263}
264
265double SDMemTable::getInterval(Int whichRow) const
266{
267 Double intval;
268 integrCol_.get(whichRow, intval);
269 return intval;
270}
271
272bool SDMemTable::setIF(Int whichIF)
273{
274 if ( whichIF >= 0 && whichIF < nIF()) {
275 IFSel_ = whichIF;
276 return true;
277 }
278 return false;
279}
280
281bool SDMemTable::setBeam(Int whichBeam)
282{
283 if ( whichBeam >= 0 && whichBeam < nBeam()) {
284 beamSel_ = whichBeam;
285 return true;
286 }
287 return false;
288}
289
290bool SDMemTable::setPol(Int whichPol)
291{
292 if ( whichPol >= 0 && whichPol < nPol()) {
293 polSel_ = whichPol;
294 return true;
295 }
296 return false;
297}
298
299void SDMemTable::resetCursor()
300{
301 polSel_ = 0;
302 IFSel_ = 0;
303 beamSel_ = 0;
304}
305
306bool SDMemTable::setMask(std::vector<int> whichChans)
307{
308 std::vector<int>::iterator it;
309 uInt n = flagsCol_.shape(0)(3);
310 if (whichChans.empty()) {
311 chanMask_ = std::vector<bool>(n,true);
312 return true;
313 }
314 chanMask_.resize(n,true);
315 for (it = whichChans.begin(); it != whichChans.end(); ++it) {
316 if (*it < n) {
317 chanMask_[*it] = false;
318 }
319 }
320 return true;
321}
322
323std::vector<bool> SDMemTable::getMask(Int whichRow) const
324{
325
326 std::vector<bool> mask;
327
328 Array<uChar> arr;
329 flagsCol_.get(whichRow, arr);
330
331 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
332 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
333 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
334 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
335 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
336 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
337
338 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
339
340 std::vector<bool> tmp;
341 tmp = chanMask_; // WHY the fxxx do I have to make a copy here
342 std::vector<bool>::iterator miter;
343 miter = tmp.begin();
344
345 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
346 bool out =!static_cast<bool>(*i);
347 if (useUserMask) {
348 out = out && (*miter);
349 miter++;
350 }
351 mask.push_back(out);
352 }
353 return mask;
354}
355
356
357
358std::vector<float> SDMemTable::getSpectrum(Int whichRow) const
359{
360
361 Array<Float> arr;
362 specCol_.get(whichRow, arr);
363 return getFloatSpectrum(arr);
364}
365
366std::vector<float> SDMemTable::getStokesSpectrum(Int whichRow, Bool doPol,
367 Float paOffset) const
368 // Gets
369 // doPol=False : I,Q,U,V
370 // doPol=True : I,P,PA,V ; P = sqrt(Q**2+U**2), PA = 0.5*atan2(Q,U)
371 //
372{
373 AlwaysAssert(asap::nAxes==4,AipsError);
374 if (nPol()!=1 && nPol()!=2 && nPol()!=4) {
375 throw (AipsError("You must have 1,2 or 4 polarizations to get the Stokes parameters"));
376 }
377 Array<Float> arr;
378 stokesCol_.get(whichRow, arr);
379
380 if (doPol && (polSel_==1 || polSel_==2)) { // Q,U --> P, P.A.
381
382 // Set current cursor location
383 const IPosition& shape = arr.shape();
384 IPosition start, end;
385 setCursorSlice (start, end, shape);
386
387 // Get Q and U slices
388
389 Array<Float> Q = SDPolUtil::getStokesSlice(arr,start,end,"Q");
390 Array<Float> U = SDPolUtil::getStokesSlice(arr,start,end,"U");
391
392 // Compute output
393
394 Array<Float> out;
395 if (polSel_==1) { // P
396 out = SDPolUtil::polarizedIntensity(Q,U);
397 } else if (polSel_==2) { // P.A.
398 out = SDPolUtil::positionAngle(Q,U) + paOffset;
399 }
400
401 // Copy to output
402
403 IPosition vecShape(1,shape(asap::ChanAxis));
404 Vector<Float> outV = out.reform(vecShape);
405 std::vector<float> stlout;
406 outV.tovector(stlout);
407 return stlout;
408
409 } else {
410
411 // Selects at the cursor location
412 return getFloatSpectrum(arr);
413 }
414}
415
416std::vector<float> SDMemTable::getCircularSpectrum(Int whichRow,
417 Bool doRR) const
418 // Gets
419 // RR = I + V
420 // LL = I - V
421{
422 AlwaysAssert(asap::nAxes==4,AipsError);
423 if (nPol()!=4) {
424 throw (AipsError("You must have 4 polarizations to get RR or LL"));
425 }
426 Array<Float> arr;
427 stokesCol_.get(whichRow, arr);
428
429 // Set current cursor location
430
431 const IPosition& shape = arr.shape();
432 IPosition start, end;
433 setCursorSlice(start, end, shape);
434
435 // Get I and V slices
436
437 Array<Float> I = SDPolUtil::getStokesSlice(arr,start,end,"I");
438 Array<Float> V = SDPolUtil::getStokesSlice(arr,start,end,"V");
439
440 // Compute output
441
442 Array<Float> out = SDPolUtil::circularPolarizationFromStokes(I, V, doRR);
443
444 // Copy to output
445
446 IPosition vecShape(1,shape(asap::ChanAxis));
447 Vector<Float> outV = out.reform(vecShape);
448 std::vector<float> stlout;
449 outV.tovector(stlout);
450
451 return stlout;
452}
453
454
455std::vector<string> SDMemTable::getCoordInfo() const
456{
457 String un;
458 Table t = table_.keywordSet().asTable("FREQUENCIES");
459 String sunit;
460 t.keywordSet().get("UNIT",sunit);
461 String dpl;
462 t.keywordSet().get("DOPPLER",dpl);
463 if (dpl == "") dpl = "RADIO";
464 String rfrm;
465 t.keywordSet().get("REFFRAME",rfrm);
466 std::vector<string> inf;
467 inf.push_back(sunit);
468 inf.push_back(rfrm);
469 inf.push_back(dpl);
470 t.keywordSet().get("BASEREFFRAME",rfrm);
471 inf.push_back(rfrm);
472 return inf;
473}
474
475void SDMemTable::setCoordInfo(std::vector<string> theinfo)
476{
477 std::vector<string>::iterator it;
478 String un,rfrm, brfrm,dpl;
479 un = theinfo[0]; // Abcissa unit
480 rfrm = theinfo[1]; // Active (or conversion) frame
481 dpl = theinfo[2]; // Doppler
482 brfrm = theinfo[3]; // Base frame
483
484 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
485
486 Vector<Double> rstf;
487 t.keywordSet().get("RESTFREQS",rstf);
488
489 Bool canDo = True;
490 Unit u1("km/s");Unit u2("Hz");
491 if (Unit(un) == u1) {
492 Vector<Double> rstf;
493 t.keywordSet().get("RESTFREQS",rstf);
494 if (rstf.nelements() == 0) {
495 throw(AipsError("Can't set unit to km/s if no restfrequencies are specified"));
496 }
497 } else if (Unit(un) != u2 && un != "") {
498 throw(AipsError("Unit not conformant with Spectral Coordinates"));
499 }
500 t.rwKeywordSet().define("UNIT", un);
501
502 MFrequency::Types mdr;
503 if (!MFrequency::getType(mdr, rfrm)) {
504
505 Int a,b;const uInt* c;
506 const String* valid = MFrequency::allMyTypes(a, b, c);
507 String pfix = "Please specify a legal frame type. Types are\n";
508 throw(AipsError(pfix+(*valid)));
509 } else {
510 t.rwKeywordSet().define("REFFRAME",rfrm);
511 }
512
513 MDoppler::Types dtype;
514 dpl.upcase();
515 if (!MDoppler::getType(dtype, dpl)) {
516 throw(AipsError("Doppler type unknown"));
517 } else {
518 t.rwKeywordSet().define("DOPPLER",dpl);
519 }
520
521 if (!MFrequency::getType(mdr, brfrm)) {
522 Int a,b;const uInt* c;
523 const String* valid = MFrequency::allMyTypes(a, b, c);
524 String pfix = "Please specify a legal frame type. Types are\n";
525 throw(AipsError(pfix+(*valid)));
526 } else {
527 t.rwKeywordSet().define("BASEREFFRAME",brfrm);
528 }
529}
530
531
532std::vector<double> SDMemTable::getAbcissa(Int whichRow) const
533{
534 std::vector<double> abc(nChan());
535
536 // Get header units keyword
537 Table t = table_.keywordSet().asTable("FREQUENCIES");
538 String sunit;
539 t.keywordSet().get("UNIT",sunit);
540 if (sunit == "") sunit = "pixel";
541 Unit u(sunit);
542
543 // Easy if just wanting pixels
544 if (sunit==String("pixel")) {
545 // assume channels/pixels
546 std::vector<double>::iterator it;
547 uInt i=0;
548 for (it = abc.begin(); it != abc.end(); ++it) {
549 (*it) = Double(i++);
550 }
551 return abc;
552 }
553
554 // Continue with km/s or Hz. Get FreqIDs
555 Vector<uInt> freqIDs;
556 freqidCol_.get(whichRow, freqIDs);
557 uInt freqID = freqIDs(IFSel_);
558 restfreqidCol_.get(whichRow, freqIDs);
559 uInt restFreqID = freqIDs(IFSel_);
560
561 // Get SpectralCoordinate, set reference frame conversion,
562 // velocity conversion, and rest freq state
563
564 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
565 Vector<Double> pixel(nChan());
566 indgen(pixel);
567
568 if (u == Unit("km/s")) {
569 Vector<Double> world;
570 spc.pixelToVelocity(world,pixel);
571 std::vector<double>::iterator it;
572 uInt i = 0;
573 for (it = abc.begin(); it != abc.end(); ++it) {
574 (*it) = world[i];
575 i++;
576 }
577 } else if (u == Unit("Hz")) {
578
579 // Set world axis units
580 Vector<String> wau(1); wau = u.getName();
581 spc.setWorldAxisUnits(wau);
582
583 std::vector<double>::iterator it;
584 Double tmp;
585 uInt i = 0;
586 for (it = abc.begin(); it != abc.end(); ++it) {
587 spc.toWorld(tmp,pixel[i]);
588 (*it) = tmp;
589 i++;
590 }
591 }
592 return abc;
593}
594
595std::string SDMemTable::getAbcissaString(Int whichRow) const
596{
597 Table t = table_.keywordSet().asTable("FREQUENCIES");
598
599 String sunit;
600 t.keywordSet().get("UNIT",sunit);
601 if (sunit == "") sunit = "pixel";
602 Unit u(sunit);
603
604 Vector<uInt> freqIDs;
605 freqidCol_.get(whichRow, freqIDs);
606 uInt freqID = freqIDs(IFSel_);
607 restfreqidCol_.get(whichRow, freqIDs);
608 uInt restFreqID = freqIDs(IFSel_);
609
610 // Get SpectralCoordinate, with frame, velocity, rest freq state set
611 SpectralCoordinate spc = getSpectralCoordinate(freqID, restFreqID, whichRow);
612
613 String s = "Channel";
614 if (u == Unit("km/s")) {
615 s = CoordinateUtil::axisLabel(spc,0,True,True,True);
616 } else if (u == Unit("Hz")) {
617 Vector<String> wau(1);wau = u.getName();
618 spc.setWorldAxisUnits(wau);
619
620 s = CoordinateUtil::axisLabel(spc,0,True,True,False);
621 }
622 return s;
623}
624
625void SDMemTable::setSpectrum(std::vector<float> spectrum, int whichRow)
626{
627 Array<Float> arr;
628 specCol_.get(whichRow, arr);
629 if (spectrum.size() != arr.shape()(asap::ChanAxis)) {
630 throw(AipsError("Attempting to set spectrum with incorrect length."));
631 }
632
633 // Setup accessors
634 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
635 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection
636 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
637 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection
638 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
639 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection
640
641 // Iterate
642 std::vector<float>::iterator it = spectrum.begin();
643 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
644 (*i) = Float(*it);
645 it++;
646 }
647 specCol_.put(whichRow, arr);
648}
649
650void SDMemTable::getSpectrum(Vector<Float>& spectrum, Int whichRow) const
651{
652 Array<Float> arr;
653 specCol_.get(whichRow, arr);
654
655 // Iterate and extract
656 spectrum.resize(arr.shape()(3));
657 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
658 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
659 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
660 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
661 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
662 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
663
664 ArrayAccessor<Float, Axis<asap::BeamAxis> > va(spectrum);
665 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
666 (*va) = (*i);
667 va++;
668 }
669}
670
671
672/*
673void SDMemTable::getMask(Vector<Bool>& mask, Int whichRow) const {
674 Array<uChar> arr;
675 flagsCol_.get(whichRow, arr);
676 mask.resize(arr.shape()(3));
677
678 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
679 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
680 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
681 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
682 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
683 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
684
685 Bool useUserMask = ( chanMask_.size() == arr.shape()(3) );
686
687 ArrayAccessor<Bool, Axis<asap::BeamAxis> > va(mask);
688 std::vector<bool> tmp;
689 tmp = chanMask_; // WHY the fxxx do I have to make a copy here. The
690 // iterator should work on chanMask_??
691 std::vector<bool>::iterator miter;
692 miter = tmp.begin();
693
694 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
695 bool out =!static_cast<bool>(*i);
696 if (useUserMask) {
697 out = out && (*miter);
698 miter++;
699 }
700 (*va) = out;
701 va++;
702 }
703}
704*/
705
706MaskedArray<Float> SDMemTable::rowAsMaskedArray(uInt whichRow,
707 Bool toStokes) const
708{
709 // Get flags
710 Array<uChar> farr;
711 flagsCol_.get(whichRow, farr);
712
713 // Get data and convert mask to Bool
714 Array<Float> arr;
715 Array<Bool> mask;
716 if (toStokes) {
717 stokesCol_.get(whichRow, arr);
718
719 Array<Bool> tMask(farr.shape());
720 convertArray(tMask, farr);
721 mask = SDPolUtil::stokesData (tMask, True);
722 } else {
723 specCol_.get(whichRow, arr);
724 mask.resize(farr.shape());
725 convertArray(mask, farr);
726 }
727
728 return MaskedArray<Float>(arr,!mask);
729}
730
731Float SDMemTable::getTsys(Int whichRow) const
732{
733 Array<Float> arr;
734 tsCol_.get(whichRow, arr);
735 Float out;
736
737 IPosition ip(arr.shape());
738 ip(asap::BeamAxis) = beamSel_;
739 ip(asap::IFAxis) = IFSel_;
740 ip(asap::PolAxis) = polSel_;
741 ip(asap::ChanAxis)=0; // First channel only
742
743 out = arr(ip);
744 return out;
745}
746
747MDirection SDMemTable::getDirection(Int whichRow, Bool refBeam) const
748{
749 MDirection::Types mdr = getDirectionReference();
750 Array<Double> posit;
751 dirCol_.get(whichRow,posit);
752 Vector<Double> wpos(2);
753 Int rb;
754 rbeamCol_.get(whichRow,rb);
755 wpos[0] = posit(IPosition(2,beamSel_,0));
756 wpos[1] = posit(IPosition(2,beamSel_,1));
757 if (refBeam && rb > -1) { // use refBeam instead if it exists
758 wpos[0] = posit(IPosition(2,rb,0));
759 wpos[1] = posit(IPosition(2,rb,1));
760 }
761
762 Quantum<Double> lon(wpos[0],Unit(String("rad")));
763 Quantum<Double> lat(wpos[1],Unit(String("rad")));
764 return MDirection(lon, lat, mdr);
765}
766
767MEpoch SDMemTable::getEpoch(Int whichRow) const
768{
769 MEpoch::Types met = getTimeReference();
770
771 Double obstime;
772 timeCol_.get(whichRow,obstime);
773 MVEpoch tm2(Quantum<Double>(obstime, Unit(String("d"))));
774 return MEpoch(tm2, met);
775}
776
777MPosition SDMemTable::getAntennaPosition () const
778{
779 Vector<Double> antpos;
780 table_.keywordSet().get("AntennaPosition", antpos);
781 MVPosition mvpos(antpos(0),antpos(1),antpos(2));
782 return MPosition(mvpos);
783}
784
785
786SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID) const
787{
788
789 Table t = table_.keywordSet().asTable("FREQUENCIES");
790 if (freqID> t.nrow() ) {
791 throw(AipsError("SDMemTable::getSpectralCoordinate - freqID out of range"));
792 }
793
794 Double rp,rv,inc;
795 String rf;
796 ROScalarColumn<Double> rpc(t, "REFPIX");
797 ROScalarColumn<Double> rvc(t, "REFVAL");
798 ROScalarColumn<Double> incc(t, "INCREMENT");
799 t.keywordSet().get("BASEREFFRAME",rf);
800
801 // Create SpectralCoordinate (units Hz)
802 MFrequency::Types mft;
803 if (!MFrequency::getType(mft, rf)) {
804 cerr << "Frequency type unknown assuming TOPO" << endl;
805 mft = MFrequency::TOPO;
806 }
807 rpc.get(freqID, rp);
808 rvc.get(freqID, rv);
809 incc.get(freqID, inc);
810
811 SpectralCoordinate spec(mft,rv,inc,rp);
812 return spec;
813}
814
815
816SpectralCoordinate SDMemTable::getSpectralCoordinate(uInt freqID,
817 uInt restFreqID,
818 uInt whichRow) const
819{
820
821 // Create basic SC
822 SpectralCoordinate spec = getSpectralCoordinate (freqID);
823
824 Table t = table_.keywordSet().asTable("FREQUENCIES");
825
826 // Set rest frequency
827 Vector<Double> restFreqIDs;
828 t.keywordSet().get("RESTFREQS",restFreqIDs);
829 if (restFreqID < restFreqIDs.nelements()) { // Should always be true
830 spec.setRestFrequency(restFreqIDs(restFreqID),True);
831 }
832
833 // Set up frame conversion layer
834 String frm;
835 t.keywordSet().get("REFFRAME",frm);
836 if (frm == "") frm = "TOPO";
837 MFrequency::Types mtype;
838 if (!MFrequency::getType(mtype, frm)) {
839 // Should never happen
840 cerr << "Frequency type unknown assuming TOPO" << endl;
841 mtype = MFrequency::TOPO;
842 }
843
844 // Set reference frame conversion (requires row)
845 MDirection dir = getDirection(whichRow);
846 MEpoch epoch = getEpoch(whichRow);
847 MPosition pos = getAntennaPosition();
848
849 if (!spec.setReferenceConversion(mtype,epoch,pos,dir)) {
850 throw(AipsError("Couldn't convert frequency frame."));
851 }
852
853 // Now velocity conversion if appropriate
854 String unitStr;
855 t.keywordSet().get("UNIT",unitStr);
856
857 String dpl;
858 t.keywordSet().get("DOPPLER",dpl);
859 MDoppler::Types dtype;
860 MDoppler::getType(dtype, dpl);
861
862 // Only set velocity unit if non-blank and non-Hz
863 if (!unitStr.empty()) {
864 Unit unitU(unitStr);
865 if (unitU==Unit("Hz")) {
866 } else {
867 spec.setVelocity(unitStr, dtype);
868 }
869 }
870
871 return spec;
872}
873
874
875Bool SDMemTable::setCoordinate(const SpectralCoordinate& speccord,
876 uInt freqID) {
877 Table t = table_.rwKeywordSet().asTable("FREQUENCIES");
878 if (freqID > t.nrow() ) {
879 throw(AipsError("SDMemTable::setCoordinate - coord no out of range"));
880 }
881 ScalarColumn<Double> rpc(t, "REFPIX");
882 ScalarColumn<Double> rvc(t, "REFVAL");
883 ScalarColumn<Double> incc(t, "INCREMENT");
884
885 rpc.put(freqID, speccord.referencePixel()[0]);
886 rvc.put(freqID, speccord.referenceValue()[0]);
887 incc.put(freqID, speccord.increment()[0]);
888
889 return True;
890}
891
892Int SDMemTable::nCoordinates() const
893{
894 return table_.keywordSet().asTable("FREQUENCIES").nrow();
895}
896
897
898std::vector<double> SDMemTable::getRestFreqs() const
899{
900 Table t = table_.keywordSet().asTable("FREQUENCIES");
901 Vector<Double> tvec;
902 t.keywordSet().get("RESTFREQS",tvec);
903 std::vector<double> stlout;
904 tvec.tovector(stlout);
905 return stlout;
906}
907
908bool SDMemTable::putSDFreqTable(const SDFrequencyTable& sdft)
909{
910 TableDesc td("", "1", TableDesc::Scratch);
911 td.addColumn(ScalarColumnDesc<Double>("REFPIX"));
912 td.addColumn(ScalarColumnDesc<Double>("REFVAL"));
913 td.addColumn(ScalarColumnDesc<Double>("INCREMENT"));
914 SetupNewTable aNewTab("freqs", td, Table::New);
915 Table aTable (aNewTab, Table::Memory, sdft.length());
916 ScalarColumn<Double> sc0(aTable, "REFPIX");
917 ScalarColumn<Double> sc1(aTable, "REFVAL");
918 ScalarColumn<Double> sc2(aTable, "INCREMENT");
919 for (uInt i=0; i < sdft.length(); ++i) {
920 sc0.put(i,sdft.referencePixel(i));
921 sc1.put(i,sdft.referenceValue(i));
922 sc2.put(i,sdft.increment(i));
923 }
924 String rf = sdft.refFrame();
925 if (rf.contains("TOPO")) rf = "TOPO";
926
927 aTable.rwKeywordSet().define("BASEREFFRAME", rf);
928 aTable.rwKeywordSet().define("REFFRAME", rf);
929 aTable.rwKeywordSet().define("EQUINOX", sdft.equinox());
930 aTable.rwKeywordSet().define("UNIT", String(""));
931 aTable.rwKeywordSet().define("DOPPLER", String("RADIO"));
932 Vector<Double> rfvec;
933 String rfunit;
934 sdft.restFrequencies(rfvec,rfunit);
935 Quantum<Vector<Double> > q(rfvec, rfunit);
936 rfvec.resize();
937 rfvec = q.getValue("Hz");
938 aTable.rwKeywordSet().define("RESTFREQS", rfvec);
939 table_.rwKeywordSet().defineTable("FREQUENCIES", aTable);
940 return true;
941}
942
943bool SDMemTable::putSDFitTable(const SDFitTable& sdft)
944{
945 TableDesc td("", "1", TableDesc::Scratch);
946 td.addColumn(ArrayColumnDesc<String>("FUNCTIONS"));
947 td.addColumn(ArrayColumnDesc<Int>("COMPONENTS"));
948 td.addColumn(ArrayColumnDesc<Double>("PARAMETERS"));
949 td.addColumn(ArrayColumnDesc<Bool>("PARMASK"));
950 td.addColumn(ArrayColumnDesc<String>("FRAMEINFO"));
951 SetupNewTable aNewTab("fits", td, Table::New);
952 Table aTable(aNewTab, Table::Memory);
953 ArrayColumn<String> sc0(aTable, "FUNCTIONS");
954 ArrayColumn<Int> sc1(aTable, "COMPONENTS");
955 ArrayColumn<Double> sc2(aTable, "PARAMETERS");
956 ArrayColumn<Bool> sc3(aTable, "PARMASK");
957 ArrayColumn<String> sc4(aTable, "FRAMEINFO");
958 for (uInt i; i<sdft.length(); ++i) {
959 const Vector<Double>& parms = sdft.getParameters(i);
960 const Vector<Bool>& parmask = sdft.getParameterMask(i);
961 const Vector<String>& funcs = sdft.getFunctions(i);
962 const Vector<Int>& comps = sdft.getComponents(i);
963 const Vector<String>& finfo = sdft.getFrameInfo(i);
964 sc0.put(i,funcs);
965 sc1.put(i,comps);
966 sc3.put(i,parmask);
967 sc2.put(i,parms);
968 sc4.put(i,finfo);
969 }
970 table_.rwKeywordSet().defineTable("FITS", aTable);
971 return true;
972}
973
974SDFitTable SDMemTable::getSDFitTable() const
975{
976 const Table& t = table_.keywordSet().asTable("FITS");
977 Vector<Double> parms;
978 Vector<Bool> parmask;
979 Vector<String> funcs;
980 Vector<String> finfo;
981 Vector<Int> comps;
982 ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
983 ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
984 ROArrayColumn<Int> compsCol(t, "COMPONENTS");
985 ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
986 ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
987 uInt n = t.nrow();
988 SDFitTable sdft;
989 for (uInt i=0; i<n; ++i) {
990 parmaskCol.get(i, parmask);
991 parmsCol.get(i, parms);
992 funcsCol.get(i, funcs);
993 compsCol.get(i, comps);
994 finfoCol.get(i, finfo);
995 sdft.addFit(parms, parmask, funcs, comps, finfo);
996 }
997 return sdft;
998}
999
1000SDFitTable SDMemTable::getSDFitTable(uInt whichRow) const {
1001 Array<Int> fitid;
1002 fitCol_.get(whichRow, fitid);
1003 if (fitid.nelements() == 0) return SDFitTable();
1004
1005 IPosition shp = fitid.shape();
1006 IPosition start(4, beamSel_, IFSel_, polSel_,0);
1007 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
1008
1009 // reform the output array slice to be of dim=1
1010 Vector<Int> tmp = (fitid(start, end)).reform(IPosition(1,shp[3]));
1011
1012 const Table& t = table_.keywordSet().asTable("FITS");
1013 Vector<Double> parms;
1014 Vector<Bool> parmask;
1015 Vector<String> funcs;
1016 Vector<String> finfo;
1017 Vector<Int> comps;
1018 ROArrayColumn<Double> parmsCol(t, "PARAMETERS");
1019 ROArrayColumn<Bool> parmaskCol(t, "PARMASK");
1020 ROArrayColumn<Int> compsCol(t, "COMPONENTS");
1021 ROArrayColumn<String> funcsCol(t, "FUNCTIONS");
1022 ROArrayColumn<String> finfoCol(t, "FRAMEINFO");
1023 if (t.nrow() == 0) return SDFitTable();
1024 SDFitTable sdft;
1025 Int k=-1;
1026 for (uInt i=0; i< tmp.nelements(); ++i) {
1027 k = tmp[i];
1028 if ( k > -1 && k < t.nrow() ) {
1029 parms.resize();
1030 parmsCol.get(k, parms);
1031 parmask.resize();
1032 parmaskCol.get(k, parmask);
1033 funcs.resize();
1034 funcsCol.get(k, funcs);
1035 comps.resize();
1036 compsCol.get(k, comps);
1037 finfo.resize();
1038 finfoCol.get(k, finfo);
1039 sdft.addFit(parms, parmask, funcs, comps, finfo);
1040 }
1041 }
1042 return sdft;
1043}
1044
1045void SDMemTable::addFit(uInt whichRow,
1046 const Vector<Double>& p, const Vector<Bool>& m,
1047 const Vector<String>& f, const Vector<Int>& c)
1048{
1049 if (whichRow >= nRow()) {
1050 throw(AipsError("Specified row out of range"));
1051 }
1052 Table t = table_.keywordSet().asTable("FITS");
1053 uInt nrow = t.nrow();
1054 t.addRow();
1055 ArrayColumn<Double> parmsCol(t, "PARAMETERS");
1056 ArrayColumn<Bool> parmaskCol(t, "PARMASK");
1057 ArrayColumn<Int> compsCol(t, "COMPONENTS");
1058 ArrayColumn<String> funcsCol(t, "FUNCTIONS");
1059 ArrayColumn<String> finfoCol(t, "FRAMEINFO");
1060 parmsCol.put(nrow, p);
1061 parmaskCol.put(nrow, m);
1062 compsCol.put(nrow, c);
1063 funcsCol.put(nrow, f);
1064 Vector<String> fi = mathutil::toVectorString(getCoordInfo());
1065 finfoCol.put(nrow, fi);
1066
1067 Array<Int> fitarr;
1068 fitCol_.get(whichRow, fitarr);
1069
1070 Array<Int> newarr; // The new Array containing the fitid
1071 Int pos =-1; // The fitid position in the array
1072 if ( fitarr.nelements() == 0 ) { // no fits at all in this row
1073 Array<Int> arr(IPosition(4,nBeam(),nIF(),nPol(),1));
1074 arr = -1;
1075 newarr.reference(arr);
1076 pos = 0;
1077 } else {
1078 IPosition shp = fitarr.shape();
1079 IPosition start(4, beamSel_, IFSel_, polSel_,0);
1080 IPosition end(4, beamSel_, IFSel_, polSel_, shp[3]-1);
1081 // reform the output array slice to be of dim=1
1082 Array<Int> tmp = (fitarr(start, end)).reform(IPosition(1,shp[3]));
1083 const Vector<Int>& fits = tmp;
1084 VectorSTLIterator<Int> it(fits);
1085 Int i = 0;
1086 while (it != fits.end()) {
1087 if (*it == -1) {
1088 pos = i;
1089 break;
1090 }
1091 ++i;
1092 ++it;
1093 };
1094 }
1095 if (pos == -1) {
1096 mathutil::extendLastArrayAxis(newarr,fitarr, -1);
1097 pos = fitarr.shape()[3]; // the new element position
1098 } else {
1099 if (fitarr.nelements() > 0)
1100 newarr = fitarr;
1101 }
1102 newarr(IPosition(4, beamSel_, IFSel_, polSel_, pos)) = Int(nrow);
1103 fitCol_.put(whichRow, newarr);
1104
1105}
1106
1107SDFrequencyTable SDMemTable::getSDFreqTable() const
1108{
1109 const Table& t = table_.keywordSet().asTable("FREQUENCIES");
1110 SDFrequencyTable sdft;
1111
1112 // Add refpix/refval/incr. What are the units ? Hz I suppose
1113 // but it's nowhere described...
1114 Vector<Double> refPix, refVal, incr;
1115 ScalarColumn<Double> refPixCol(t, "REFPIX");
1116 ScalarColumn<Double> refValCol(t, "REFVAL");
1117 ScalarColumn<Double> incrCol(t, "INCREMENT");
1118 refPix = refPixCol.getColumn();
1119 refVal = refValCol.getColumn();
1120 incr = incrCol.getColumn();
1121
1122 uInt n = refPix.nelements();
1123 for (uInt i=0; i<n; i++) {
1124 sdft.addFrequency(refPix[i], refVal[i], incr[i]);
1125 }
1126
1127 // Frequency reference frame. I don't know if this
1128 // is the correct frame. It might be 'REFFRAME'
1129 // rather than 'BASEREFFRAME' ?
1130 String baseFrame;
1131 t.keywordSet().get("BASEREFFRAME",baseFrame);
1132 sdft.setRefFrame(baseFrame);
1133
1134 // Equinox
1135 Float equinox;
1136 t.keywordSet().get("EQUINOX", equinox);
1137 sdft.setEquinox(equinox);
1138
1139 // Rest Frequency
1140 Vector<Double> restFreqs;
1141 t.keywordSet().get("RESTFREQS", restFreqs);
1142 for (uInt i=0; i<restFreqs.nelements(); i++) {
1143 sdft.addRestFrequency(restFreqs[i]);
1144 }
1145 sdft.setRestFrequencyUnit(String("Hz"));
1146
1147 return sdft;
1148}
1149
1150bool SDMemTable::putSDContainer(const SDContainer& sdc)
1151{
1152 uInt rno = table_.nrow();
1153 table_.addRow();
1154
1155 timeCol_.put(rno, sdc.timestamp);
1156 srcnCol_.put(rno, sdc.sourcename);
1157 fldnCol_.put(rno, sdc.fieldname);
1158 specCol_.put(rno, sdc.getSpectrum());
1159 flagsCol_.put(rno, sdc.getFlags());
1160 tsCol_.put(rno, sdc.getTsys());
1161 scanCol_.put(rno, sdc.scanid);
1162 integrCol_.put(rno, sdc.interval);
1163 freqidCol_.put(rno, sdc.getFreqMap());
1164 restfreqidCol_.put(rno, sdc.getRestFreqMap());
1165 dirCol_.put(rno, sdc.getDirection());
1166 rbeamCol_.put(rno, sdc.refbeam);
1167 tcalCol_.put(rno, sdc.tcal);
1168 tcaltCol_.put(rno, sdc.tcaltime);
1169 azCol_.put(rno, sdc.azimuth);
1170 elCol_.put(rno, sdc.elevation);
1171 paraCol_.put(rno, sdc.parangle);
1172 histCol_.put(rno, sdc.getHistory());
1173 fitCol_.put(rno, sdc.getFitMap());
1174
1175 return true;
1176}
1177
1178SDContainer SDMemTable::getSDContainer(uInt whichRow) const
1179{
1180 SDContainer sdc(nBeam(),nIF(),nPol(),nChan());
1181 timeCol_.get(whichRow, sdc.timestamp);
1182 srcnCol_.get(whichRow, sdc.sourcename);
1183 integrCol_.get(whichRow, sdc.interval);
1184 scanCol_.get(whichRow, sdc.scanid);
1185 fldnCol_.get(whichRow, sdc.fieldname);
1186 rbeamCol_.get(whichRow, sdc.refbeam);
1187 azCol_.get(whichRow, sdc.azimuth);
1188 elCol_.get(whichRow, sdc.elevation);
1189 paraCol_.get(whichRow, sdc.parangle);
1190 Vector<Float> tc;
1191 tcalCol_.get(whichRow, tc);
1192 sdc.tcal[0] = tc[0];sdc.tcal[1] = tc[1];
1193 tcaltCol_.get(whichRow, sdc.tcaltime);
1194
1195 Array<Float> spectrum;
1196 Array<Float> tsys;
1197 Array<uChar> flagtrum;
1198 Vector<uInt> fmap;
1199 Array<Double> direction;
1200 Vector<String> histo;
1201 Array<Int> fits;
1202
1203 specCol_.get(whichRow, spectrum);
1204 sdc.putSpectrum(spectrum);
1205 flagsCol_.get(whichRow, flagtrum);
1206 sdc.putFlags(flagtrum);
1207 tsCol_.get(whichRow, tsys);
1208 sdc.putTsys(tsys);
1209 freqidCol_.get(whichRow, fmap);
1210 sdc.putFreqMap(fmap);
1211 restfreqidCol_.get(whichRow, fmap);
1212 sdc.putRestFreqMap(fmap);
1213 dirCol_.get(whichRow, direction);
1214 sdc.putDirection(direction);
1215 histCol_.get(whichRow, histo);
1216 sdc.putHistory(histo);
1217 fitCol_.get(whichRow, fits);
1218 sdc.putFitMap(fits);
1219 return sdc;
1220}
1221
1222bool SDMemTable::putSDHeader(const SDHeader& sdh)
1223{
1224 table_.rwKeywordSet().define("nIF", sdh.nif);
1225 table_.rwKeywordSet().define("nBeam", sdh.nbeam);
1226 table_.rwKeywordSet().define("nPol", sdh.npol);
1227 table_.rwKeywordSet().define("nChan", sdh.nchan);
1228 table_.rwKeywordSet().define("Observer", sdh.observer);
1229 table_.rwKeywordSet().define("Project", sdh.project);
1230 table_.rwKeywordSet().define("Obstype", sdh.obstype);
1231 table_.rwKeywordSet().define("AntennaName", sdh.antennaname);
1232 table_.rwKeywordSet().define("AntennaPosition", sdh.antennaposition);
1233 table_.rwKeywordSet().define("Equinox", sdh.equinox);
1234 table_.rwKeywordSet().define("FreqRefFrame", sdh.freqref);
1235 table_.rwKeywordSet().define("FreqRefVal", sdh.reffreq);
1236 table_.rwKeywordSet().define("Bandwidth", sdh.bandwidth);
1237 table_.rwKeywordSet().define("UTC", sdh.utc);
1238 table_.rwKeywordSet().define("FluxUnit", sdh.fluxunit);
1239 table_.rwKeywordSet().define("Epoch", sdh.epoch);
1240 return true;
1241}
1242
1243SDHeader SDMemTable::getSDHeader() const
1244{
1245 SDHeader sdh;
1246 table_.keywordSet().get("nBeam",sdh.nbeam);
1247 table_.keywordSet().get("nIF",sdh.nif);
1248 table_.keywordSet().get("nPol",sdh.npol);
1249 table_.keywordSet().get("nChan",sdh.nchan);
1250 table_.keywordSet().get("Observer", sdh.observer);
1251 table_.keywordSet().get("Project", sdh.project);
1252 table_.keywordSet().get("Obstype", sdh.obstype);
1253 table_.keywordSet().get("AntennaName", sdh.antennaname);
1254 table_.keywordSet().get("AntennaPosition", sdh.antennaposition);
1255 table_.keywordSet().get("Equinox", sdh.equinox);
1256 table_.keywordSet().get("FreqRefFrame", sdh.freqref);
1257 table_.keywordSet().get("FreqRefVal", sdh.reffreq);
1258 table_.keywordSet().get("Bandwidth", sdh.bandwidth);
1259 table_.keywordSet().get("UTC", sdh.utc);
1260 table_.keywordSet().get("FluxUnit", sdh.fluxunit);
1261 table_.keywordSet().get("Epoch", sdh.epoch);
1262 return sdh;
1263}
1264void SDMemTable::makePersistent(const std::string& filename)
1265{
1266 table_.deepCopy(filename,Table::New);
1267}
1268
1269Int SDMemTable::nScan() const {
1270 Int n = 0;
1271 Int previous = -1;Int current=0;
1272 for (uInt i=0; i< scanCol_.nrow();i++) {
1273 scanCol_.getScalar(i,current);
1274 if (previous != current) {
1275 previous = current;
1276 n++;
1277 }
1278 }
1279 return n;
1280}
1281
1282String SDMemTable::formatSec(Double x) const
1283{
1284 Double xcop = x;
1285 MVTime mvt(xcop/24./3600.); // make days
1286
1287 if (x < 59.95)
1288 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_HM, 7)+"s";
1289 else if (x < 3599.95)
1290 return String(" ") + mvt.string(MVTime::TIME_CLEAN_NO_H,7)+" ";
1291 else {
1292 ostringstream oss;
1293 oss << setw(2) << std::right << setprecision(1) << mvt.hour();
1294 oss << ":" << mvt.string(MVTime::TIME_CLEAN_NO_H,7) << " ";
1295 return String(oss);
1296 }
1297};
1298
1299String SDMemTable::formatDirection(const MDirection& md) const
1300{
1301 Vector<Double> t = md.getAngle(Unit(String("rad"))).getValue();
1302 Int prec = 7;
1303
1304 MVAngle mvLon(t[0]);
1305 String sLon = mvLon.string(MVAngle::TIME,prec);
1306 MVAngle mvLat(t[1]);
1307 String sLat = mvLat.string(MVAngle::ANGLE+MVAngle::DIG2,prec);
1308 return sLon + String(" ") + sLat;
1309}
1310
1311
1312std::string SDMemTable::getFluxUnit() const
1313{
1314 String tmp;
1315 table_.keywordSet().get("FluxUnit", tmp);
1316 return tmp;
1317}
1318
1319void SDMemTable::setFluxUnit(const std::string& unit)
1320{
1321 String tmp(unit);
1322 Unit tU(tmp);
1323 if (tU==Unit("K") || tU==Unit("Jy")) {
1324 table_.rwKeywordSet().define(String("FluxUnit"), tmp);
1325 } else {
1326 throw AipsError("Illegal unit - must be compatible with Jy or K");
1327 }
1328}
1329
1330
1331void SDMemTable::setInstrument(const std::string& name)
1332{
1333 Bool throwIt = True;
1334 Instrument ins = convertInstrument(name, throwIt);
1335 String nameU(name);
1336 nameU.upcase();
1337 table_.rwKeywordSet().define(String("AntennaName"), nameU);
1338}
1339
1340std::string SDMemTable::summary(bool verbose) const {
1341
1342 // Format header info
1343 ostringstream oss;
1344 oss << endl;
1345 oss << "--------------------------------------------------------------------------------" << endl;
1346 oss << " Scan Table Summary" << endl;
1347 oss << "--------------------------------------------------------------------------------" << endl;
1348 oss.flags(std::ios_base::left);
1349 oss << setw(15) << "Beams:" << setw(4) << nBeam() << endl
1350 << setw(15) << "IFs:" << setw(4) << nIF() << endl
1351 << setw(15) << "Polarisations:" << setw(4) << nPol() << endl
1352 << setw(15) << "Channels:" << setw(4) << nChan() << endl;
1353 oss << endl;
1354 String tmp;
1355 table_.keywordSet().get("Observer", tmp);
1356 oss << setw(15) << "Observer:" << tmp << endl;
1357 oss << setw(15) << "Obs Date:" << getTime(-1,True) << endl;
1358 table_.keywordSet().get("Project", tmp);
1359 oss << setw(15) << "Project:" << tmp << endl;
1360 table_.keywordSet().get("Obstype", tmp);
1361 oss << setw(15) << "Obs. Type:" << tmp << endl;
1362 table_.keywordSet().get("AntennaName", tmp);
1363 oss << setw(15) << "Antenna Name:" << tmp << endl;
1364 table_.keywordSet().get("FluxUnit", tmp);
1365 oss << setw(15) << "Flux Unit:" << tmp << endl;
1366 Table t = table_.keywordSet().asTable("FREQUENCIES");
1367 Vector<Double> vec;
1368 t.keywordSet().get("RESTFREQS",vec);
1369 oss << setw(15) << "Rest Freqs:";
1370 if (vec.nelements() > 0) {
1371 oss << setprecision(0) << vec << " [Hz]" << endl;
1372 } else {
1373 oss << "None set" << endl;
1374 }
1375 oss << setw(15) << "Abcissa:" << getAbcissaString() << endl;
1376 oss << setw(15) << "Cursor:" << "Beam[" << getBeam() << "] "
1377 << "IF[" << getIF() << "] " << "Pol[" << getPol() << "]" << endl;
1378 oss << endl;
1379
1380 String dirtype ="Position ("+ MDirection::showType(getDirectionReference()) + ")";
1381 oss << setw(5) << "Scan"
1382 << setw(15) << "Source"
1383 << setw(24) << dirtype
1384 << setw(10) << "Time"
1385 << setw(18) << "Integration"
1386 << setw(7) << "FreqIDs" << endl;
1387 oss << "--------------------------------------------------------------------------------" << endl;
1388
1389 // Generate list of scan start and end integrations
1390 Vector<Int> scanIDs = scanCol_.getColumn();
1391 Vector<uInt> startInt, endInt;
1392 mathutil::scanBoundaries(startInt, endInt, scanIDs);
1393
1394 const uInt nScans = startInt.nelements();
1395 String name;
1396 Vector<uInt> freqIDs, listFQ;
1397 uInt nInt;
1398
1399 for (uInt i=0; i<nScans; i++) {
1400 // Get things from first integration of scan
1401 String time = getTime(startInt(i),False);
1402 String tInt = formatSec(Double(getInterval(startInt(i))));
1403 String posit = formatDirection(getDirection(startInt(i),True));
1404 srcnCol_.getScalar(startInt(i),name);
1405
1406 // Find all the FreqIDs in this scan
1407 listFQ.resize(0);
1408 for (uInt j=startInt(i); j<endInt(i)+1; j++) {
1409 freqidCol_.get(j, freqIDs);
1410 for (uInt k=0; k<freqIDs.nelements(); k++) {
1411 mathutil::addEntry(listFQ, freqIDs(k));
1412 }
1413 }
1414
1415 nInt = endInt(i) - startInt(i) + 1;
1416 oss << setw(3) << std::right << i << std::left << setw(2) << " "
1417 << setw(15) << name
1418 << setw(24) << posit
1419 << setw(10) << time
1420 << setw(3) << std::right << nInt << setw(3) << " x " << std::left
1421 << setw(6) << tInt
1422 << " " << listFQ << endl;
1423 }
1424 oss << endl;
1425 oss << "Table contains " << table_.nrow() << " integration(s) in "
1426 << nScans << " scan(s)." << endl;
1427
1428 // Frequency Table
1429 if (verbose) {
1430 std::vector<string> info = getCoordInfo();
1431 SDFrequencyTable sdft = getSDFreqTable();
1432 oss << endl << endl;
1433 oss << "FreqID Frame RefFreq(Hz) RefPix Increment(Hz)" << endl;
1434 oss << "--------------------------------------------------------------------------------" << endl;
1435 for (uInt i=0; i<sdft.length(); i++) {
1436 oss << setw(8) << i << setw(8)
1437 << info[3] << setw(16) << setprecision(8)
1438 << sdft.referenceValue(i) << setw(10)
1439 << sdft.referencePixel(i) << setw(12)
1440 << sdft.increment(i) << endl;
1441 }
1442 oss << "--------------------------------------------------------------------------------" << endl;
1443 }
1444 return String(oss);
1445}
1446
1447Int SDMemTable::nBeam() const
1448{
1449 Int n;
1450 table_.keywordSet().get("nBeam",n);
1451 return n;
1452}
1453
1454Int SDMemTable::nIF() const {
1455 Int n;
1456 table_.keywordSet().get("nIF",n);
1457 return n;
1458}
1459
1460Int SDMemTable::nPol() const {
1461 Int n;
1462 table_.keywordSet().get("nPol",n);
1463 return n;
1464}
1465
1466Int SDMemTable::nChan() const {
1467 Int n;
1468 table_.keywordSet().get("nChan",n);
1469 return n;
1470}
1471
1472bool SDMemTable::appendHistory(const std::string& hist, int whichRow)
1473{
1474 Vector<String> history;
1475 histCol_.get(whichRow, history);
1476 history.resize(history.nelements()+1,True);
1477 history[history.nelements()-1] = hist;
1478 histCol_.put(whichRow, history);
1479}
1480
1481std::vector<std::string> SDMemTable::history(int whichRow) const
1482{
1483 Vector<String> history;
1484 histCol_.get(whichRow, history);
1485 std::vector<std::string> stlout = mathutil::tovectorstring(history);
1486 return stlout;
1487}
1488/*
1489void SDMemTable::maskChannels(const std::vector<Int>& whichChans ) {
1490
1491 std::vector<int>::iterator it;
1492 ArrayAccessor<uChar, Axis<asap::PolAxis> > j(flags_);
1493 for (it = whichChans.begin(); it != whichChans.end(); it++) {
1494 j.reset(j.begin(uInt(*it)));
1495 for (ArrayAccessor<uChar, Axis<asap::BeamAxis> > i(j); i != i.end(); ++i) {
1496 for (ArrayAccessor<uChar, Axis<asap::IFAxis> > ii(i); ii != ii.end(); ++ii) {
1497 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > iii(ii);
1498 iii != iii.end(); ++iii) {
1499 (*iii) =
1500 }
1501 }
1502 }
1503 }
1504
1505}
1506*/
1507void SDMemTable::flag(int whichRow)
1508{
1509 Array<uChar> arr;
1510 flagsCol_.get(whichRow, arr);
1511
1512 ArrayAccessor<uChar, Axis<asap::BeamAxis> > aa0(arr);
1513 aa0.reset(aa0.begin(uInt(beamSel_)));//go to beam
1514 ArrayAccessor<uChar, Axis<asap::IFAxis> > aa1(aa0);
1515 aa1.reset(aa1.begin(uInt(IFSel_)));// go to IF
1516 ArrayAccessor<uChar, Axis<asap::PolAxis> > aa2(aa1);
1517 aa2.reset(aa2.begin(uInt(polSel_)));// go to pol
1518
1519 for (ArrayAccessor<uChar, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1520 (*i) = uChar(True);
1521 }
1522
1523 flagsCol_.put(whichRow, arr);
1524}
1525
1526MDirection::Types SDMemTable::getDirectionReference() const
1527{
1528 Float eq;
1529 table_.keywordSet().get("Equinox",eq);
1530 std::map<float,string> mp;
1531 mp[2000.0] = "J2000";
1532 mp[1950.0] = "B1950";
1533 MDirection::Types mdr;
1534 if (!MDirection::getType(mdr, mp[eq])) {
1535 mdr = MDirection::J2000;
1536 cerr << "Unknown equinox using J2000" << endl;
1537 }
1538
1539 return mdr;
1540}
1541
1542MEpoch::Types SDMemTable::getTimeReference() const
1543{
1544 MEpoch::Types met;
1545 String ep;
1546 table_.keywordSet().get("Epoch",ep);
1547 if (!MEpoch::getType(met, ep)) {
1548 cerr << "Epoch type unknown - using UTC" << endl;
1549 met = MEpoch::UTC;
1550 }
1551
1552 return met;
1553}
1554
1555
1556Instrument SDMemTable::convertInstrument(const String& instrument,
1557 Bool throwIt)
1558{
1559 String t(instrument);
1560 t.upcase();
1561
1562 // The strings are what SDReader returns, after cunning interrogation
1563 // of station names... :-(
1564 Instrument inst = asap::UNKNOWN;
1565 if (t==String("DSS-43")) {
1566 inst = TIDBINBILLA;
1567 } else if (t==String("ATPKSMB")) {
1568 inst = ATPKSMB;
1569 } else if (t==String("ATPKSHOH")) {
1570 inst = ATPKSHOH;
1571 } else if (t==String("ATMOPRA")) {
1572 inst = ATMOPRA;
1573 } else if (t==String("CEDUNA")) {
1574 inst = CEDUNA;
1575 } else if (t==String("HOBART")) {
1576 inst = HOBART;
1577 } else {
1578 if (throwIt) {
1579 throw AipsError("Unrecognized instrument - use function scan.set_instrument to set");
1580 }
1581 }
1582 return inst;
1583}
1584
1585Bool SDMemTable::setRestFreqs(const Vector<Double>& restFreqsIn,
1586 const String& sUnit,
1587 const vector<string>& lines,
1588 const String& source,
1589 Int whichIF)
1590{
1591 const Int nIFs = nIF();
1592 if (whichIF>=nIFs) {
1593 throw(AipsError("Illegal IF index"));
1594 }
1595
1596 // Find vector of restfrequencies
1597 // Double takes precedence over String
1598 Unit unit;
1599 Vector<Double> restFreqs;
1600 if (restFreqsIn.nelements()>0) {
1601 restFreqs.resize(restFreqsIn.nelements());
1602 restFreqs = restFreqsIn;
1603 unit = Unit(sUnit);
1604 } else if (lines.size()>0) {
1605 const uInt nLines = lines.size();
1606 unit = Unit("Hz");
1607 restFreqs.resize(nLines);
1608 MFrequency lineFreq;
1609 for (uInt i=0; i<nLines; i++) {
1610 String tS(lines[i]);
1611 tS.upcase();
1612 if (MeasTable::Line(lineFreq, tS)) {
1613 restFreqs[i] = lineFreq.getValue().getValue(); // Hz
1614 } else {
1615 String s = String(lines[i]) +
1616 String(" is an unrecognized spectral line");
1617 throw(AipsError(s));
1618 }
1619 }
1620 } else {
1621 throw(AipsError("You have not specified any rest frequencies or lines"));
1622 }
1623
1624 // If multiple restfreqs, must be length nIF. In this
1625 // case we will just replace the rest frequencies
1626 const uInt nRestFreqs = restFreqs.nelements();
1627 Int idx = -1;
1628 SDFrequencyTable sdft = getSDFreqTable();
1629
1630 if (nRestFreqs>1) {
1631 // Replace restFreqs, one per IF
1632 if (nRestFreqs != nIFs) {
1633 throw (AipsError("Number of rest frequencies must be equal to the number of IFs"));
1634 }
1635 cout << "Replacing rest frequencies with given list, one per IF" << endl;
1636 sdft.deleteRestFrequencies();
1637 for (uInt i=0; i<nRestFreqs; i++) {
1638 Quantum<Double> rf(restFreqs[i], unit);
1639 sdft.addRestFrequency(rf.getValue("Hz"));
1640 }
1641 } else {
1642
1643 // Add new rest freq
1644 Quantum<Double> rf(restFreqs[0], unit);
1645 idx = sdft.addRestFrequency(rf.getValue("Hz"));
1646 cout << "Selecting given rest frequency" << endl;
1647 }
1648
1649 // Replace
1650 Bool empty = source.empty();
1651 Bool ok = False;
1652 if (putSDFreqTable(sdft)) {
1653 const uInt nRow = table_.nrow();
1654 String srcName;
1655 Vector<uInt> restFreqIDs;
1656 for (uInt i=0; i<nRow; i++) {
1657 srcnCol_.get(i, srcName);
1658 restfreqidCol_.get(i,restFreqIDs);
1659
1660 if (idx==-1) {
1661 // Replace vector of restFreqs; one per IF.
1662 // No selection possible
1663 for (uInt i=0; i<nIFs; i++) restFreqIDs[i] = i;
1664 } else {
1665 // Set RestFreqID for selected data
1666 if (empty || source==srcName) {
1667 if (whichIF<0) {
1668 restFreqIDs = idx;
1669 } else {
1670 restFreqIDs[whichIF] = idx;
1671 }
1672 }
1673 }
1674 restfreqidCol_.put(i,restFreqIDs);
1675 }
1676 ok = True;
1677 } else {
1678 ok = False;
1679 }
1680
1681 return ok;
1682}
1683
1684void SDMemTable::spectralLines() const
1685{
1686 Vector<String> lines = MeasTable::Lines();
1687 MFrequency lineFreq;
1688 Double freq;
1689
1690 cout.flags(std::ios_base::left);
1691 cout << "Line Frequency (Hz)" << endl;
1692 cout << "-----------------------" << endl;
1693 for (uInt i=0; i<lines.nelements(); i++) {
1694 MeasTable::Line(lineFreq, lines[i]);
1695 freq = lineFreq.getValue().getValue(); // Hz
1696 cout << setw(11) << lines[i] << setprecision(10) << freq << endl;
1697 }
1698}
1699
1700void SDMemTable::renumber()
1701{
1702 uInt nRow = scanCol_.nrow();
1703 Int newscanid = 0;
1704 Int cIdx;// the current scanid
1705 // get the first scanid
1706 scanCol_.getScalar(0,cIdx);
1707 Int pIdx = cIdx;// the scanid of the previous row
1708 for (uInt i=0; i<nRow;++i) {
1709 scanCol_.getScalar(i,cIdx);
1710 if (pIdx == cIdx) {
1711 // renumber
1712 scanCol_.put(i,newscanid);
1713 } else {
1714 ++newscanid;
1715 pIdx = cIdx; // store scanid
1716 --i; // don't increment next loop
1717 }
1718 }
1719}
1720
1721
1722void SDMemTable::setCursorSlice(IPosition& start, IPosition& end,
1723 const IPosition& shape) const
1724{
1725 const uInt nDim = shape.nelements();
1726 start.resize(nDim);
1727 end.resize(nDim);
1728
1729 start(asap::BeamAxis) = beamSel_;
1730 end(asap::BeamAxis) = beamSel_;
1731 start(asap::IFAxis) = IFSel_;
1732 end(asap::IFAxis) = IFSel_;
1733
1734 start(asap::PolAxis) = polSel_;
1735 end(asap::PolAxis) = polSel_;
1736
1737 start(asap::ChanAxis) = 0;
1738 end(asap::ChanAxis) = shape(asap::ChanAxis) - 1;
1739}
1740
1741
1742std::vector<float> SDMemTable::getFloatSpectrum(const Array<Float>& arr) const
1743 // Get spectrum at cursor location
1744{
1745
1746 // Setup accessors
1747 ArrayAccessor<Float, Axis<asap::BeamAxis> > aa0(arr);
1748 aa0.reset(aa0.begin(uInt(beamSel_))); // Beam selection
1749
1750 ArrayAccessor<Float, Axis<asap::IFAxis> > aa1(aa0);
1751 aa1.reset(aa1.begin(uInt(IFSel_))); // IF selection
1752
1753 ArrayAccessor<Float, Axis<asap::PolAxis> > aa2(aa1);
1754 aa2.reset(aa2.begin(uInt(polSel_))); // Pol selection
1755
1756 std::vector<float> spectrum;
1757 for (ArrayAccessor<Float, Axis<asap::ChanAxis> > i(aa2); i != i.end(); ++i) {
1758 spectrum.push_back(*i);
1759 }
1760 return spectrum;
1761}
1762
Note: See TracBrowser for help on using the repository browser.