source: trunk/src/SDReader.cc@ 409

Last change on this file since 409 was 405, checked in by mar637, 20 years ago

fixed reder rest behaviour. renamed all pyhton_SDReader functions to "_"

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.4 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDReader.cc: A class to read single dish spectra from SDFITS, RPFITS
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#include <casa/iostream.h>
32#include <casa/iomanip.h>
33
34#include <casa/Exceptions.h>
35#include <casa/OS/Path.h>
36#include <casa/OS/File.h>
37#include <casa/Quanta/Unit.h>
38#include <atnf/PKSIO/PKSreader.h>
39
40#include "SDReader.h"
41#include "SDDefs.h"
42
43using namespace casa;
44using namespace asap;
45
46SDReader::SDReader() :
47 reader_(0),
48 header_(0),
49 frequencies_(0),
50 table_(new SDMemTable()) {
51 cursor_ = 0;
52}
53SDReader::SDReader(const std::string& filename,
54 int whichIF, int whichBeam) :
55 reader_(0),
56 header_(0),
57 frequencies_(0),
58 table_(new SDMemTable()) {
59 cursor_ = 0;
60 open(filename, whichIF, whichBeam);
61}
62
63SDReader::SDReader(CountedPtr<SDMemTable> tbl) :
64 reader_(0),
65 header_(0),
66 table_(tbl) {
67 cursor_ = 0;
68}
69
70SDReader::~SDReader() {
71 if (reader_) delete reader_;
72 if (header_) delete header_;
73 if (frequencies_) delete frequencies_;
74}
75
76void SDReader::reset() {
77 cursor_ = 0;
78 table_ = new SDMemTable();
79 open(filename_,ifOffset_, beamOffset_);
80}
81
82
83void SDReader::close() {
84 cerr << "disabled" << endl;
85}
86
87void SDReader::open(const std::string& filename,
88 int whichIF, int whichBeam) {
89 if (reader_) delete reader_; reader_ = 0;
90 Bool haveBase, haveSpectra, haveXPol;
91
92 String inName(filename);
93 Path path(inName);
94 inName = path.expandedName();
95
96 File file(inName);
97 if (!file.exists()) {
98 throw(AipsError("File does not exist"));
99 }
100 filename_ = inName;
101
102 // Create reader and fill in values for arguments
103 String format;
104 Vector<Bool> beams;
105 if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
106 nChan_, nPol_, haveBase, haveSpectra,
107 haveXPol)) == 0) {
108 throw(AipsError("Creation of PKSreader failed"));
109 }
110 if (!haveSpectra) {
111 delete reader_;
112 reader_ = 0;
113 throw(AipsError("No spectral data in file."));
114 return;
115 }
116
117 nBeam_ = beams.nelements();
118 // Get basic parameters.
119 if (haveXPol) {
120 cout << "Warning. Ignoring cross polarisation." << endl;
121 nPol_--;
122 }
123 if (header_) delete header_;
124 header_ = new SDHeader();
125 header_->nchan = nChan_;
126 header_->npol = nPol_;
127 header_->nbeam = nBeam_;
128
129 Int status = reader_->getHeader(header_->observer, header_->project,
130 header_->antennaname, header_->antennaposition,
131 header_->obstype,header_->equinox,
132 header_->freqref,
133 header_->utc, header_->reffreq,
134 header_->bandwidth);
135
136 if (status) {
137 delete reader_;
138 reader_ = 0;
139 throw(AipsError("Failed to get header."));
140 return;
141 }
142 if ((header_->obstype).matches("*SW*")) {
143 // need robust way here - probably read ahead of next timestamp
144 cout << "Header indicates frequency switched observation.\n"
145 "setting # of IFs = 1 " << endl;
146 nIF_ = 1;
147 }
148
149 // Determine Telescope and set brightness unit
150
151 Bool throwIt = False;
152 Instrument inst = SDMemTable::convertInstrument(header_->antennaname,
153 throwIt);
154 header_->fluxunit = "Jy";
155 if (inst==ATMOPRA || inst==TIDBINBILLA) {
156 header_->fluxunit = "K";
157 }
158
159 header_->nif = nIF_;
160 header_->epoch = "UTC";
161
162 // Apply selection criteria.
163
164 Vector<Int> ref;
165 Vector<Bool> beamSel(nBeam_,True);
166 Vector<Bool> IFsel(nIF_,True);
167
168 ifOffset_ = 0;
169 if (whichIF>=0) {
170 if (whichIF>=0 && whichIF<nIF_) {
171 IFsel = False;
172 IFsel(whichIF) = True;
173 header_->nif = 1;
174 nIF_ = 1;
175 ifOffset_ = whichIF;
176 } else {
177 throw(AipsError("Illegal IF selection"));
178 }
179 }
180
181 beamOffset_ = 0;
182 if (whichBeam>=0) {
183 if (whichBeam>=0 && whichBeam<nBeam_) {
184 beamSel = False;
185 beamSel(whichBeam) = True;
186 header_->nbeam = 1;
187 nBeam_ = 1;
188 beamOffset_ = whichBeam;
189 } else {
190 throw(AipsError("Illegal Beam selection"));
191 }
192 }
193 Vector<Int> start(nIF_, 1);
194 Vector<Int> end(nIF_, 0);
195 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol);
196 table_->putSDHeader(*header_);
197
198 if (frequencies_) delete frequencies_;
199 frequencies_ = new SDFrequencyTable();
200 frequencies_->setRefFrame(header_->freqref);
201 frequencies_->setRestFrequencyUnit("Hz");
202 frequencies_->setEquinox(header_->equinox);
203}
204
205int SDReader::read(const std::vector<int>& seq) {
206 int status = 0;
207
208 Int beamNo, IFno, refBeam, scanNo, cycleNo;
209 Float azimuth, elevation, focusAxi, focusRot, focusTan,
210 humidity, parAngle, pressure, temperature, windAz, windSpeed;
211 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
212 String fieldName, srcName, tcalTime;
213 Vector<Float> calFctr, sigma, tcal, tsys;
214 Matrix<Float> baseLin, baseSub;
215 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
216 Matrix<Float> spectra;
217 Matrix<uChar> flagtra;
218 Complex xCalFctr;
219 Vector<Complex> xPol;
220 uInt n = seq.size();
221
222 uInt stepsize = header_->nif*header_->nbeam;
223 uInt seqi = 0;
224 Bool getAll = False;
225
226 if (seq[n-1] == -1) getAll = True;
227 while ( (cursor_ <= seq[n-1]) || getAll) {
228 // only needs to be create if in seq
229 SDContainer sc(header_->nbeam,header_->nif,header_->npol,header_->nchan);
230 // iterate over one correlator integration cycle = nBeam*nIF
231 for (uInt row=0; row < stepsize; row++) {
232 // stepsize as well
233 // spectra(nChan,nPol)!!!
234 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
235 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
236 bandwidth, freqInc, restFreq, tcal, tcalTime,
237 azimuth, elevation, parAngle, focusAxi,
238 focusTan, focusRot, temperature, pressure,
239 humidity, windSpeed, windAz, refBeam,
240 beamNo, direction, scanRate,
241 tsys, sigma, calFctr, baseLin, baseSub,
242 spectra, flagtra, xCalFctr, xPol);
243
244 // Make sure beam/IF numbers are 0-relative - dealing with
245 // possible IF or Beam selection
246 beamNo = beamNo - beamOffset_ - 1;
247 IFno = IFno - ifOffset_ - 1;
248
249 if (status) {
250 if (status == -1) {
251 // EOF.
252 if (row > 0 && row < stepsize-1)
253 cerr << "incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan are present." << endl;
254
255 // flush frequency table
256 table_->putSDFreqTable(*frequencies_);
257 return status;
258 }
259 }
260 // if in the given list
261 if (cursor_ == seq[seqi] || getAll) {
262 // add integration cycle
263 if (row==0) {
264 //add invariant info: scanNo, mjd, interval, fieldName,
265 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
266 //focusRot, temperature, pressure, humidity, windSpeed,
267 //windAz srcDir, srcPM, srcVel
268 sc.timestamp = mjd;
269 sc.interval = interval;
270 sc.sourcename = srcName;
271 sc.fieldname = fieldName;
272 sc.azimuth = azimuth;
273 sc.elevation = elevation;
274 }
275 // add specific info
276 // refPix = nChan/2+1 in 1-rel Integer arith.!
277 Int refPix = header_->nchan/2; // 0-rel
278 uInt freqID = frequencies_->addFrequency(refPix, refFreq, freqInc);
279 uInt restFreqID = frequencies_->addRestFrequency(restFreq);
280
281 sc.setFrequencyMap(freqID, IFno);
282 sc.setRestFrequencyMap(restFreqID, IFno);
283
284 sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1];
285 sc.tcaltime = tcalTime;
286 sc.parangle = parAngle;
287 sc.refbeam = refBeam-1;//make it 0-based; -1 if nbeams == 1
288 sc.scanid = scanNo-1;//make it 0-based
289 sc.setSpectrum(spectra, beamNo, IFno);
290 sc.setFlags(flagtra, beamNo, IFno);
291 sc.setTsys(tsys, beamNo, IFno);
292 sc.setDirection(direction, beamNo);
293 }
294 }
295
296 if (cursor_ == seq[seqi] || getAll) {
297 // insert container into table/list
298 table_->putSDContainer(sc);
299 seqi++;// next in list
300 }
301 cursor_++;// increment position in file
302
303 }
304 table_->putSDFreqTable(*frequencies_);
305 return status;
306}
Note: See TracBrowser for help on using the repository browser.