source: trunk/src/SDReader.cc@ 567

Last change on this file since 567 was 477, checked in by kil064, 20 years ago

function convertInstrument moved to SDAttr

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