source: trunk/src/SDReader.cc@ 50

Last change on this file since 50 was 46, checked in by mmarquar, 20 years ago

Minor edits to verbose info. added empty close function. Switched pseudoHeader pol<->chan

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