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
RevLine 
[2]1//#---------------------------------------------------------------------------
2//# SDContainer.cc: A container class for single dish integrations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
[125]5//# ATNF
[2]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//#---------------------------------------------------------------------------
[349]31
[125]32#include <casa/aips.h>
[327]33#include <casa/iostream.h>
34#include <casa/iomanip.h>
[104]35#include <casa/Exceptions.h>
[349]36#include <casa/Utilities/Assert.h>
[81]37#include <tables/Tables/Table.h>
38#include <casa/Arrays/IPosition.h>
[125]39#include <casa/Arrays/Matrix.h>
[349]40#include <casa/Arrays/VectorIter.h>
[81]41#include <casa/Arrays/ArrayAccessor.h>
[417]42#include <casa/Arrays/ArrayMath.h>
[308]43#include <casa/BasicMath/Math.h>
[81]44#include <casa/Quanta/MVTime.h>
[2]45
[308]46
[349]47
[215]48#include "SDDefs.h"
[2]49#include "SDContainer.h"
50
[125]51using namespace casa;
[83]52using namespace asap;
[2]53
[18]54void SDHeader::print() const {
55 MVTime mvt(this->utc);
[47]56 mvt.setFormat(MVTime::YMD);
[18]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): "
[47]67 << mvt
[18]68 << endl;
69 //setprecision(10) << this->utc << endl;
70}
71
72
[16]73SDContainer::SDContainer(uInt nBeam, uInt nIF, uInt nPol, uInt nChan)
[2]74 : nBeam_(nBeam),
[16]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)),
[34]80 tsys_(IPosition(4,nBeam,nIF,nPol,nChan)),
[79]81 freqidx_(nIF),
[388]82 restfreqidx_(nIF),
[79]83 direction_(IPosition(2,nBeam,2)) {
[2]84 uChar x = 0;
85 flags_ = ~x;
[104]86 tcal.resize(2);
[2]87}
88
[16]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),
[34]96 tsys_(shp),
[388]97 freqidx_(shp(1)),
98 restfreqidx_(shp(1)) {
[79]99 IPosition ip(2,shp(0),2);
100 direction_.resize(ip);
[16]101 uChar x = 0;
102 flags_ = ~x;
[104]103 tcal.resize(2);
[16]104}
105
[2]106SDContainer::~SDContainer() {
107}
108
[47]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));
[388]118 restfreqidx_.resize(shp(1));
[79]119 IPosition ip(2,shp(0),2);
120 direction_.resize(ip);
[47]121}
122
[2]123Bool SDContainer::putSpectrum(const Array<Float>& spec) {
124 spectrum_ = spec;
125}
[7]126Bool SDContainer::putFlags(const Array<uChar>& flag) {
127 flags_ = flag;
[2]128}
[7]129Bool SDContainer::putTsys(const Array<Float>& tsys) {
130 tsys_ = tsys;
[2]131}
132
133Bool SDContainer::setSpectrum(const Matrix<Float>& spec,
[417]134 const Vector<Complex>& cSpec,
[349]135 uInt whichBeam, uInt whichIF)
136//
137// spec is [nChan,nPol]
[417]138// cspec is [nChan]
[349]139// spectrum_ is [,,,nChan]
140//
[417]141// nPOl_ = 4 - xx, yy, real(xy), imag(xy)
142//
[349]143{
[417]144 AlwaysAssert(nPol_==4,AipsError);
[2]145
[417]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
[349]195// Get slice and check dim
196
197 IPosition start, end;
198 setSlice (start, end, spec.shape(), spectrum_.shape(),
[417]199 whichBeam, whichIF, False, False);
[349]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);
[417]209 while (!outIt.pastEnd()) {
[349]210 outIt.vector() = inIt.vector();
211//
212 inIt.next();
213 outIt.next();
[14]214 }
[349]215
[2]216 // unset flags for this spectrum, they might be set again by the
217 // setFlags method
[14]218
[2]219 Array<uChar> arr(flags_(start,end));
220 arr = uChar(0);
[336]221//
222 return True;
[2]223}
224
225
[417]226
227
[349]228Bool SDContainer::setFlags (const Matrix<uChar>& flags,
[417]229 uInt whichBeam, uInt whichIF,
230 Bool hasXPol)
[349]231//
232// flags is [nChan,nPol]
233// flags_ is [,,,nChan]
234// How annoying.
[417]235// there are no separate flags for XY so make
236// them up from X and Y
[349]237//
238{
[417]239 if (hasXPol) AlwaysAssert(nPol_==4,AipsError);
240
[349]241// Get slice and check dim
[14]242
[349]243 IPosition start, end;
244 setSlice (start, end, flags.shape(), flags_.shape(),
[417]245 whichBeam, whichIF, False, hasXPol);
[349]246
247// Get a reference to the Pol/Chan slice we are interested in
248
249 Array<uChar> subArr = flags_(start,end);
250
[417]251// Iterate through pol/chan plane and fill
[349]252
253 ReadOnlyVectorIterator<uChar> inIt(flags,0);
254 VectorIterator<uChar> outIt(subArr,asap::ChanAxis);
255//
[417]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//
[349]282 outIt.next();
[2]283 }
[349]284//
[16]285 return True;
[2]286}
287
[349]288
[2]289Bool SDContainer::setTsys(const Vector<Float>& tsys,
[417]290 uInt whichBeam, uInt whichIF,
291 Bool hasXpol)
[349]292//
293// Tsys does not depend upon channel but is replicated
[417]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)
[349]297//
298{
299
300// Get slice and check dim
301
302 IPosition start, end;
303 setSlice (start, end, tsys.shape(), tsys_.shape(),
[417]304 whichBeam, whichIF, True, hasXpol);
[349]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()) {
[417]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//
[349]323 outIt.next();
[7]324 }
[2]325}
[27]326
[125]327Array<Float> SDContainer::getSpectrum(uInt whichBeam, uInt whichIF) const
328{
[27]329 Matrix<Float> spectra(nChan_, nPol_);
330
331 // Beam.
[204]332 ArrayAccessor<Float, Axis<asap::BeamAxis> > i0(spectrum_);
[27]333 i0.reset(i0.begin(whichBeam));
334
335 // IF.
[204]336 ArrayAccessor<Float, Axis<asap::IFAxis> > i1(i0);
[27]337 i1.reset(i1.begin(whichIF));
338
339 // Polarization.
[204]340 ArrayAccessor<Float, Axis<asap::PolAxis> > i2(i1);
341 ArrayAccessor<Float, Axis<asap::IFAxis> > o1(spectra);
[27]342
343 while (i2 != i2.end()) {
344 // Channel.
[204]345 ArrayAccessor<Float, Axis<asap::ChanAxis> > i3(i2);
346 ArrayAccessor<Float, Axis<asap::BeamAxis> > o0(o1);
[27]347
348 while (i3 != i3.end()) {
349 *o0 = *i3;
350
351 i3++;
352 o0++;
353 }
354
355 i2++;
356 o1++;
357 }
358
[67]359 return spectra.copy();
[27]360}
361
[125]362
[67]363Array<uChar> SDContainer::getFlags(uInt whichBeam, uInt whichIF) const
[27]364{
365 Matrix<uChar> flagtra(nChan_, nPol_);
366
367 // Beam.
[204]368 ArrayAccessor<uChar, Axis<asap::BeamAxis> > i0(flags_);
[27]369 i0.reset(i0.begin(whichBeam));
370
371 // IF.
[204]372 ArrayAccessor<uChar, Axis<asap::IFAxis> > i1(i0);
[27]373 i1.reset(i1.begin(whichIF));
374
375 // Polarization.
[204]376 ArrayAccessor<uChar, Axis<asap::PolAxis> > i2(i1);
377 ArrayAccessor<uChar, Axis<asap::IFAxis> > o1(flagtra);
[27]378
379 while (i2 != i2.end()) {
380 // Channel.
[204]381 ArrayAccessor<uChar, Axis<asap::ChanAxis> > i3(i2);
382 ArrayAccessor<uChar, Axis<asap::BeamAxis> > o0(o1);
[27]383
384 while (i3 != i3.end()) {
385 *o0 = *i3;
386
387 i3++;
388 o0++;
389 }
390
391 i2++;
392 o1++;
393 }
394
[67]395 return flagtra.copy();
[27]396}
397
[67]398Array<Float> SDContainer::getTsys(uInt whichBeam, uInt whichIF) const
[27]399{
400 Vector<Float> tsys(nPol_);
401
402 // Beam.
[204]403 ArrayAccessor<Float, Axis<asap::BeamAxis> > i0(tsys_);
[27]404 i0.reset(i0.begin(whichBeam));
405
406 // IF.
[204]407 ArrayAccessor<Float, Axis<asap::IFAxis> > i1(i0);
[27]408 i1.reset(i1.begin(whichIF));
409
410 // Channel.
[204]411 ArrayAccessor<Float, Axis<asap::ChanAxis> > i3(i1);
[27]412
413 // Polarization.
[204]414 ArrayAccessor<Float, Axis<asap::PolAxis> > i2(i3);
415 ArrayAccessor<Float, Axis<asap::BeamAxis> > o0(tsys);
[27]416
417 while (i2 != i2.end()) {
418 *o0 = *i2;
419
420 i2++;
421 o0++;
422 }
[67]423 return tsys.copy();
[27]424}
[34]425
[79]426Array<Double> SDContainer::getDirection(uInt whichBeam) const {
427 Vector<Double> direct(2);
[204]428 ArrayAccessor<Double, Axis<asap::BeamAxis> > i0(direction_);
[79]429 i0.reset(i0.begin(whichBeam));
[204]430 ArrayAccessor<Double, Axis<asap::BeamAxis> > o0(direct);
431 ArrayAccessor<Double, Axis<asap::IFAxis> > i1(i0);
[79]432 while (i1 != i1.end()) {
433 *o0 = *i1;
434 i1++;
435 o0++;
436 }
437 return direct.copy();
438}
439
440
[388]441Bool SDContainer::setFrequencyMap(uInt freqID, uInt whichIF) {
442 freqidx_[whichIF] = freqID;
[34]443 return True;
444}
445
446Bool SDContainer::putFreqMap(const Vector<uInt>& freqs) {
447 freqidx_.resize();
448 freqidx_ = freqs;
449 return True;
450}
451
[388]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
[79]463Bool SDContainer::setDirection(const Vector<Double>& point, uInt whichBeam) {
464 if (point.nelements() != 2) return False;
[204]465 ArrayAccessor<Double, Axis<asap::BeamAxis> > aa0(direction_);
[79]466 aa0.reset(aa0.begin(whichBeam));
[204]467 ArrayAccessor<Double, Axis<asap::BeamAxis> > jj(point);
468 for (ArrayAccessor<Double, Axis<asap::IFAxis> > i(aa0);i != i.end(); ++i) {
[79]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
[204]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
[349]495void SDContainer::setSlice (IPosition& start, IPosition& end,
496 const IPosition& shpIn, const IPosition& shpOut,
[417]497 uInt whichBeam, uInt whichIF, Bool tSys,
498 Bool xPol) const
[349]499//
500// tSYs
501// shpIn [nPol]
502// else
503// shpIn [nCHan,nPol]
504//
[417]505// if xPol present, the output nPol = 4 but
506// the input spectra are nPol=2 (tSys) or nPol=2 (spectra)
507//
[349]508{
509 AlwaysAssert(asap::nAxes==4,AipsError);
[417]510 uInt pOff = 0;
511 if (xPol) pOff += 2;
[349]512 if (tSys) {
[417]513 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(0)+pOff,AipsError); // pol
[349]514 } else {
[417]515 AlwaysAssert(shpOut(asap::ChanAxis)==shpIn(0),AipsError); // chan
516 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(1)+pOff,AipsError); // pol
[349]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
[388]531uInt SDFrequencyTable::addFrequency(Double refPix, Double refVal, Double inc)
532{
[34]533 if (length() > 0) {
[388]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;
[34]539 }
540 }
541 }
[388]542
543// Not found - add it
544
[34]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;
[388]552 return nFreq_-1;
[34]553}
554
[388]555uInt SDFrequencyTable::addRestFrequency(Double val)
[204]556{
[388]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;
[204]562 }
563 }
564 }
[388]565
566// Not found - add it
567
568 nFreq += 1;
569 restFreqs_.resize(nFreq,True);
570 restFreqs_[nFreq-1] = val;
571 return nFreq-1;
[204]572}
[388]573
574
[204]575void SDFrequencyTable::restFrequencies(Vector<Double>& rfs,
576 String& rfunit ) const
577{
578 rfs.resize(restFreqs_.nelements());
579 rfs = restFreqs_;
580 rfunit = restFreqUnit_;
581}
[326]582
583
[349]584
585// SDDataDesc
586
[396]587uInt SDDataDesc::addEntry (const String& source, uInt ID,
588 const MDirection& dir, uInt secID)
[326]589{
590
591// See if already exists
592
593 if (n_ > 0) {
594 for (uInt i=0; i<n_; i++) {
[396]595 if (source==source_[i] && ID==ID_[i]) {
[326]596 return i;
597 }
598 }
599 }
600
601// Not found - add it
602
603 n_ += 1;
604 source_.resize(n_,True);
[396]605 ID_.resize(n_,True);
606 secID_.resize(n_,True);
607 secDir_.resize(n_,True,True);
[326]608//
609 source_[n_-1] = source;
[396]610 ID_[n_-1] = ID;
611 secID_[n_-1] = secID;
612 secDir_[n_-1] = dir;
[326]613//
614 return n_-1;
615}
616
617void SDDataDesc::summary() const
618{
[396]619 if (n_>0) {
620 cerr << "Source ID" << endl;
621 for (uInt i=0; i<n_; i++) {
[414]622 cout << setw(11) << source_(i) << ID_(i) << endl;
[396]623 }
[326]624 }
625}
Note: See TracBrowser for help on using the repository browser.