source: branches/Release12/src/SDReader.cc@ 2679

Last change on this file since 2679 was 754, checked in by mar637, 19 years ago

Bug fix: SDFreqTable did not copy user settings for unit and frame

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.0 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 pushLog("Cross polarization present");
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 pushLog("Header indicates frequency switched observation.\n"
153 "setting # of IFs = 1 ");
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_->setBaseRefFrame(header_->freqref);
209 frequencies_->setRestFrequencyUnit("Hz");
210 frequencies_->setEquinox(header_->equinox);
211}
212
213int SDReader::read(const std::vector<int>& seq) {
214 int status = 0;
215
216 Int beamNo, IFno, refBeam, scanNo, cycleNo;
217 Float azimuth, elevation, focusAxi, focusRot, focusTan,
218 humidity, parAngle, pressure, temperature, windAz, windSpeed;
219 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
220 String fieldName, srcName, tcalTime;
221 Vector<Float> calFctr, sigma, tcal, tsys;
222 Matrix<Float> baseLin, baseSub;
223 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
224 Matrix<Float> spectra;
225 Matrix<uChar> flagtra;
226 Complex xCalFctr;
227 Vector<Complex> xPol;
228 uInt n = seq.size();
229
230
231 uInt stepsize = header_->nif*header_->nbeam;
232 uInt seqi = 0;
233 Bool getAll = False;
234
235 if (seq[n-1] == -1) getAll = True;
236 while ( (cursor_ <= seq[n-1]) || getAll) {
237 // only needs to be create if in seq
238 SDContainer sc(header_->nbeam,header_->nif,header_->npol,header_->nchan);
239 // iterate over one correlator integration cycle = nBeam*nIF
240 for (uInt row=0; row < stepsize; row++) {
241 // stepsize as well
242 // spectra(nChan,nPol)!!!
243 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
244 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
245 bandwidth, freqInc, restFreq, tcal, tcalTime,
246 azimuth, elevation, parAngle, focusAxi,
247 focusTan, focusRot, temperature, pressure,
248 humidity, windSpeed, windAz, refBeam,
249 beamNo, direction, scanRate,
250 tsys, sigma, calFctr, baseLin, baseSub,
251 spectra, flagtra, xCalFctr, xPol);
252
253 // Make sure beam/IF numbers are 0-relative - dealing with
254 // possible IF or Beam selection
255 beamNo = beamNo - beamOffset_ - 1;
256 IFno = IFno - ifOffset_ - 1;
257
258 if (status) {
259 if (status == -1) {
260 // EOF.
261 if (row > 0 && row < stepsize-1)
262 pushLog("incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan are present.");
263
264 // flush frequency table
265 table_->putSDFreqTable(*frequencies_);
266 return status;
267 }
268 }
269 // if in the given list
270 if (cursor_ == seq[seqi] || getAll) {
271 // add integration cycle
272 if (row==0) {
273 //add invariant info: scanNo, mjd, interval, fieldName,
274 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
275 //focusRot, temperature, pressure, humidity, windSpeed,
276 //windAz srcDir, srcPM, srcVel
277 sc.timestamp = mjd;
278 sc.interval = interval;
279 sc.sourcename = srcName;
280 sc.fieldname = fieldName;
281 sc.azimuth = azimuth;
282 sc.elevation = elevation;
283 }
284 // add specific info
285 // refPix = nChan/2+1 in 1-rel Integer arith.!
286 Int refPix = header_->nchan/2; // 0-rel
287 uInt freqID = frequencies_->addFrequency(refPix, refFreq, freqInc);
288 uInt restFreqID = frequencies_->addRestFrequency(restFreq);
289
290 sc.setFrequencyMap(freqID, IFno);
291 sc.setRestFrequencyMap(restFreqID, IFno);
292
293 sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1];
294 sc.tcaltime = tcalTime;
295 sc.parangle = parAngle;
296 sc.refbeam = -1; //nbeams == 1
297 if (nBeam_ > 1) // circumvent a bug "asap0000" in read which
298 // returns a random refbema number on multiple
299 // reads
300 sc.refbeam = refBeam-1;//make it 0-based;
301 sc.scanid = scanNo-1;//make it 0-based
302 if (haveXPol_) {
303 sc.setSpectrum(spectra, xPol, beamNo, IFno);
304 sc.setFlags(flagtra, beamNo, IFno, True);
305 } else {
306 sc.setSpectrum(spectra, beamNo, IFno);
307 sc.setFlags(flagtra, beamNo, IFno, False);
308 }
309 sc.setTsys(tsys, beamNo, IFno, haveXPol_);
310 sc.setDirection(direction, beamNo);
311 }
312 }
313
314 if (cursor_ == seq[seqi] || getAll) {
315 // insert container into table/list
316 table_->putSDContainer(sc);
317 seqi++;// next in list
318 }
319 cursor_++;// increment position in file
320
321 }
322 table_->putSDFreqTable(*frequencies_);
323 return status;
324}
Note: See TracBrowser for help on using the repository browser.