source: branches/Release-2-fixes/src/SDContainer.cc@ 740

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