source: trunk/src/SDReader.cc@ 467

Last change on this file since 467 was 420, checked in by kil064, 20 years ago

handle cross polarization

  • 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
43using namespace casa;
44using namespace asap;
45
46SDReader::SDReader() :
47 reader_(0),
48 header_(0),
49 frequencies_(0),
50 table_(new SDMemTable()),
51 haveXPol_(False)
52{
53 cursor_ = 0;
54}
55SDReader::SDReader(const std::string& filename,
56 int whichIF, int whichBeam) :
57 reader_(0),
58 header_(0),
59 frequencies_(0),
60 table_(new SDMemTable()),
61 haveXPol_(False)
62{
63 cursor_ = 0;
64 open(filename, whichIF, whichBeam);
65}
66
67SDReader::SDReader(CountedPtr<SDMemTable> tbl) :
68 reader_(0),
69 header_(0),
70 table_(tbl),
71 haveXPol_(False)
72{
73 cursor_ = 0;
74}
75
76SDReader::~SDReader() {
77 if (reader_) delete reader_;
78 if (header_) delete header_;
79 if (frequencies_) delete frequencies_;
80}
81
82void SDReader::reset() {
83 cursor_ = 0;
84 table_ = new SDMemTable();
85 open(filename_,ifOffset_, beamOffset_);
86}
87
88
89void SDReader::close() {
90 cerr << "disabled" << endl;
91}
92
93void SDReader::open(const std::string& filename,
94 int whichIF, int whichBeam) {
95 if (reader_) delete reader_; reader_ = 0;
96 Bool haveBase, haveSpectra;
97
98 String inName(filename);
99 Path path(inName);
100 inName = path.expandedName();
101
102 File file(inName);
103 if (!file.exists()) {
104 throw(AipsError("File does not exist"));
105 }
106 filename_ = inName;
107
108 // Create reader and fill in values for arguments
109 String format;
110 Vector<Bool> beams;
111 if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
112 nChan_, nPol_, haveBase, haveSpectra,
113 haveXPol_)) == 0) {
114 throw(AipsError("Creation of PKSreader failed"));
115 }
116 if (!haveSpectra) {
117 delete reader_;
118 reader_ = 0;
119 throw(AipsError("No spectral data in file."));
120 return;
121 }
122
123 nBeam_ = beams.nelements();
124 // Get basic parameters.
125 if (haveXPol_) {
126 cout << "Cross polarization present" << endl;
127 nPol_ += 2; // Convert Complex -> 2 Floats
128 }
129
130 if (header_) delete header_;
131 header_ = new SDHeader();
132 header_->nchan = nChan_;
133 header_->npol = nPol_;
134 header_->nbeam = nBeam_;
135
136 Int status = reader_->getHeader(header_->observer, header_->project,
137 header_->antennaname, header_->antennaposition,
138 header_->obstype,header_->equinox,
139 header_->freqref,
140 header_->utc, header_->reffreq,
141 header_->bandwidth);
142
143 if (status) {
144 delete reader_;
145 reader_ = 0;
146 throw(AipsError("Failed to get header."));
147 return;
148 }
149 if ((header_->obstype).matches("*SW*")) {
150 // need robust way here - probably read ahead of next timestamp
151 cout << "Header indicates frequency switched observation.\n"
152 "setting # of IFs = 1 " << endl;
153 nIF_ = 1;
154 }
155
156 // Determine Telescope and set brightness unit
157
158 Bool throwIt = False;
159 Instrument inst = SDMemTable::convertInstrument(header_->antennaname,
160 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.