source: trunk/src/SDContainer.cc@ 466

Last change on this file since 466 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: 14.6 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDContainer.cc: A container class 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 <casa/aips.h>
33#include <casa/iostream.h>
34#include <casa/iomanip.h>
35#include <casa/Exceptions.h>
36#include <casa/Utilities/Assert.h>
37#include <tables/Tables/Table.h>
38#include <casa/Arrays/IPosition.h>
39#include <casa/Arrays/Matrix.h>
40#include <casa/Arrays/VectorIter.h>
41#include <casa/Arrays/ArrayMath.h>
42#include <casa/BasicMath/Math.h>
43#include <casa/Quanta/MVTime.h>
44
45
46
47#include "SDDefs.h"
48#include "SDContainer.h"
49
50using namespace casa;
51using namespace asap;
52
53void SDHeader::print() const {
54 MVTime mvt(this->utc);
55 mvt.setFormat(MVTime::YMD);
56 cout << "Observer: " << this->observer << endl
57 << "Project: " << this->project << endl
58 << "Obstype: " << this->obstype << endl
59 << "Antenna: " << this->antennaname << endl
60 << "Ant. Position: " << this->antennaposition << endl
61 << "Equinox: " << this->equinox << endl
62 << "Freq. ref.: " << this->freqref << endl
63 << "Ref. frequency: " << this->reffreq << endl
64 << "Bandwidth: " << this->bandwidth << endl
65 << "Time (utc): "
66 << mvt
67 << endl;
68 //setprecision(10) << this->utc << endl;
69}
70
71
72SDContainer::SDContainer(uInt nBeam, uInt nIF, uInt nPol, uInt nChan)
73 : nBeam_(nBeam),
74 nIF_(nIF),
75 nPol_(nPol),
76 nChan_(nChan),
77 spectrum_(IPosition(4,nBeam,nIF,nPol,nChan)),
78 flags_(IPosition(4,nBeam,nIF,nPol,nChan)),
79 tsys_(IPosition(4,nBeam,nIF,nPol,nChan)),
80 freqidx_(nIF),
81 restfreqidx_(nIF),
82 direction_(IPosition(2,nBeam,2)) {
83 uChar x = 0;
84 flags_ = ~x;
85 tcal.resize(2);
86}
87
88SDContainer::SDContainer(IPosition shp)
89 : nBeam_(shp(0)),
90 nIF_(shp(1)),
91 nPol_(shp(2)),
92 nChan_(shp(3)),
93 spectrum_(shp),
94 flags_(shp),
95 tsys_(shp),
96 freqidx_(shp(1)),
97 restfreqidx_(shp(1)) {
98 IPosition ip(2,shp(0),2);
99 direction_.resize(ip);
100 uChar x = 0;
101 flags_ = ~x;
102 tcal.resize(2);
103}
104
105SDContainer::~SDContainer() {
106}
107
108Bool SDContainer::resize(IPosition shp) {
109 nBeam_ = shp(0);
110 nIF_ = shp(1);
111 nPol_ = shp(2);
112 nChan_ = shp(3);
113 spectrum_.resize(shp);
114 flags_.resize(shp);
115 tsys_.resize(shp);
116 freqidx_.resize(shp(1));
117 restfreqidx_.resize(shp(1));
118 IPosition ip(2,shp(0),2);
119 direction_.resize(ip);
120}
121
122Bool SDContainer::putSpectrum(const Array<Float>& spec) {
123 spectrum_ = spec;
124}
125Bool SDContainer::putFlags(const Array<uChar>& flag) {
126 flags_ = flag;
127}
128Bool SDContainer::putTsys(const Array<Float>& tsys) {
129 tsys_ = tsys;
130}
131
132Bool SDContainer::setSpectrum(const Matrix<Float>& spec,
133 const Vector<Complex>& cSpec,
134 uInt whichBeam, uInt whichIF)
135//
136// spec is [nChan,nPol]
137// cspec is [nChan]
138// spectrum_ is [,,,nChan]
139//
140// nPOl_ = 4 - xx, yy, real(xy), imag(xy)
141//
142{
143 AlwaysAssert(nPol_==4,AipsError);
144
145// Get slice and check dim
146
147 Bool tSys = False;
148 Bool xPol = True;
149 IPosition start, end;
150 setSlice (start, end, spec.shape(), spectrum_.shape(),
151 whichBeam, whichIF, tSys, xPol);
152
153// Get a reference to the Pol/Chan slice we are interested in
154
155 Array<Float> subArr = spectrum_(start,end);
156
157// Iterate through pol-chan plane and fill
158
159 ReadOnlyVectorIterator<Float> inIt(spec,0);
160 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
161 while (!outIt.pastEnd()) {
162 const IPosition& pos = outIt.pos();
163 if (pos(asap::PolAxis)<2) {
164 outIt.vector() = inIt.vector();
165 inIt.next();
166 } else if (pos(asap::PolAxis)==2) {
167 outIt.vector() = real(cSpec);
168 } else if (pos(asap::PolAxis)==3) {
169 outIt.vector() = imag(cSpec);
170 }
171//
172 outIt.next();
173 }
174
175 // unset flags for this spectrum, they might be set again by the
176 // setFlags method
177
178 Array<uChar> arr(flags_(start,end));
179 arr = uChar(0);
180//
181 return True;
182}
183
184
185Bool SDContainer::setSpectrum(const Matrix<Float>& spec,
186 uInt whichBeam, uInt whichIF)
187//
188// spec is [nChan,nPol]
189// spectrum_ is [,,,nChan]
190// How annoying.
191// nPol==2
192{
193
194// Get slice and check dim
195
196 IPosition start, end;
197 setSlice (start, end, spec.shape(), spectrum_.shape(),
198 whichBeam, whichIF, False, False);
199
200// Get a reference to the Pol/Chan slice we are interested in
201
202 Array<Float> subArr = spectrum_(start,end);
203
204// Iterate through it and fill
205
206 ReadOnlyVectorIterator<Float> inIt(spec,0);
207 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
208 while (!outIt.pastEnd()) {
209 outIt.vector() = inIt.vector();
210//
211 inIt.next();
212 outIt.next();
213 }
214
215 // unset flags for this spectrum, they might be set again by the
216 // setFlags method
217
218 Array<uChar> arr(flags_(start,end));
219 arr = uChar(0);
220//
221 return True;
222}
223
224
225
226
227Bool SDContainer::setFlags(const Matrix<uChar>& flags,
228 uInt whichBeam, uInt whichIF,
229 Bool hasXPol)
230//
231// flags is [nChan,nPol]
232// flags_ is [,,,nChan]
233// How annoying.
234// there are no separate flags for XY so make
235// them up from X and Y
236//
237{
238 if (hasXPol) AlwaysAssert(nPol_==4,AipsError);
239
240// Get slice and check dim
241
242 IPosition start, end;
243 setSlice (start, end, flags.shape(), flags_.shape(),
244 whichBeam, whichIF, False, hasXPol);
245
246// Get a reference to the Pol/Chan slice we are interested in
247
248 Array<uChar> subArr = flags_(start,end);
249
250// Iterate through pol/chan plane and fill
251
252 ReadOnlyVectorIterator<uChar> inIt(flags,0);
253 VectorIterator<uChar> outIt(subArr,asap::ChanAxis);
254//
255 Vector<uChar> maskPol0;
256 Vector<uChar> maskPol1;
257 Vector<uChar> maskPol01;
258 while (!outIt.pastEnd()) {
259 const IPosition& pos = outIt.pos();
260 if (pos(asap::PolAxis)<2) {
261 outIt.vector() = inIt.vector();
262//
263 if (hasXPol) {
264 if (pos(asap::PolAxis)==0) {
265 maskPol0 = inIt.vector();
266 } else if (pos(asap::PolAxis)==1) {
267 maskPol1 = inIt.vector();
268//
269 maskPol01.resize(maskPol0.nelements());
270 for (uInt i=0; i<maskPol01.nelements(); i++) maskPol01[i] = maskPol0[i]&maskPol1[i];
271 }
272 }
273//
274 inIt.next();
275 } else if (pos(asap::PolAxis)==2) {
276 outIt.vector() = maskPol01;
277 } else if (pos(asap::PolAxis)==3) {
278 outIt.vector() = maskPol01;
279 }
280//
281 outIt.next();
282 }
283//
284 return True;
285}
286
287
288Bool SDContainer::setTsys(const Vector<Float>& tsys,
289 uInt whichBeam, uInt whichIF,
290 Bool hasXpol)
291//
292// Tsys does not depend upon channel but is replicated
293// for simplicity of use.
294// There is no Tsys measurement for the cross polarization
295// so I have set TSys for XY to sqrt(Tx*Ty)
296//
297{
298
299// Get slice and check dim
300
301 IPosition start, end;
302 setSlice (start, end, tsys.shape(), tsys_.shape(),
303 whichBeam, whichIF, True, hasXpol);
304
305// Get a reference to the Pol/Chan slice we are interested in
306
307 Array<Float> subArr = tsys_(start,end);
308
309// Iterate through it and fill
310
311 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
312 uInt i=0;
313 while (!outIt.pastEnd()) {
314 const IPosition& pos = outIt.pos();
315//
316 if (pos(asap::PolAxis)<2) {
317 outIt.vector() = tsys(i++);
318 } else {
319 outIt.vector() = sqrt(tsys[0]*tsys[1]);
320 }
321//
322 outIt.next();
323 }
324}
325
326Array<Float> SDContainer::getSpectrum(uInt whichBeam, uInt whichIF)
327//
328// non-const function because of Array(start,end) slicer
329//
330// Input [nBeam,nIF,nPol,nChan]
331// Output [nChan,nPol]
332//
333{
334
335// Get reference to slice and check dim
336
337 IPosition start, end;
338 setSlice (start, end, spectrum_.shape(), whichBeam, whichIF);
339//
340 Array<Float> dataIn = spectrum_(start,end);
341 Array<Float> dataOut(IPosition(2, nChan_, nPol_));
342//
343 ReadOnlyVectorIterator<Float> itIn(dataIn, asap::ChanAxis);
344 VectorIterator<Float> itOut(dataOut, 0);
345 while (!itOut.pastEnd()) {
346 itOut.vector() = itIn.vector();
347//
348 itIn.next();
349 itOut.next();
350 }
351//
352 return dataOut.copy();
353}
354
355Array<uChar> SDContainer::getFlags(uInt whichBeam, uInt whichIF)
356//
357// non-const function because of Array(start,end) slicer
358//
359// Input [nBeam,nIF,nPol,nChan]
360// Output [nChan,nPol]
361//
362{
363
364// Get reference to slice and check dim
365
366 IPosition start, end;
367 setSlice (start, end, flags_.shape(), whichBeam, whichIF);
368//
369 Array<uChar> dataIn = flags_(start,end);
370 Array<uChar> dataOut(IPosition(2, nChan_, nPol_));
371//
372 ReadOnlyVectorIterator<uChar> itIn(dataIn, asap::ChanAxis);
373 VectorIterator<uChar> itOut(dataOut, 0);
374 while (!itOut.pastEnd()) {
375 itOut.vector() = itIn.vector();
376//
377 itIn.next();
378 itOut.next();
379 }
380//
381 return dataOut.copy();
382}
383
384Array<Float> SDContainer::getTsys(uInt whichBeam, uInt whichIF)
385//
386// Input [nBeam,nIF,nPol,nChan]
387// Output [nPol] (drop channel dependency and select first value)
388//
389{
390// Get reference to slice and check dim
391
392 IPosition start, end;
393 setSlice (start, end, spectrum_.shape(), whichBeam, whichIF);
394//
395 Array<Float> dataIn = tsys_(start,end);
396 Vector<Float> dataOut(nPol_);
397//
398 ReadOnlyVectorIterator<Float> itIn(dataIn, asap::ChanAxis);
399 VectorIterator<Float> itOut(dataOut, 0);
400 uInt i = 0;
401 while (!itIn.pastEnd()) {
402 dataOut[i++] = itIn.vector()[0];
403 itIn.next();
404 }
405//
406 return dataOut.copy();
407}
408
409
410
411Array<Double> SDContainer::getDirection(uInt whichBeam) const
412//
413// Input [nBeam,2]
414// Output [nBeam]
415//
416{
417 Vector<Double> dataOut(2);
418 dataOut(0) = direction_(IPosition(2,whichBeam,0));
419 dataOut(1) = direction_(IPosition(2,whichBeam,1));
420 return dataOut.copy();
421}
422
423
424Bool SDContainer::setFrequencyMap(uInt freqID, uInt whichIF) {
425 freqidx_[whichIF] = freqID;
426 return True;
427}
428
429Bool SDContainer::putFreqMap(const Vector<uInt>& freqs) {
430 freqidx_.resize();
431 freqidx_ = freqs;
432 return True;
433}
434
435Bool SDContainer::setRestFrequencyMap(uInt freqID, uInt whichIF) {
436 restfreqidx_[whichIF] = freqID;
437 return True;
438}
439
440Bool SDContainer::putRestFreqMap(const Vector<uInt>& freqs) {
441 restfreqidx_.resize();
442 restfreqidx_ = freqs;
443 return True;
444}
445
446Bool SDContainer::setDirection(const Vector<Double>& point, uInt whichBeam)
447//
448// Input [2]
449// Output [nBeam,2]
450//
451{
452 if (point.nelements() != 2) return False;
453//
454 Vector<Double> dataOut(2);
455 direction_(IPosition(2,whichBeam,0)) = point[0];
456 direction_(IPosition(2,whichBeam,1)) = point[1];
457 return True;
458}
459
460
461Bool SDContainer::putDirection(const Array<Double>& dir) {
462 direction_.resize();
463 direction_ = dir;
464 return True;
465}
466
467Bool SDContainer::appendHistory(const String& hist)
468{
469 history_.resize(history_.nelements()+1,True);
470 history_[history_.nelements()-1] = hist;
471 return True;
472}
473Bool SDContainer::putHistory(const Vector<String>& hist)
474{
475 history_.resize();
476 history_ = hist;
477 return True;
478}
479
480Bool SDContainer::putFitMap(const Array<Int>& arr) {
481 fitIDMap_.resize();
482 fitIDMap_ = arr;
483 return True;
484}
485
486void SDContainer::setSlice(IPosition& start, IPosition& end,
487 const IPosition& shpIn, const IPosition& shpOut,
488 uInt whichBeam, uInt whichIF, Bool tSys,
489 Bool xPol) const
490//
491// tSYs
492// shpIn [nPol]
493// else
494// shpIn [nCHan,nPol]
495//
496// if xPol present, the output nPol = 4 but
497// the input spectra are nPol=2 (tSys) or nPol=2 (spectra)
498//
499{
500 AlwaysAssert(asap::nAxes==4,AipsError);
501 uInt pOff = 0;
502 if (xPol) pOff += 2;
503 if (tSys) {
504 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(0)+pOff,AipsError); // pol
505 } else {
506 AlwaysAssert(shpOut(asap::ChanAxis)==shpIn(0),AipsError); // chan
507 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(1)+pOff,AipsError); // pol
508 }
509//
510 setSlice(start, end, shpOut, whichBeam, whichIF);
511}
512
513
514void SDContainer::setSlice(IPosition& start, IPosition& end,
515 const IPosition& shape,
516 uInt whichBeam, uInt whichIF) const
517{
518 AlwaysAssert(asap::nAxes==4,AipsError);
519 //
520 start.resize(asap::nAxes);
521 start = 0;
522 start(asap::BeamAxis) = whichBeam;
523 start(asap::IFAxis) = whichIF;
524//
525 end.resize(asap::nAxes);
526 end = shape-1;
527 end(asap::BeamAxis) = whichBeam;
528 end(asap::IFAxis) = whichIF;
529}
530
531
532uInt SDFrequencyTable::addFrequency(Double refPix, Double refVal, Double inc)
533{
534 if (length() > 0) {
535 for (uInt i=0; i< length();i++) {
536 if (near(refVal,refVal_[i]) &&
537 near(refPix,refPix_[i]) &&
538 near(inc,increment_[i])) {
539 return i;
540 }
541 }
542 }
543
544// Not found - add it
545
546 nFreq_ += 1;
547 refPix_.resize(nFreq_,True);
548 refVal_.resize(nFreq_,True);
549 increment_.resize(nFreq_,True);
550 refPix_[nFreq_-1] = refPix;
551 refVal_[nFreq_-1] = refVal;
552 increment_[nFreq_-1] = inc;
553 return nFreq_-1;
554}
555
556uInt SDFrequencyTable::addRestFrequency(Double val)
557{
558 uInt nFreq = restFreqs_.nelements();
559 if (nFreq>0) {
560 for (uInt i=0; i<nFreq;i++) {
561 if (near(restFreqs_[i],val)) {
562 return i;
563 }
564 }
565 }
566
567// Not found - add it
568
569 nFreq += 1;
570 restFreqs_.resize(nFreq,True);
571 restFreqs_[nFreq-1] = val;
572 return nFreq-1;
573}
574
575
576void SDFrequencyTable::restFrequencies(Vector<Double>& rfs,
577 String& rfunit ) const
578{
579 rfs.resize(restFreqs_.nelements());
580 rfs = restFreqs_;
581 rfunit = restFreqUnit_;
582}
583
584// SDDataDesc
585
586uInt SDDataDesc::addEntry(const String& source, uInt ID,
587 const MDirection& dir, uInt secID)
588{
589
590// See if already exists
591
592 if (n_ > 0) {
593 for (uInt i=0; i<n_; i++) {
594 if (source==source_[i] && ID==ID_[i]) {
595 return i;
596 }
597 }
598 }
599
600// Not found - add it
601
602 n_ += 1;
603 source_.resize(n_,True);
604 ID_.resize(n_,True);
605 secID_.resize(n_,True);
606 secDir_.resize(n_,True,True);
607//
608 source_[n_-1] = source;
609 ID_[n_-1] = ID;
610 secID_[n_-1] = secID;
611 secDir_[n_-1] = dir;
612//
613 return n_-1;
614}
615
616void SDDataDesc::summary() const
617{
618 if (n_>0) {
619 cerr << "Source ID" << endl;
620 for (uInt i=0; i<n_; i++) {
621 cout << setw(11) << source_(i) << ID_(i) << endl;
622 }
623 }
624}
625
Note: See TracBrowser for help on using the repository browser.