source: trunk/src/SDReader.cc@ 333

Last change on this file since 333 was 332, checked in by kil064, 20 years ago

add IF and beam selection

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