source: trunk/src/SDReader.cc@ 402

Last change on this file since 402 was 390, checked in by kil064, 20 years ago

support rest freqs being a column in SDmemTable

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.4 KB
RevLine 
[2]1//#---------------------------------------------------------------------------
2//# SDReader.cc: A class to read single dish spectra from SDFITS, RPFITS
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
[125]5//# ATNF
[2]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//#---------------------------------------------------------------------------
[125]31#include <casa/Exceptions.h>
[284]32#include <casa/OS/Path.h>
[290]33#include <casa/OS/File.h>
[338]34#include <casa/Quanta/Unit.h>
[2]35#include <atnf/PKSIO/PKSreader.h>
[125]36
[2]37#include "SDReader.h"
[337]38#include "SDDefs.h"
[2]39
[125]40using namespace casa;
[83]41using namespace asap;
[2]42
43SDReader::SDReader() :
44 reader_(0),
45 table_(new SDMemTable()) {
46 cursor_ = 0;
47}
[342]48SDReader::SDReader(const std::string& filename,
[338]49 int whichIF, int whichBeam) :
[46]50 reader_(0),
51 table_(new SDMemTable()) {
52 cursor_ = 0;
[342]53 open(filename, whichIF, whichBeam);
[46]54}
[2]55
56SDReader::SDReader(CountedPtr<SDMemTable> tbl) :
57 reader_(0),
58 table_(tbl) {
59 cursor_ = 0;
60}
61
62SDReader::~SDReader() {
[46]63 if (reader_) delete reader_; reader_ = 0;
[2]64}
65
66void SDReader::reset() {
[46]67 cursor_ = 0;
[87]68 table_ = new SDMemTable();
69 if (reader_) delete reader_; reader_ = 0;
70 // re-create the reader
71 Bool haveBase, haveSpectra, haveXPol;
72 String format;
73 Vector<Bool> beams;
74 if ((reader_ = getPKSreader(filename_, 0, False, format, beams, nIF_,
75 nChan_, nPol_, haveBase, haveSpectra,
76 haveXPol)) == 0) {
77 throw(AipsError("PKSreader failed"));
78 }
79 // re-enter the header
80 table_->putSDHeader(header_);
[2]81}
82
[332]83
[46]84void SDReader::close() {
85 cerr << "disabled" << endl;
86}
87
[332]88
89
[342]90void SDReader::open(const std::string& filename,
[338]91 int whichIF, int whichBeam) {
[2]92 if (reader_) delete reader_; reader_ = 0;
93 Bool haveBase, haveSpectra, haveXPol;
[284]94//
[2]95 String inName(filename);
[284]96 Path path(inName);
97 inName = path.expandedName();
[290]98//
99 File file(inName);
100 if (!file.exists()) {
101 throw(AipsError("File does not exist"));
102 }
[87]103 filename_ = inName;
[332]104
105// Create reader and fill in values for arguments
106
[2]107 String format;
108 Vector<Bool> beams;
109 if ((reader_ = getPKSreader(inName, 0, False, format, beams, nIF_,
110 nChan_, nPol_, haveBase, haveSpectra,
111 haveXPol)) == 0) {
[87]112 throw(AipsError("PKSreader failed"));
[2]113 }
114 if (!haveSpectra) {
115 delete reader_;
116 reader_ = 0;
[87]117 throw(AipsError("No spectral data in file."));
[2]118 return;
119 }
120 nBeam_ = beams.nelements();
121 // Get basic parameters.
[110]122 if (haveXPol) {
123 cout << "Warning. Ignoring cross polarisation." << endl;
124 nPol_--;
125 }
[2]126 header_ = SDHeader();
127 header_.nchan = nChan_;
128 header_.npol = nPol_;
129 header_.nbeam = nBeam_;
[87]130 Int status = reader_->getHeader(header_.observer, header_.project,
131 header_.antennaname, header_.antennaposition,
132 header_.obstype,header_.equinox,
133 header_.freqref,
134 header_.utc, header_.reffreq,
[2]135 header_.bandwidth);
136 if (status) {
137 delete reader_;
138 reader_ = 0;
[87]139 throw(AipsError("Failed to get header."));
[2]140 return;
141 }
142 if ((header_.obstype).matches("*SW*")) {
143 // need robust way here - probably read ahead of next timestamp
144 cout << "Header indicates frequency switched observation.\n"
145 "setting # of IFs = 1 " << endl;
146 nIF_ = 1;
147 }
[337]148
[342]149// Determine Telescope and set brightness unit
150
151 Bool throwIt = False;
152 Instrument inst = SDMemTable::convertInstrument (header_.antennaname, throwIt);
153 header_.fluxunit = "Jy";
154 if (inst==ATMOPRA || inst==TIDBINBILLA) {
155 header_.fluxunit = "K";
156 }
[337]157//
[2]158 header_.nif = nIF_;
[205]159 header_.epoch = "UTC";
[332]160
[2]161 // Apply selection criteria.
[332]162
[2]163 Vector<Int> ref;
164 Vector<Bool> beamSel(nBeam_,True);
165 Vector<Bool> IFsel(nIF_,True);
[332]166//
167 ifOffset_ = 0;
168 if (whichIF>=0) {
169 if (whichIF>=0 && whichIF<nIF_) {
170 IFsel = False;
171 IFsel(whichIF) = True;
172 header_.nif = 1;
173 nIF_ = 1;
174 ifOffset_ = whichIF;
175 } else {
176 throw(AipsError("Illegal IF selection"));
177 }
178 }
179//
180 beamOffset_ = 0;
181 if (whichBeam>=0) {
182 if (whichBeam>=0 && whichBeam<nBeam_) {
183 beamSel = False;
184 beamSel(whichBeam) = True;
185 header_.nbeam = 1;
186 nBeam_ = 1;
187 beamOffset_ = whichBeam;
188 } else {
189 throw(AipsError("Illegal Beam selection"));
190 }
191 }
192//
193 Vector<Int> start(nIF_, 1);
194 Vector<Int> end(nIF_, 0);
[2]195 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol);
[18]196 table_->putSDHeader(header_);
[33]197 frequencies_.setRefFrame(header_.freqref);
[205]198 frequencies_.setRestFrequencyUnit("Hz");
[33]199 frequencies_.setEquinox(header_.equinox);
[2]200}
201
202int SDReader::read(const std::vector<int>& seq) {
[87]203 int status = 0;
204
[17]205 Int beamNo, IFno, refBeam, scanNo, cycleNo;
[87]206 Float azimuth, elevation, focusAxi, focusRot, focusTan,
[2]207 humidity, parAngle, pressure, temperature, windAz, windSpeed;
[73]208 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
[2]209 String fieldName, srcName, tcalTime;
210 Vector<Float> calFctr, sigma, tcal, tsys;
211 Matrix<Float> baseLin, baseSub;
212 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
213 Matrix<Float> spectra;
214 Matrix<uChar> flagtra;
215 Complex xCalFctr;
216 Vector<Complex> xPol;
217 uInt n = seq.size();
[25]218
[2]219 uInt stepsize = header_.nif*header_.nbeam;
220 uInt seqi = 0;
221 Bool getAll = False;
222 if (seq[n-1] == -1) getAll = True;
223 while ( (cursor_ <= seq[n-1]) || getAll) {
224 // only needs to be create if in seq
[16]225 SDContainer sc(header_.nbeam,header_.nif,header_.npol,header_.nchan);
[2]226 // iterate over one correlator integration cycle = nBeam*nIF
227 for (uInt row=0; row < stepsize; row++) {
228 // stepsize as well
229 // spectra(nChan,nPol)!!!
[87]230 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
231 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
232 bandwidth, freqInc, restFreq, tcal, tcalTime,
233 azimuth, elevation, parAngle, focusAxi,
234 focusTan, focusRot, temperature, pressure,
235 humidity, windSpeed, windAz, refBeam,
236 beamNo, direction, scanRate,
237 tsys, sigma, calFctr, baseLin, baseSub,
238 spectra, flagtra, xCalFctr, xPol);
[332]239
[337]240// Make sure beam/IF numbers are 0-relative - dealing with possible IF or Beam selection
[332]241
242 beamNo = beamNo - beamOffset_ - 1;
243 IFno = IFno - ifOffset_ - 1;
244//
[2]245 if (status) {
[87]246 if (status == -1) {
247 // EOF.
[205]248 if (row > 0 && row < stepsize-1)
249 cerr << "incomplete scan data.\n Probably means not all Beams/IFs/Pols within a scan are present." << endl;
[87]250 //cerr << "EOF" << endl;
251 table_->putSDFreqTable(frequencies_);
252 return status;
253 }
254 }
[2]255 // if in the given list
[33]256 if (cursor_ == seq[seqi] || getAll) {
[87]257 // add integration cycle
258 if (row==0) {
259 //add invariant info: scanNo, mjd, interval, fieldName,
260 //srcName, azimuth, elevation, parAngle, focusAxi, focusTan,
261 //focusRot, temperature, pressure, humidity, windSpeed,
262 //windAz srcDir, srcPM, srcVel
263 sc.timestamp = mjd;
264 sc.interval = interval;
265 sc.sourcename = srcName;
[110]266 sc.fieldname = fieldName;
267 sc.azimuth = azimuth;
268 sc.elevation = elevation;
[87]269 }
270 // add specific info
271 // refPix = nChan/2+1 in Integer arith.!
272 Int refPix = header_.nchan/2+1;
[390]273 uInt freqID = frequencies_.addFrequency(refPix, refFreq, freqInc);
274 uInt restFreqID = frequencies_.addRestFrequency(restFreq);
[332]275//
[390]276 sc.setFrequencyMap(freqID, IFno);
277 sc.setRestFrequencyMap(restFreqID, IFno);
[110]278 sc.tcal[0] = tcal[0];sc.tcal[1] = tcal[1];
279 sc.tcaltime = tcalTime;
280 sc.parangle = parAngle;
[374]281 sc.refbeam = refBeam-1;//make it 0-based; -1 if nbeams == 1
[87]282 sc.scanid = scanNo-1;//make it 0-based
[332]283 sc.setSpectrum(spectra, beamNo, IFno);
284 sc.setFlags(flagtra, beamNo, IFno);
285 sc.setTsys(tsys, beamNo, IFno);
286 sc.setDirection(direction, beamNo);
[2]287 }
288 }
[33]289 if (cursor_ == seq[seqi] || getAll) {
[2]290 // insert container into table/list
291 table_->putSDContainer(sc);
292 seqi++;// next in list
293 }
[16]294 cursor_++;// increment position in file
[2]295 }
[33]296 table_->putSDFreqTable(frequencies_);
[2]297 return status;
298}
Note: See TracBrowser for help on using the repository browser.