source: trunk/src/STFiller.cpp@ 849

Last change on this file since 849 was 847, checked in by mar637, 19 years ago

numerous changes before move to new svn repository sourcecode.atnf.csiro.au
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.4 KB
Line 
1//
2// C++ Implementation: STFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iostream.h>
13#include <casa/iomanip.h>
14
15#include <casa/Exceptions.h>
16#include <casa/OS/Path.h>
17#include <casa/OS/File.h>
18#include <casa/Quanta/Unit.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Utilities/Regex.h>
21
22#include <casa/Containers/RecordField.h>
23
24#include <tables/Tables/TableRow.h>
25
26#include <atnf/PKSIO/PKSreader.h>
27
28#include "STDefs.h"
29#include "SDAttr.h"
30
31#include "STFiller.h"
32
33using namespace casa;
34
35namespace asap {
36
37STFiller::STFiller() :
38 reader_(0),
39 header_(0),
40 table_(0)
41{
42}
43
44STFiller::STFiller( CountedPtr< Scantable > stbl ) :
45 reader_(0),
46 header_(0),
47 table_(stbl)
48{
49}
50
51STFiller::STFiller(const std::string & filename,
52 int whichIF, int whichBeam ) :
53 reader_(0),
54 header_(0),
55 table_(0)
56{
57 open(filename, whichIF, whichBeam);
58}
59
60STFiller::~STFiller()
61{
62 close();
63}
64
65void STFiller::open( const std::string & filename,
66 int whichIF, int whichBeam )
67{
68 if (table_.null()) {
69 table_ = new Scantable();
70 }
71 if (reader_) { delete reader_; reader_ = 0; }
72 Bool haveBase, haveSpectra;
73
74 String inName(filename);
75 Path path(inName);
76 inName = path.expandedName();
77
78 File file(inName);
79 if ( !file.exists() ) {
80 throw(AipsError("File does not exist"));
81 }
82 filename_ = inName;
83
84 // Create reader and fill in values for arguments
85 String format;
86 Vector<Bool> beams;
87 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, nIF_,
88 nChan_, nPol_, haveBase, haveSpectra,
89 haveXPol_)) == 0 ) {
90 throw(AipsError("Creation of PKSreader failed"));
91 }
92 if (!haveSpectra) {
93 delete reader_;
94 reader_ = 0;
95 throw(AipsError("No spectral data in file."));
96 return;
97 }
98
99 nBeam_ = beams.nelements();
100 // Get basic parameters.
101 if ( haveXPol_ ) {
102 pushLog("Cross polarization present");
103 nPol_ += 2; // Convert Complex -> 2 Floats
104 }
105
106 if (header_) delete header_;
107 header_ = new SDHeader();
108 header_->nchan = nChan_;
109 header_->npol = nPol_;
110 header_->nbeam = nBeam_;
111
112 Int status = reader_->getHeader(header_->observer, header_->project,
113 header_->antennaname, header_->antennaposition,
114 header_->obstype,header_->equinox,
115 header_->freqref,
116 header_->utc, header_->reffreq,
117 header_->bandwidth);
118
119 if (status) {
120 delete reader_;
121 reader_ = 0;
122 delete header_;
123 header_ = 0;
124 throw(AipsError("Failed to get header."));
125 }
126 if ((header_->obstype).matches("*SW*")) {
127 // need robust way here - probably read ahead of next timestamp
128 pushLog("Header indicates frequency switched observation.\n"
129 "setting # of IFs = 1 ");
130 nIF_ = 1;
131 header_->obstype = String("fswitch");
132 }
133
134 // Determine Telescope and set brightness unit
135
136 Bool throwIt = False;
137 Instrument inst = SDAttr::convertInstrument(header_->antennaname, throwIt);
138 header_->fluxunit = "Jy";
139 if (inst==ATMOPRA || inst==TIDBINBILLA) {
140 header_->fluxunit = "K";
141 }
142
143 header_->nif = nIF_;
144 header_->epoch = "UTC";
145 // *** header_->frequnit = "Hz"
146 // Apply selection criteria.
147
148 Vector<Int> ref;
149 Vector<Bool> beamSel(nBeam_,True);
150 Vector<Bool> IFsel(nIF_,True);
151
152 ifOffset_ = 0;
153 if (whichIF>=0) {
154 if (whichIF>=0 && whichIF<nIF_) {
155 IFsel = False;
156 IFsel(whichIF) = True;
157 header_->nif = 1;
158 nIF_ = 1;
159 ifOffset_ = whichIF;
160 } else {
161 delete reader_;
162 reader_ = 0;
163 delete header_;
164 header_ = 0;
165 throw(AipsError("Illegal IF selection"));
166 }
167 }
168
169 beamOffset_ = 0;
170 if (whichBeam>=0) {
171 if (whichBeam>=0 && whichBeam<nBeam_) {
172 beamSel = False;
173 beamSel(whichBeam) = True;
174 header_->nbeam = 1;
175 nBeam_ = 1;
176 beamOffset_ = whichBeam;
177 } else {
178 delete reader_;
179 reader_ = 0;
180 delete header_;
181 header_ = 0;
182 throw(AipsError("Illegal Beam selection"));
183 }
184 }
185 Vector<Int> start(nIF_, 1);
186 Vector<Int> end(nIF_, 0);
187 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol_);
188 table_->setHeader(*header_);
189}
190
191void STFiller::close( )
192{
193 delete reader_;reader_=0;
194 delete header_;header_=0;
195 table_ = 0;
196}
197
198int asap::STFiller::read( )
199{
200 int status = 0;
201
202 Int beamNo, IFno, refBeam, scanNo, cycleNo;
203 Float azimuth, elevation, focusAxi, focusRot, focusTan,
204 humidity, parAngle, pressure, temperature, windAz, windSpeed;
205 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
206 String fieldName, srcName, tcalTime;
207 Vector<Float> calFctr, sigma, tcal, tsys;
208 Matrix<Float> baseLin, baseSub;
209 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
210 Matrix<Float> spectra;
211 Matrix<uChar> flagtra;
212 Complex xCalFctr;
213 Vector<Complex> xPol;
214 while ( status == 0 ) {
215 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
216 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
217 bandwidth, freqInc, restFreq, tcal, tcalTime,
218 azimuth, elevation, parAngle, focusAxi,
219 focusTan, focusRot, temperature, pressure,
220 humidity, windSpeed, windAz, refBeam,
221 beamNo, direction, scanRate,
222 tsys, sigma, calFctr, baseLin, baseSub,
223 spectra, flagtra, xCalFctr, xPol);
224 if ( status != 0 ) break;
225 TableRow row(table_->table());
226 TableRecord& rec = row.record();
227 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
228 *scanoCol = scanNo-1;
229 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
230 *cyclenoCol = cycleNo-1;
231 RecordFieldPtr<Double> mjdCol(rec, "TIME");
232 *mjdCol = mjd;
233 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
234 *intCol = interval;
235 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
236 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
237 // try to auto-identify if it is on or off.
238 Regex rx(".*_[e|w|R]$");
239 Regex rx2("_[e|w|R|S]$");
240 Int match = srcName.matches(rx);
241 if (match) {
242 *srcnCol = srcName;
243 } else {
244 *srcnCol = srcName.before(rx2);
245 }
246 //*srcnCol = srcName;//.before(rx2);
247 *srctCol = match;
248 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
249 *beamCol = beamNo-beamOffset_-1;
250 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
251 Int rb = -1;
252 if (nBeam_ > 1 ) rb = refBeam-1;
253 *rbCol = rb;
254 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
255 *ifCol = IFno-ifOffset_- 1;
256 uInt id;
257 /// @todo this has to change when nchan isn't global anymore
258 id = table_->frequencies().addEntry(Double(header_->nchan/2),
259 refFreq, freqInc);
260 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
261 *mfreqidCol = id;
262
263 id = table_->molecules().addEntry(restFreq);
264 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
265 *molidCol = id;
266
267 id = table_->tcal().addEntry(tcalTime, tcal);
268 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
269 *mcalidCol = id;
270 id = table_->weather().addEntry(temperature, pressure, humidity,
271 windSpeed, windAz);
272 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
273 *mweatheridCol = id;
274 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
275 id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
276 *mfocusidCol = id;
277 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
278 *dirCol = direction;
279 RecordFieldPtr<Double> azCol(rec, "AZIMUTH");
280 *azCol = azimuth;
281 RecordFieldPtr<Double> elCol(rec, "ELEVATION");
282 *elCol = elevation;
283
284 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
285 *parCol = parAngle;
286
287 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
288 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
289 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
290
291 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
292 // Turn the (nchan,npol) matrix and possible complex xPol vector
293 // into 2-4 rows in the scantable
294 Vector<Float> tsysvec(1);
295 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
296 uInt npol = (spectra.ncolumn()==1 ? 1: 2);
297 for ( uInt i=0; i< npol; ++i ) {
298 tsysvec = tsys(i);
299 *tsysCol = tsysvec;
300 *polnoCol = i;
301
302 *specCol = spectra.column(i);
303 *flagCol = flagtra.column(i);
304 table_->table().addRow();
305 row.put(table_->table().nrow()-1, rec);
306 }
307 if ( haveXPol_ ) {
308 // no tsys given for xpol, so emulate it
309 tsysvec = sqrt(tsys[0]*tsys[1]);
310 *tsysCol = tsysvec;
311 // add real part of cross pol
312 *polnoCol = 2;
313 Vector<Float> r(real(xPol));
314 *specCol = r;
315 // make up flags from linears
316 /// @fixme this has to be a bitwise or of both pols
317 *flagCol = flagtra.column(0);// | flagtra.column(1);
318 table_->table().addRow();
319 row.put(table_->table().nrow()-1, rec);
320 // ad imaginary part of cross pol
321 *polnoCol = 3;
322 Vector<Float> im(imag(xPol));
323 *specCol = im;
324 table_->table().addRow();
325 row.put(table_->table().nrow()-1, rec);
326 }
327 }
328 if (status > 0) {
329 close();
330 throw(AipsError("Reading error occured, data possibly corrupted."));
331 }
332
333 return status;
334}
335
336}//namespace asap
Note: See TracBrowser for help on using the repository browser.