source: trunk/src/SDReader.cc@ 16

Last change on this file since 16 was 16, checked in by mmarquar, 21 years ago

Updated data container. Changed the axis order in the spectrum/flag arrays to [nBeam,nIF,nPol,nChan]
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.6 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDReader.cc: A class to read single dish spectra from SDFITS, RPFITS
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# Malte Marquarding, 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#include <atnf/PKSIO/PKSreader.h>
32#include <aips/Quanta/MVTime.h>
33#include "SDContainer.h"
34#include "SDReader.h"
35
36using namespace atnf_sd;
37
38void SDHeader::print() const {
39 MVTime mvt(this->utc);
40
41 cout << "Observer: " << this->observer << endl
42 << "Project: " << this->project << endl
43 << "Obstype: " << this->obstype << endl
44 << "Antenna: " << this->antennaname << endl
45 << "Ant. Position: " << this->antennaposition << endl
46 << "Equinox: " << this->equinox << endl
47 << "Freq. ref.: " << this->freqref << endl
48 << "Ref. frequency: " << this->reffreq << endl
49 << "Bandwidth: " << this->bandwidth << endl
50 << "Time (utc): "
51 << mvt.string()
52 << endl;
53 //setprecision(10) << this->utc << endl;
54}
55
56
57SDReader::SDReader() :
58 reader_(0),
59 table_(new SDMemTable()) {
60 cursor_ = 0;
61 timestamp_ = 0;
62}
63
64SDReader::SDReader(CountedPtr<SDMemTable> tbl) :
65 reader_(0),
66 table_(tbl) {
67 cursor_ = 0;
68 timestamp_ = 0;
69}
70
71SDReader::~SDReader() {
72 delete reader_;reader_ =0;
73}
74
75void SDReader::reset() {
76 cursor_ = 0;
77 timestamp_ = 0;
78 open(filename_);
79}
80
81void SDReader::open(const std::string& filename) {
82 if (reader_) delete reader_; reader_ = 0;
83 Bool haveBase, haveSpectra, haveXPol;
84 //Int nChan, nIF, nPol;
85 String inName(filename);
86 String format;
87 Vector<Bool> beams;
88 if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
89 nChan_, nPol_, haveBase, haveSpectra,
90 haveXPol)) == 0) {
91 cerr << "PKSreader failed" << endl;
92 }
93 if (!haveSpectra) {
94 delete reader_;
95 reader_ = 0;
96 cerr << "Spectral data absent." << endl;
97 return;
98 }
99 nBeam_ = beams.nelements();
100 cout << "Reading " + format + " format from " + inName << endl;
101 cout << "nChannels = " << nChan_ << ", " << "nPol = " << nPol_ << endl
102 << "nIF = " << nIF_ << endl
103 << "nBeams = " << nBeam_ << endl;
104
105 // Get basic parameters.
106 header_ = SDHeader();
107 header_.nchan = nChan_;
108 header_.npol = nPol_;
109 header_.nbeam = nBeam_;
110
111 //cerr << " getting header ..." << endl;
112 Int status = reader_->getHeader(header_.observer, header_.project,
113 header_.antennaname, header_.antennaposition,
114 header_.obstype,header_.equinox,
115 header_.freqref,
116 header_.utc, header_.reffreq,
117 header_.bandwidth);
118 if (status) {
119 delete reader_;
120 reader_ = 0;
121 cerr << "Failed to get data description." << endl;
122 return;
123 }
124 header_.print();
125 if ((header_.obstype).matches("*SW*")) {
126 // need robust way here - probably read ahead of next timestamp
127 cout << "Header indicates frequency switched observation.\n"
128 "setting # of IFs = 1 " << endl;
129 nIF_ = 1;
130 }
131 header_.nif = nIF_;
132 // Apply selection criteria.
133 cerr << "applying selection criteria..." << endl;
134 Vector<Int> start(nIF_, 1);
135 Vector<Int> end(nIF_, 0);
136 Vector<Int> ref;
137 Vector<Bool> beamSel(nBeam_,True);
138 Vector<Bool> IFsel(nIF_,True);
139 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol);
140 cerr << "open finished" << endl;
141}
142
143int SDReader::read(const std::vector<int>& seq) {
144 cerr << "SDReader::read" << endl;
145 int status = 0;
146
147 Int beamNo, IFno, refBeam, scanNo;
148 Float azimuth, elevation, focusAxi, focusRot, focusTan,
149 humidity, parAngle, pressure, temperature, windAz, windSpeed;
150 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
151 String fieldName, srcName, tcalTime;
152 Vector<Float> calFctr, sigma, tcal, tsys;
153 Matrix<Float> baseLin, baseSub;
154 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
155 Matrix<Float> spectra;
156 Matrix<uChar> flagtra;
157 Complex xCalFctr;
158 Vector<Complex> xPol;
159 mjd = 0;
160 uInt n = seq.size();
161 cerr << header_.nif << ", " << header_.nbeam << endl;
162 uInt stepsize = header_.nif*header_.nbeam;
163 cerr << "SDReader stepsize = " << stepsize << endl;
164 uInt scanid = 0;
165 String prevName = "_noname_";// temporary until scanid is present
166 uInt seqi = 0;
167 //SDFrequencyTable sdft();
168 Bool getAll = False;
169 if (seq[n-1] == -1) getAll = True;
170 while ( (cursor_ <= seq[n-1]) || getAll) {
171 // only needs to be create if in seq
172 SDContainer sc(header_.nbeam,header_.nif,header_.npol,header_.nchan);
173
174 // iterate over one correlator integration cycle = nBeam*nIF
175 for (uInt row=0; row < stepsize; row++) {
176 // add scanid from GROUP field -- this will remove the need for
177 // stepsize as well
178 // spectra(nChan,nPol)!!!
179 status = reader_->read(scanNo, mjd, interval, fieldName, srcName,
180 srcDir, srcPM, srcVel, IFno, refFreq,
181 bandwidth, freqInc, tcal, tcalTime,
182 azimuth, elevation, parAngle, focusAxi,
183 focusTan, focusRot, temperature, pressure,
184 humidity, windSpeed, windAz, refBeam,
185 beamNo, direction, scanRate,
186 tsys, sigma, calFctr, baseLin, baseSub,
187 spectra, flagtra, xCalFctr, xPol);
188 if (status) {
189 if (status == -1) {
190 // EOF.
191 if (row < stepsize-1) cerr << "incomplete integration data." << endl;
192 cerr << "EOF" << endl;
193 return status;
194 }
195 }
196 // if in the given list
197 if (cursor_ == seq[seqi]) {
198 // add integration cycle
199 if (row==0) {
200 //add invariant info: scanNo, mjd, interval, fieldName,
201 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
202 //focusRot, temperature, pressure, humidity, windSpeed,
203 //windAz srcDir, srcPM, srcVel
204 sc.timestamp = mjd;
205 sc.interval = interval;
206 sc.sourcename = srcName;
207 }
208 // add specific info
209 // IFno beamNo are 1-relative
210 // refPix = nChan/2+1 in Integer arith.!
211 //uInt refPix = header_.nchan/2+1;
212 //uInt frqslot = sdft.addFrequency(refPix, refFreq, freqInc);
213
214 if ( srcName != prevName ) {//temp
215 scanid++;//temp
216 prevName = srcName;//temp
217 }//temp
218 sc.scanid = scanid;
219 //sc.setFrequencyMap(frqslot,IFno-1);
220 sc.setSpectrum(spectra, beamNo-1, IFno-1);
221 sc.setFlags(flagtra, beamNo-1, IFno-1);
222 sc.setTsys(tsys, beamNo-1, IFno-1);
223 //sc.addPointing(direction, beamNo-1);
224 }
225 }
226 if (cursor_ == seq[seqi]) {
227 // insert container into table/list
228 table_->putSDContainer(sc);
229 seqi++;// next in list
230 }
231 cursor_++;// increment position in file
232 }
233 return status;
234}
Note: See TracBrowser for help on using the repository browser.