source: trunk/src/SDReader.cc@ 38

Last change on this file since 38 was 33, checked in by mmarquar, 20 years ago

Added reading of frequency information. Also added the getAll flag to get all integrations in the file.

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