source: trunk/src/SDContainer.cc@ 425

Last change on this file since 425 was 417, checked in by kil064, 20 years ago

handle the cross correlation

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