source: trunk/src/SDMemTable.cc@ 491

Last change on this file since 491 was 490, checked in by kil064, 20 years ago

add functions

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