source: trunk/src/SDContainer.cc@ 402

Last change on this file since 402 was 396, checked in by kil064, 20 years ago

add a secondary index argument to the SDDataDesc container

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