source: trunk/src/SDContainer.cc@ 379

Last change on this file since 379 was 349, checked in by kil064, 20 years ago

rewrite setSpectrum, setFLags and setTsys using Array(start,end)
slice and VectorIterator for access, replacing the ArrayAccessor stuff which
is harder to read. I wouldn't have bothered, but to my surprise,
the new implementation is about 10-20% faster. Need to do the get* methods
too.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.9 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/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 direction_(IPosition(2,nBeam,2)) {
82 uChar x = 0;
83 flags_ = ~x;
84 tcal.resize(2);
85}
86
87SDContainer::SDContainer(IPosition shp)
88 : nBeam_(shp(0)),
89 nIF_(shp(1)),
90 nPol_(shp(2)),
91 nChan_(shp(3)),
92 spectrum_(shp),
93 flags_(shp),
94 tsys_(shp),
95 freqidx_(shp(1)) {
96 IPosition ip(2,shp(0),2);
97 direction_.resize(ip);
98 uChar x = 0;
99 flags_ = ~x;
100 tcal.resize(2);
101}
102
103SDContainer::~SDContainer() {
104}
105
106Bool SDContainer::resize(IPosition shp) {
107 nBeam_ = shp(0);
108 nIF_ = shp(1);
109 nPol_ = shp(2);
110 nChan_ = shp(3);
111 spectrum_.resize(shp);
112 flags_.resize(shp);
113 tsys_.resize(shp);
114 freqidx_.resize(shp(1));
115 IPosition ip(2,shp(0),2);
116 direction_.resize(ip);
117}
118
119Bool SDContainer::putSpectrum(const Array<Float>& spec) {
120 spectrum_ = spec;
121}
122Bool SDContainer::putFlags(const Array<uChar>& flag) {
123 flags_ = flag;
124}
125Bool SDContainer::putTsys(const Array<Float>& tsys) {
126 tsys_ = tsys;
127}
128
129Bool SDContainer::setSpectrum(const Matrix<Float>& spec,
130 uInt whichBeam, uInt whichIF)
131//
132// spec is [nChan,nPol]
133// spectrum_ is [,,,nChan]
134// How annoying.
135//
136{
137
138// Get slice and check dim
139
140 IPosition start, end;
141 setSlice (start, end, spec.shape(), spectrum_.shape(),
142 whichBeam, whichIF, False);
143
144// Get a reference to the Pol/Chan slice we are interested in
145
146 Array<Float> subArr = spectrum_(start,end);
147
148// Iterate through it and fill
149
150 ReadOnlyVectorIterator<Float> inIt(spec,0);
151 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
152 while (!inIt.pastEnd()) {
153 outIt.vector() = inIt.vector();
154//
155 inIt.next();
156 outIt.next();
157 }
158
159 // unset flags for this spectrum, they might be set again by the
160 // setFlags method
161
162 Array<uChar> arr(flags_(start,end));
163 arr = uChar(0);
164//
165 return True;
166}
167
168
169Bool SDContainer::setFlags (const Matrix<uChar>& flags,
170 uInt whichBeam, uInt whichIF)
171//
172// flags is [nChan,nPol]
173// flags_ is [,,,nChan]
174// How annoying.
175//
176{
177// Get slice and check dim
178
179 IPosition start, end;
180 setSlice (start, end, flags.shape(), flags_.shape(),
181 whichBeam, whichIF, False);
182
183// Get a reference to the Pol/Chan slice we are interested in
184
185 Array<uChar> subArr = flags_(start,end);
186
187// Iterate through it and fill
188
189 ReadOnlyVectorIterator<uChar> inIt(flags,0);
190 VectorIterator<uChar> outIt(subArr,asap::ChanAxis);
191 while (!inIt.pastEnd()) {
192 outIt.vector() = inIt.vector();
193//
194 inIt.next();
195 outIt.next();
196 }
197//
198 return True;
199}
200
201
202Bool SDContainer::setTsys(const Vector<Float>& tsys,
203 uInt whichBeam, uInt whichIF)
204//
205// Tsys does not depend upon channel but is replicated
206// for simplicity of use
207//
208{
209
210// Get slice and check dim
211
212 IPosition start, end;
213 setSlice (start, end, tsys.shape(), tsys_.shape(),
214 whichBeam, whichIF, True);
215
216// Get a reference to the Pol/Chan slice we are interested in
217
218 Array<Float> subArr = tsys_(start,end);
219
220// Iterate through it and fill
221
222 VectorIterator<Float> outIt(subArr,asap::ChanAxis);
223 uInt i=0;
224 while (!outIt.pastEnd()) {
225 outIt.vector() = tsys(i++);
226 outIt.next();
227 }
228}
229
230Array<Float> SDContainer::getSpectrum(uInt whichBeam, uInt whichIF) const
231{
232 Matrix<Float> spectra(nChan_, nPol_);
233
234 // Beam.
235 ArrayAccessor<Float, Axis<asap::BeamAxis> > i0(spectrum_);
236 i0.reset(i0.begin(whichBeam));
237
238 // IF.
239 ArrayAccessor<Float, Axis<asap::IFAxis> > i1(i0);
240 i1.reset(i1.begin(whichIF));
241
242 // Polarization.
243 ArrayAccessor<Float, Axis<asap::PolAxis> > i2(i1);
244 ArrayAccessor<Float, Axis<asap::IFAxis> > o1(spectra);
245
246 while (i2 != i2.end()) {
247 // Channel.
248 ArrayAccessor<Float, Axis<asap::ChanAxis> > i3(i2);
249 ArrayAccessor<Float, Axis<asap::BeamAxis> > o0(o1);
250
251 while (i3 != i3.end()) {
252 *o0 = *i3;
253
254 i3++;
255 o0++;
256 }
257
258 i2++;
259 o1++;
260 }
261
262 return spectra.copy();
263}
264
265
266Array<uChar> SDContainer::getFlags(uInt whichBeam, uInt whichIF) const
267{
268 Matrix<uChar> flagtra(nChan_, nPol_);
269
270 // Beam.
271 ArrayAccessor<uChar, Axis<asap::BeamAxis> > i0(flags_);
272 i0.reset(i0.begin(whichBeam));
273
274 // IF.
275 ArrayAccessor<uChar, Axis<asap::IFAxis> > i1(i0);
276 i1.reset(i1.begin(whichIF));
277
278 // Polarization.
279 ArrayAccessor<uChar, Axis<asap::PolAxis> > i2(i1);
280 ArrayAccessor<uChar, Axis<asap::IFAxis> > o1(flagtra);
281
282 while (i2 != i2.end()) {
283 // Channel.
284 ArrayAccessor<uChar, Axis<asap::ChanAxis> > i3(i2);
285 ArrayAccessor<uChar, Axis<asap::BeamAxis> > o0(o1);
286
287 while (i3 != i3.end()) {
288 *o0 = *i3;
289
290 i3++;
291 o0++;
292 }
293
294 i2++;
295 o1++;
296 }
297
298 return flagtra.copy();
299}
300
301Array<Float> SDContainer::getTsys(uInt whichBeam, uInt whichIF) const
302{
303 Vector<Float> tsys(nPol_);
304
305 // Beam.
306 ArrayAccessor<Float, Axis<asap::BeamAxis> > i0(tsys_);
307 i0.reset(i0.begin(whichBeam));
308
309 // IF.
310 ArrayAccessor<Float, Axis<asap::IFAxis> > i1(i0);
311 i1.reset(i1.begin(whichIF));
312
313 // Channel.
314 ArrayAccessor<Float, Axis<asap::ChanAxis> > i3(i1);
315
316 // Polarization.
317 ArrayAccessor<Float, Axis<asap::PolAxis> > i2(i3);
318 ArrayAccessor<Float, Axis<asap::BeamAxis> > o0(tsys);
319
320 while (i2 != i2.end()) {
321 *o0 = *i2;
322
323 i2++;
324 o0++;
325 }
326 return tsys.copy();
327}
328
329Array<Double> SDContainer::getDirection(uInt whichBeam) const {
330 Vector<Double> direct(2);
331 ArrayAccessor<Double, Axis<asap::BeamAxis> > i0(direction_);
332 i0.reset(i0.begin(whichBeam));
333 ArrayAccessor<Double, Axis<asap::BeamAxis> > o0(direct);
334 ArrayAccessor<Double, Axis<asap::IFAxis> > i1(i0);
335 while (i1 != i1.end()) {
336 *o0 = *i1;
337 i1++;
338 o0++;
339 }
340 return direct.copy();
341}
342
343
344Bool SDContainer::setFrequencyMap(uInt freqslot, uInt whichIF) {
345 freqidx_[whichIF] = freqslot;
346 return True;
347}
348
349Bool SDContainer::putFreqMap(const Vector<uInt>& freqs) {
350 freqidx_.resize();
351 freqidx_ = freqs;
352 return True;
353}
354
355Bool SDContainer::setDirection(const Vector<Double>& point, uInt whichBeam) {
356 if (point.nelements() != 2) return False;
357 ArrayAccessor<Double, Axis<asap::BeamAxis> > aa0(direction_);
358 aa0.reset(aa0.begin(whichBeam));
359 ArrayAccessor<Double, Axis<asap::BeamAxis> > jj(point);
360 for (ArrayAccessor<Double, Axis<asap::IFAxis> > i(aa0);i != i.end(); ++i) {
361
362 (*i) = (*jj);
363 jj++;
364 }
365 return True;
366}
367
368Bool SDContainer::putDirection(const Array<Double>& dir) {
369 direction_.resize();
370 direction_ = dir;
371 return True;
372}
373
374Bool SDContainer::appendHistory(const String& hist)
375{
376 history_.resize(history_.nelements()+1,True);
377 history_[history_.nelements()-1] = hist;
378 return True;
379}
380Bool SDContainer::putHistory(const Vector<String>& hist)
381{
382 history_.resize();
383 history_ = hist;
384 return True;
385}
386
387void SDContainer::setSlice (IPosition& start, IPosition& end,
388 const IPosition& shpIn, const IPosition& shpOut,
389 uInt whichBeam, uInt whichIF, Bool tSys) const
390//
391// tSYs
392// shpIn [nPol]
393// else
394// shpIn [nCHan,nPol]
395//
396{
397 AlwaysAssert(asap::nAxes==4,AipsError);
398 if (tSys) {
399 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(0),AipsError); // pol
400 } else {
401 AlwaysAssert(shpOut(asap::ChanAxis)==shpIn(0),AipsError); // chan
402 AlwaysAssert(shpOut(asap::PolAxis)==shpIn(1),AipsError); // pol
403 }
404//
405 start.resize(asap::nAxes);
406 start = 0;
407 start(asap::BeamAxis) = whichBeam;
408 start(asap::IFAxis) = whichIF;
409//
410 end.resize(asap::nAxes);
411 end = shpOut-1;
412 end(asap::BeamAxis) = whichBeam;
413 end(asap::IFAxis) = whichIF;
414}
415
416
417
418// SDFrequenctTable
419
420
421Int SDFrequencyTable::addFrequency(Double refPix, Double refVal, Double inc) {
422 Int idx = -1;
423 Bool addit = False;
424 if (length() > 0) {
425 for (uInt i=0; i< length();++i) {
426 if ( near(refVal,refVal_[i]) ) {
427 if (near(refPix,refPix_[i]) )
428 if ( near(inc,increment_[i]) )
429 idx = Int(i);
430 }
431 }
432 if (idx >= 0) {
433 return idx;
434 }
435 }
436 nFreq_ += 1;
437 refPix_.resize(nFreq_,True);
438 refVal_.resize(nFreq_,True);
439 increment_.resize(nFreq_,True);
440 refPix_[nFreq_-1] = refPix;
441 refVal_[nFreq_-1] = refVal;
442 increment_[nFreq_-1] = inc;
443 idx = nFreq_-1;
444 return idx;
445}
446
447void SDFrequencyTable::addRestFrequency(Double val)
448{
449 if (restFreqs_.nelements() == 0) {
450 restFreqs_.resize(1);
451 restFreqs_[0] = val;
452 } else {
453 Bool found = False;
454 for (uInt i=0;i<restFreqs_.nelements();++i) {
455 if (restFreqs_[i] == val) {
456 found = True;
457 return;
458 }
459 }
460 if (!found) {
461 restFreqs_.resize(restFreqs_.nelements()+1,True);
462 restFreqs_[restFreqs_.nelements()-1] = val;
463 }
464 }
465}
466void SDFrequencyTable::restFrequencies(Vector<Double>& rfs,
467 String& rfunit ) const
468{
469 rfs.resize(restFreqs_.nelements());
470 rfs = restFreqs_;
471 rfunit = restFreqUnit_;
472}
473
474
475
476// SDDataDesc
477
478uInt SDDataDesc::addEntry (const String& source, uInt freqID, const MDirection& dir)
479{
480
481// See if already exists
482
483 if (n_ > 0) {
484 for (uInt i=0; i<n_; i++) {
485 if (source==source_[i] && freqID==freqID_[i]) {
486 return i;
487 }
488 }
489 }
490
491// Not found - add it
492
493 n_ += 1;
494 source_.resize(n_,True);
495 freqID_.resize(n_,True);
496 dir_.resize(n_,True,True);
497//
498 source_[n_-1] = source;
499 freqID_[n_-1] = freqID;
500 dir_[n_-1] = dir;
501//
502 return n_-1;
503}
504
505
506void SDDataDesc::summary() const
507{
508 cerr << "Source FreqID" << endl;
509 for (uInt i=0; i<n_; i++) {
510 cerr << setw(11) << source_(i) << freqID_(i) << endl;
511 }
512}
513
514
Note: See TracBrowser for help on using the repository browser.