source: trunk/src/SDReader.cc@ 277

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