source: trunk/src/SDReader.cc@ 21

Last change on this file since 21 was 18, checked in by mmarquar, 20 years ago

Moved SDHeader from SDReader to SDConatiner. Added header to SDMemTable. Added access funtions to nif,nbema,npol,nchan

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.8 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 //Int nChan, nIF, nPol;
65 String inName(filename);
66 String format;
67 Vector<Bool> beams;
68 if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
69 nChan_, nPol_, haveBase, haveSpectra,
70 haveXPol)) == 0) {
71 cerr << "PKSreader failed" << endl;
72 }
73 if (!haveSpectra) {
74 delete reader_;
75 reader_ = 0;
76 cerr << "Spectral data absent." << endl;
77 return;
78 }
79 nBeam_ = beams.nelements();
80 //cout << "Reading " + format + " format from " + inName << endl;
81 //cout << "nChannels = " << nChan_ << ", " << "nPol = " << nPol_ << endl
82 // << "nIF = " << nIF_ << endl
83 // << "nBeams = " << nBeam_ << endl;
84
85 // Get basic parameters.
86 header_ = SDHeader();
87 header_.nchan = nChan_;
88 header_.npol = nPol_;
89 header_.nbeam = nBeam_;
90
91 //cerr << " getting header ..." << endl;
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}
121
122int SDReader::read(const std::vector<int>& seq) {
123 int status = 0;
124
125 Int beamNo, IFno, refBeam, scanNo, cycleNo;
126 Float azimuth, elevation, focusAxi, focusRot, focusTan,
127 humidity, parAngle, pressure, temperature, windAz, windSpeed;
128 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
129 String fieldName, srcName, tcalTime;
130 Vector<Float> calFctr, sigma, tcal, tsys;
131 Matrix<Float> baseLin, baseSub;
132 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
133 Matrix<Float> spectra;
134 Matrix<uChar> flagtra;
135 Complex xCalFctr;
136 Vector<Complex> xPol;
137 mjd = 0;
138 uInt n = seq.size();
139 //cerr << header_.nif << ", " << header_.nbeam << endl;
140 uInt stepsize = header_.nif*header_.nbeam;
141 cerr << "SDReader stepsize = " << stepsize << endl;
142 uInt scanid = 0;
143 String prevName = "_noname_";// temporary until scanid is present
144 uInt seqi = 0;
145 //SDFrequencyTable sdft();
146 Bool getAll = False;
147 if (seq[n-1] == -1) getAll = True;
148 while ( (cursor_ <= seq[n-1]) || getAll) {
149 // only needs to be create if in seq
150 SDContainer sc(header_.nbeam,header_.nif,header_.npol,header_.nchan);
151
152 // iterate over one correlator integration cycle = nBeam*nIF
153 for (uInt row=0; row < stepsize; row++) {
154 // stepsize as well
155 // spectra(nChan,nPol)!!!
156 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
157 srcName,
158 srcDir, srcPM, srcVel, IFno, refFreq,
159 bandwidth, freqInc, tcal, tcalTime,
160 azimuth, elevation, parAngle, focusAxi,
161 focusTan, focusRot, temperature, pressure,
162 humidity, windSpeed, windAz, refBeam,
163 beamNo, direction, scanRate,
164 tsys, sigma, calFctr, baseLin, baseSub,
165 spectra, flagtra, xCalFctr, xPol);
166 if (status) {
167 if (status == -1) {
168 // EOF.
169 if (row < stepsize-1) cerr << "incomplete integration data." << endl;
170 cerr << "EOF" << endl;
171 return status;
172 }
173 }
174 // if in the given list
175 if (cursor_ == seq[seqi]) {
176 // add integration cycle
177 if (row==0) {
178 //add invariant info: scanNo, mjd, interval, fieldName,
179 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
180 //focusRot, temperature, pressure, humidity, windSpeed,
181 //windAz srcDir, srcPM, srcVel
182 sc.timestamp = mjd;
183 sc.interval = interval;
184 sc.sourcename = srcName;
185 }
186 // add specific info
187 // IFno beamNo are 1-relative
188 // refPix = nChan/2+1 in Integer arith.!
189 //uInt refPix = header_.nchan/2+1;
190 //uInt frqslot = sdft.addFrequency(refPix, refFreq, freqInc);
191
192 sc.scanid = scanNo;
193 //sc.setFrequencyMap(frqslot,IFno-1);
194 sc.setSpectrum(spectra, beamNo-1, IFno-1);
195 sc.setFlags(flagtra, beamNo-1, IFno-1);
196 sc.setTsys(tsys, beamNo-1, IFno-1);
197 //sc.addPointing(direction, beamNo-1);
198 }
199 }
200 if (cursor_ == seq[seqi]) {
201 // insert container into table/list
202 table_->putSDContainer(sc);
203 seqi++;// next in list
204 }
205 cursor_++;// increment position in file
206 }
207 return status;
208}
Note: See TracBrowser for help on using the repository browser.