source: trunk/src/SDContainer.cc@ 474

Last change on this file since 474 was 473, checked in by kil064, 20 years ago

fix bug in setFlags to do with Xpol and iterating
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 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::appendHistory(const String& hist)
482{
483 history_.resize(history_.nelements()+1,True);
484 history_[history_.nelements()-1] = hist;
485 return True;
486}
487Bool SDContainer::putHistory(const Vector<String>& hist)
488{
489 history_.resize();
490 history_ = hist;
491 return True;
492}
493
494Bool SDContainer::putFitMap(const Array<Int>& arr) {
495 fitIDMap_.resize();
496 fitIDMap_ = arr;
497 return True;
498}
499
500void SDContainer::setSlice(IPosition& start, IPosition& end,
501 const IPosition& shpIn, const IPosition& shpOut,
502 uInt whichBeam, uInt whichIF, Bool tSys,
503 Bool xPol) const
504//
505// tSYs
506// shpIn [nPol]
507// else
508// shpIn [nCHan,nPol]
509//
510// if xPol present, the output nPol = 4 but
511// the input spectra are nPol=2 (tSys) or nPol=2 (spectra)
512//
513{
514 AlwaysAssert(asap::nAxes==4,AipsError);
515 uInt pOff = 0;
516 if (xPol) pOff += 2;
517 if (tSys) {
518 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(0)+pOff,AipsError); // pol
519 } else {
520 AlwaysAssert(shpOut(asap::ChanAxis)==shpIn(0),AipsError); // chan
521 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(1)+pOff,AipsError); // pol
522 }
523//
524 setSlice(start, end, shpOut, whichBeam, whichIF);
525}
526
527
528void SDContainer::setSlice(IPosition& start, IPosition& end,
529 const IPosition& shape,
530 uInt whichBeam, uInt whichIF) const
531{
532 AlwaysAssert(asap::nAxes==4,AipsError);
533 //
534 start.resize(asap::nAxes);
535 start = 0;
536 start(asap::BeamAxis) = whichBeam;
537 start(asap::IFAxis) = whichIF;
538//
539 end.resize(asap::nAxes);
540 end = shape-1;
541 end(asap::BeamAxis) = whichBeam;
542 end(asap::IFAxis) = whichIF;
543}
544
545
546uInt SDFrequencyTable::addFrequency(Double refPix, Double refVal, Double inc)
547{
548 if (length() > 0) {
549 for (uInt i=0; i< length();i++) {
550 if (near(refVal,refVal_[i]) &&
551 near(refPix,refPix_[i]) &&
552 near(inc,increment_[i])) {
553 return i;
554 }
555 }
556 }
557
558// Not found - add it
559
560 nFreq_ += 1;
561 refPix_.resize(nFreq_,True);
562 refVal_.resize(nFreq_,True);
563 increment_.resize(nFreq_,True);
564 refPix_[nFreq_-1] = refPix;
565 refVal_[nFreq_-1] = refVal;
566 increment_[nFreq_-1] = inc;
567 return nFreq_-1;
568}
569
570uInt SDFrequencyTable::addRestFrequency(Double val)
571{
572 uInt nFreq = restFreqs_.nelements();
573 if (nFreq>0) {
574 for (uInt i=0; i<nFreq;i++) {
575 if (near(restFreqs_[i],val)) {
576 return i;
577 }
578 }
579 }
580
581// Not found - add it
582
583 nFreq += 1;
584 restFreqs_.resize(nFreq,True);
585 restFreqs_[nFreq-1] = val;
586 return nFreq-1;
587}
588
589
590void SDFrequencyTable::restFrequencies(Vector<Double>& rfs,
591 String& rfunit ) const
592{
593 rfs.resize(restFreqs_.nelements());
594 rfs = restFreqs_;
595 rfunit = restFreqUnit_;
596}
597
598// SDDataDesc
599
600uInt SDDataDesc::addEntry(const String& source, uInt ID,
601 const MDirection& dir, uInt secID)
602{
603
604// See if already exists
605
606 if (n_ > 0) {
607 for (uInt i=0; i<n_; i++) {
608 if (source==source_[i] && ID==ID_[i]) {
609 return i;
610 }
611 }
612 }
613
614// Not found - add it
615
616 n_ += 1;
617 source_.resize(n_,True);
618 ID_.resize(n_,True);
619 secID_.resize(n_,True);
620 secDir_.resize(n_,True,True);
621//
622 source_[n_-1] = source;
623 ID_[n_-1] = ID;
624 secID_[n_-1] = secID;
625 secDir_[n_-1] = dir;
626//
627 return n_-1;
628}
629
630void SDDataDesc::summary() const
631{
632 if (n_>0) {
633 cerr << "Source ID" << endl;
634 for (uInt i=0; i<n_; i++) {
635 cout << setw(11) << source_(i) << ID_(i) << endl;
636 }
637 }
638}
639
Note: See TracBrowser for help on using the repository browser.