source: branches/atnf/src/SDReader.cc@ 1510

Last change on this file since 1510 was 2, checked in by mmarquar, 20 years ago

Initial revision

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