source: trunk/src/SDReader.cc@ 181

Last change on this file since 181 was 125, checked in by mar637, 20 years ago

Moved to casa namespace.
Adjusted the copyright to be ATNF.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 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 // Apply selection criteria.
130 Vector<Int> start(nIF_, 1);
131 Vector<Int> end(nIF_, 0);
132 Vector<Int> ref;
133 Vector<Bool> beamSel(nBeam_,True);
134 Vector<Bool> IFsel(nIF_,True);
135 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol);
136 table_->putSDHeader(header_);
137 frequencies_.setRefFrame(header_.freqref);
138 frequencies_.setEquinox(header_.equinox);
139}
140
141int SDReader::read(const std::vector<int>& seq) {
142 int status = 0;
143
144 Int beamNo, IFno, refBeam, scanNo, cycleNo;
145 Float azimuth, elevation, focusAxi, focusRot, focusTan,
146 humidity, parAngle, pressure, temperature, windAz, windSpeed;
147 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
148 String fieldName, srcName, tcalTime;
149 Vector<Float> calFctr, sigma, tcal, tsys;
150 Matrix<Float> baseLin, baseSub;
151 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
152 Matrix<Float> spectra;
153 Matrix<uChar> flagtra;
154 Complex xCalFctr;
155 Vector<Complex> xPol;
156 uInt n = seq.size();
157
158 uInt stepsize = header_.nif*header_.nbeam;
159 //cerr << "SDReader stepsize = " << stepsize << endl;
160 uInt seqi = 0;
161 Bool getAll = False;
162 if (seq[n-1] == -1) getAll = True;
163 while ( (cursor_ <= seq[n-1]) || getAll) {
164 // only needs to be create if in seq
165 SDContainer sc(header_.nbeam,header_.nif,header_.npol,header_.nchan);
166
167 // iterate over one correlator integration cycle = nBeam*nIF
168 for (uInt row=0; row < stepsize; row++) {
169 // stepsize as well
170 // spectra(nChan,nPol)!!!
171 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
172 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
173 bandwidth, freqInc, restFreq, tcal, tcalTime,
174 azimuth, elevation, parAngle, focusAxi,
175 focusTan, focusRot, temperature, pressure,
176 humidity, windSpeed, windAz, refBeam,
177 beamNo, direction, scanRate,
178 tsys, sigma, calFctr, baseLin, baseSub,
179 spectra, flagtra, xCalFctr, xPol);
180 if (status) {
181 if (status == -1) {
182 // EOF.
183 if (row < stepsize-1) cerr << "incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan a present." << endl;
184 //cerr << "EOF" << endl;
185 table_->putSDFreqTable(frequencies_);
186 return status;
187 }
188 }
189 // if in the given list
190 if (cursor_ == seq[seqi] || getAll) {
191 // add integration cycle
192 if (row==0) {
193 //add invariant info: scanNo, mjd, interval, fieldName,
194 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
195 //focusRot, temperature, pressure, humidity, windSpeed,
196 //windAz srcDir, srcPM, srcVel
197 sc.timestamp = mjd;
198 sc.interval = interval;
199 sc.sourcename = srcName;
200 sc.fieldname = fieldName;
201 sc.azimuth = azimuth;
202 sc.elevation = elevation;
203 }
204 // add specific info
205 // IFno beamNo are 1-relative
206 // refPix = nChan/2+1 in Integer arith.!
207 Int refPix = header_.nchan/2+1;
208 Int frqslot = frequencies_.addFrequency(refPix, refFreq, freqInc);
209 sc.setFrequencyMap(frqslot,IFno-1);
210 sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1];
211 sc.tcaltime = tcalTime;
212 sc.parangle = parAngle;
213 sc.refbeam = refBeam;
214 sc.scanid = scanNo-1;//make it 0-based
215 sc.setSpectrum(spectra, beamNo-1, IFno-1);
216 sc.setFlags(flagtra, beamNo-1, IFno-1);
217 sc.setTsys(tsys, beamNo-1, IFno-1);
218 sc.setDirection(direction, beamNo-1);
219 }
220 }
221 if (cursor_ == seq[seqi] || getAll) {
222 // insert container into table/list
223 table_->putSDContainer(sc);
224 seqi++;// next in list
225 }
226 cursor_++;// increment position in file
227 }
228 table_->putSDFreqTable(frequencies_);
229 return status;
230}
Note: See TracBrowser for help on using the repository browser.