source: trunk/src/SDContainer.cc@ 452

Last change on this file since 452 was 449, checked in by kil064, 20 years ago

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