source: trunk/src/STFiller.cpp@ 1171

Last change on this file since 1171 was 1142, checked in by mar637, 18 years ago

undid last check-in as pksreader can handle this data now.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.3 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 beamOffset_ = 0;
170 if (whichBeam>=0) {
171 if (whichBeam>=0 && whichBeam<nBeam_) {
172 beams = False;
173 beams(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(beams, ifs, start, end, ref, True, haveXPol_[0], False);
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, obsType;
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, obsType, IFno,
217 refFreq, 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 Regex filterrx(".*[SL|PA]$");
226 Regex obsrx("^AT.+");
227 if ( header_->antennaname.matches(obsrx) && 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.