source: trunk/src/STFiller.cpp@ 1094

Last change on this file since 1094 was 1081, checked in by mar637, 18 years ago

using obsType to skip paddle scans

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