source: trunk/src/SDReader.cc@ 284

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

handle environment variable in file name

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