source: branches/Release12/src/SDContainer.cc@ 2582

Last change on this file since 2582 was 685, checked in by mar637, 19 years ago

forgot to change some hardcoded axes. Now using asap::AxisNo everywhere.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.0 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(asap::BeamAxis)),
90 nIF_(shp(asap::IFAxis)),
91 nPol_(shp(asap::PolAxis)),
92 nChan_(shp(asap::ChanAxis)),
93 spectrum_(shp),
94 flags_(shp),
95 tsys_(shp),
96 freqidx_(shp(asap::IFAxis)),
97 restfreqidx_(shp(asap::IFAxis)) {
98 IPosition ip(2,shp(asap::BeamAxis),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(asap::BeamAxis);
110 nIF_ = shp(asap::IFAxis);
111 nPol_ = shp(asap::PolAxis);
112 nChan_ = shp(asap::ChanAxis);
113 spectrum_.resize(shp);
114 flags_.resize(shp);
115 tsys_.resize(shp);
116 freqidx_.resize(shp(asap::IFAxis));
117 restfreqidx_.resize(shp(asap::IFAxis));
118 IPosition ip(2,shp(asap::BeamAxis),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{
192
193// Get slice and check dim
194
195 IPosition start, end;
196 setSlice(start, end, spec.shape(), spectrum_.shape(),
197 whichBeam, whichIF, False, False);
198
199// Get a reference to the Pol/Chan slice we are interested in
200
201 Array<Float> subArr = spectrum_(start,end);
202
203// Iterate through it and fill
204
205 ReadOnlyVectorIterator<Float> inIt(spec,0);
206 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
207 while (!outIt.pastEnd()) {
208 outIt.vector() = inIt.vector();
209//
210 inIt.next();
211 outIt.next();
212 }
213
214 // unset flags for this spectrum, they might be set again by the
215 // setFlags method
216
217 Array<uChar> arr(flags_(start,end));
218 arr = uChar(0);
219//
220 return True;
221}
222
223Bool SDContainer::setFlags(const Matrix<uChar>& flags,
224 uInt whichBeam, uInt whichIF,
225 Bool hasXPol)
226//
227// flags is [nChan,nPol]
228// flags_ is [nBeam,nIF,nPol,nChan]
229//
230// Note that there are no separate flags for the XY cross polarization so we make
231// them up from the X and Y which is all the silly code below. This means
232// that nPol==2 on input but 4 on output
233//
234{
235 if (hasXPol) AlwaysAssert(nPol_==4,AipsError);
236
237// Get slice and check dim
238
239 IPosition start, end;
240 setSlice(start, end, flags.shape(), flags_.shape(),
241 whichBeam, whichIF, False, hasXPol);
242
243// Get a reference to the Pol/Chan slice we are interested in
244
245 Array<uChar> subArr = flags_(start,end);
246
247// Iterate through pol/chan plane and fill
248
249 ReadOnlyVectorIterator<uChar> inIt(flags,0);
250 VectorIterator<uChar> outIt(subArr,asap::ChanAxis);
251//
252 Vector<uChar> maskPol0;
253 Vector<uChar> maskPol1;
254 Vector<uChar> maskPol01;
255 while (!outIt.pastEnd()) {
256 const IPosition& pos = outIt.pos();
257 if (pos(asap::PolAxis)<2) {
258 outIt.vector() = inIt.vector();
259//
260 if (hasXPol) {
261 if (pos(asap::PolAxis)==0) {
262 maskPol0 = inIt.vector();
263 } else if (pos(asap::PolAxis)==1) {
264 maskPol1 = inIt.vector();
265//
266 maskPol01.resize(maskPol0.nelements());
267 for (uInt i=0; i<maskPol01.nelements(); i++) maskPol01[i] = maskPol0[i]&maskPol1[i];
268 }
269 }
270 inIt.next();
271 } else if (pos(asap::PolAxis)==2) {
272 if (hasXPol) {
273 outIt.vector() = maskPol01;
274 } else {
275 outIt.vector() = inIt.vector();
276 inIt.next();
277 }
278
279 } else if (pos(asap::PolAxis)==3) {
280 if (hasXPol) {
281 outIt.vector() = maskPol01;
282 } else {
283 outIt.vector() = inIt.vector();
284 inIt.next();
285 }
286 }
287 outIt.next();
288 }
289//
290 return True;
291}
292
293
294Bool SDContainer::setTsys(const Vector<Float>& tsys,
295 uInt whichBeam, uInt whichIF,
296 Bool hasXpol)
297//
298// TSys shape is [nPol]
299// Tsys does not depend upon channel but is replicated
300// for simplicity of use.
301// There is no Tsys measurement for the cross polarization
302// so I have set TSys for XY to sqrt(Tx*Ty)
303//
304{
305
306// Get slice and check dim
307
308 IPosition start, end;
309 setSlice(start, end, tsys.shape(), tsys_.shape(),
310 whichBeam, whichIF, True, hasXpol);
311
312// Get a reference to the Pol/Chan slice we are interested in
313
314 Array<Float> subArr = tsys_(start,end);
315
316// Iterate through it and fill
317
318 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
319 uInt i=0;
320 while (!outIt.pastEnd()) {
321 const IPosition& pos = outIt.pos();
322//
323 if (pos(asap::PolAxis)<2) {
324 outIt.vector() = tsys(i++);
325 } else {
326 if (hasXpol) {
327 outIt.vector() = sqrt(tsys[0]*tsys[1]);
328 } else {
329 outIt.vector() = tsys(i++);
330 }
331 }
332//
333 outIt.next();
334 }
335}
336
337Array<Float> SDContainer::getSpectrum(uInt whichBeam, uInt whichIF)
338//
339// non-const function because of Array(start,end) slicer
340//
341// Input [nBeam,nIF,nPol,nChan]
342// Output [nChan,nPol]
343//
344{
345
346// Get reference to slice and check dim
347
348 IPosition start, end;
349 setSlice(start, end, spectrum_.shape(), whichBeam, whichIF);
350//
351 Array<Float> dataIn = spectrum_(start,end);
352 Array<Float> dataOut(IPosition(2, nChan_, nPol_));
353//
354 ReadOnlyVectorIterator<Float> itIn(dataIn, asap::ChanAxis);
355 VectorIterator<Float> itOut(dataOut, 0);
356 while (!itOut.pastEnd()) {
357 itOut.vector() = itIn.vector();
358//
359 itIn.next();
360 itOut.next();
361 }
362//
363 return dataOut.copy();
364}
365
366Array<uChar> SDContainer::getFlags(uInt whichBeam, uInt whichIF)
367//
368// non-const function because of Array(start,end) slicer
369//
370// Input [nBeam,nIF,nPol,nChan]
371// Output [nChan,nPol]
372//
373{
374
375// Get reference to slice and check dim
376
377 IPosition start, end;
378 setSlice(start, end, flags_.shape(), whichBeam, whichIF);
379//
380 Array<uChar> dataIn = flags_(start,end);
381 Array<uChar> dataOut(IPosition(2, nChan_, nPol_));
382//
383 ReadOnlyVectorIterator<uChar> itIn(dataIn, asap::ChanAxis);
384 VectorIterator<uChar> itOut(dataOut, 0);
385 while (!itOut.pastEnd()) {
386 itOut.vector() = itIn.vector();
387//
388 itIn.next();
389 itOut.next();
390 }
391//
392 return dataOut.copy();
393}
394
395Array<Float> SDContainer::getTsys(uInt whichBeam, uInt whichIF)
396//
397// Input [nBeam,nIF,nPol,nChan]
398// Output [nPol] (drop channel dependency and select first value only)
399//
400{
401// Get reference to slice and check dim
402
403 IPosition start, end;
404 setSlice(start, end, spectrum_.shape(), whichBeam, whichIF);
405//
406 Array<Float> dataIn = tsys_(start,end);
407 Vector<Float> dataOut(nPol_);
408//
409 ReadOnlyVectorIterator<Float> itIn(dataIn, asap::ChanAxis);
410 VectorIterator<Float> itOut(dataOut, 0);
411 uInt i = 0;
412 while (!itIn.pastEnd()) {
413 dataOut[i++] = itIn.vector()[0];
414 itIn.next();
415 }
416//
417 return dataOut.copy();
418}
419
420
421
422Array<Double> SDContainer::getDirection(uInt whichBeam) const
423//
424// Input [nBeam,2]
425// Output [nBeam]
426//
427{
428 Vector<Double> dataOut(2);
429 dataOut(0) = direction_(IPosition(2,whichBeam,0));
430 dataOut(1) = direction_(IPosition(2,whichBeam,1));
431 return dataOut.copy();
432}
433
434
435Bool SDContainer::setFrequencyMap(uInt freqID, uInt whichIF) {
436 freqidx_[whichIF] = freqID;
437 return True;
438}
439
440Bool SDContainer::putFreqMap(const Vector<uInt>& freqs) {
441 freqidx_.resize();
442 freqidx_ = freqs;
443 return True;
444}
445
446Bool SDContainer::setRestFrequencyMap(uInt freqID, uInt whichIF) {
447 restfreqidx_[whichIF] = freqID;
448 return True;
449}
450
451Bool SDContainer::putRestFreqMap(const Vector<uInt>& freqs) {
452 restfreqidx_.resize();
453 restfreqidx_ = freqs;
454 return True;
455}
456
457Bool SDContainer::setDirection(const Vector<Double>& point, uInt whichBeam)
458//
459// Input [2]
460// Output [nBeam,2]
461//
462{
463 if (point.nelements() != 2) return False;
464//
465 Vector<Double> dataOut(2);
466 direction_(IPosition(2,whichBeam,0)) = point[0];
467 direction_(IPosition(2,whichBeam,1)) = point[1];
468 return True;
469}
470
471
472Bool SDContainer::putDirection(const Array<Double>& dir) {
473 direction_.resize();
474 direction_ = dir;
475 return True;
476}
477
478Bool SDContainer::putFitMap(const Array<Int>& arr) {
479 fitIDMap_.resize();
480 fitIDMap_ = arr;
481 return True;
482}
483
484void SDContainer::setSlice(IPosition& start, IPosition& end,
485 const IPosition& shpIn, const IPosition& shpOut,
486 uInt whichBeam, uInt whichIF, Bool tSys,
487 Bool xPol) const
488//
489// tSYs
490// shpIn [nPol]
491// else
492// shpIn [nCHan,nPol]
493//
494// if xPol present, the output nPol = 4 but
495// the input spectra are nPol=2 (tSys) or nPol=2 (spectra)
496//
497{
498 AlwaysAssert(asap::nAxes==4,AipsError);
499 uInt pOff = 0;
500 if (xPol) pOff += 2;
501 if (tSys) {
502 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(0)+pOff,AipsError); // pol
503 } else {
504 AlwaysAssert(shpOut(asap::ChanAxis)==shpIn(0),AipsError); // chan
505 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(1)+pOff,AipsError); // pol
506 }
507//
508 setSlice(start, end, shpOut, whichBeam, whichIF);
509}
510
511
512void SDContainer::setSlice(IPosition& start, IPosition& end,
513 const IPosition& shape,
514 uInt whichBeam, uInt whichIF) const
515{
516 AlwaysAssert(asap::nAxes==4,AipsError);
517 //
518 start.resize(asap::nAxes);
519 start = 0;
520 start(asap::BeamAxis) = whichBeam;
521 start(asap::IFAxis) = whichIF;
522//
523 end.resize(asap::nAxes);
524 end = shape-1;
525 end(asap::BeamAxis) = whichBeam;
526 end(asap::IFAxis) = whichIF;
527}
528
529
530uInt SDFrequencyTable::addFrequency(Double refPix, Double refVal, Double inc)
531{
532 if (length() > 0) {
533 for (uInt i=0; i< length();i++) {
534 if (near(refVal,refVal_[i]) &&
535 near(refPix,refPix_[i]) &&
536 near(inc,increment_[i])) {
537 return i;
538 }
539 }
540 }
541
542// Not found - add it
543
544 nFreq_ += 1;
545 refPix_.resize(nFreq_,True);
546 refVal_.resize(nFreq_,True);
547 increment_.resize(nFreq_,True);
548 refPix_[nFreq_-1] = refPix;
549 refVal_[nFreq_-1] = refVal;
550 increment_[nFreq_-1] = inc;
551 return nFreq_-1;
552}
553
554uInt SDFrequencyTable::addRestFrequency(Double val)
555{
556 uInt nFreq = restFreqs_.nelements();
557 if (nFreq>0) {
558 for (uInt i=0; i<nFreq;i++) {
559 if (near(restFreqs_[i],val)) {
560 return i;
561 }
562 }
563 }
564
565// Not found - add it
566
567 nFreq += 1;
568 restFreqs_.resize(nFreq,True);
569 restFreqs_[nFreq-1] = val;
570 return nFreq-1;
571}
572
573
574void SDFrequencyTable::restFrequencies(Vector<Double>& rfs,
575 String& rfunit ) const
576{
577 rfs.resize(restFreqs_.nelements());
578 rfs = restFreqs_;
579 rfunit = restFreqUnit_;
580}
581
582// SDDataDesc
583
584uInt SDDataDesc::addEntry(const String& source, uInt ID,
585 const MDirection& dir, uInt secID)
586{
587
588// See if already exists
589
590 if (n_ > 0) {
591 for (uInt i=0; i<n_; i++) {
592 if (source==source_[i] && ID==ID_[i]) {
593 return i;
594 }
595 }
596 }
597
598// Not found - add it
599
600 n_ += 1;
601 source_.resize(n_,True);
602 ID_.resize(n_,True);
603 secID_.resize(n_,True);
604 secDir_.resize(n_,True,True);
605//
606 source_[n_-1] = source;
607 ID_[n_-1] = ID;
608 secID_[n_-1] = secID;
609 secDir_[n_-1] = dir;
610//
611 return n_-1;
612}
613
614void SDDataDesc::summary() const
615{
616 if (n_>0) {
617 cerr << "Source ID" << endl;
618 for (uInt i=0; i<n_; i++) {
619 cout << setw(11) << source_(i) << ID_(i) << endl;
620 }
621 }
622}
623
Note: See TracBrowser for help on using the repository browser.