source: trunk/src/STFiller.cpp@ 1412

Last change on this file since 1412 was 1410, checked in by Malte Marquarding, 17 years ago

Mark C added brightness unit to PKSreader/writer

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.3 KB
Line 
1//
2// C++ Implementation: STFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006-2007
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//#include <casa/System/ProgressMeter.h>
29
30
31#include "STDefs.h"
32#include "STAttr.h"
33
34#include "STFiller.h"
35#include "STHeader.h"
36
37using namespace casa;
38
39namespace asap {
40
41STFiller::STFiller() :
42 reader_(0),
43 header_(0),
44 table_(0)
45{
46}
47
48STFiller::STFiller( CountedPtr< Scantable > stbl ) :
49 reader_(0),
50 header_(0),
51 table_(stbl)
52{
53}
54
55STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
56 reader_(0),
57 header_(0),
58 table_(0)
59{
60 open(filename, whichIF, whichBeam);
61}
62
63STFiller::~STFiller()
64{
65 close();
66}
67
68void STFiller::open( const std::string& filename, int whichIF, int whichBeam )
69{
70 if (table_.null()) {
71 table_ = new Scantable();
72 }
73 if (reader_) { delete reader_; reader_ = 0; }
74 Bool haveBase, haveSpectra;
75
76 String inName(filename);
77 Path path(inName);
78 inName = path.expandedName();
79
80 File file(inName);
81 if ( !file.exists() ) {
82 throw(AipsError("File does not exist"));
83 }
84 filename_ = inName;
85
86 // Create reader and fill in values for arguments
87 String format;
88 Vector<Bool> beams, ifs;
89 Vector<uInt> nchans,npols;
90 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, ifs,
91 nchans, npols, haveXPol_,haveBase, haveSpectra
92 )) == 0 ) {
93 throw(AipsError("Creation of PKSreader failed"));
94 }
95 if (!haveSpectra) {
96 delete reader_;
97 reader_ = 0;
98 throw(AipsError("No spectral data in file."));
99 return;
100 }
101 nBeam_ = beams.nelements();
102 nIF_ = ifs.nelements();
103 // Get basic parameters.
104 if ( anyEQ(haveXPol_, True) ) {
105 pushLog("Cross polarization present");
106 for (uInt i=0; i< npols.nelements();++i) {
107 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
108 }
109 }
110 if (header_) delete header_;
111 header_ = new STHeader();
112 header_->nchan = max(nchans);
113 header_->npol = max(npols);
114 header_->nbeam = nBeam_;
115
116 Int status = reader_->getHeader(header_->observer, header_->project,
117 header_->antennaname, header_->antennaposition,
118 header_->obstype,
119 header_->fluxunit,
120 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
142 Bool throwIt = False;
143 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
144
145 if (inst==ATMOPRA || inst==TIDBINBILLA) {
146 header_->fluxunit = "K";
147 } else {
148 // downcase for use with Quanta
149 if (header_->fluxunit == "JY") {
150 header_->fluxunit = "Jy";
151 }
152 }
153 STAttr stattr;
154 header_->poltype = stattr.feedPolType(inst);
155 header_->nif = nIF_;
156 header_->epoch = "UTC";
157 // *** header_->frequnit = "Hz"
158 // Apply selection criteria.
159 Vector<Int> ref;
160 ifOffset_ = 0;
161 if (whichIF>=0) {
162 if (whichIF>=0 && whichIF<nIF_) {
163 ifs = False;
164 ifs(whichIF) = True;
165 header_->nif = 1;
166 nIF_ = 1;
167 ifOffset_ = whichIF;
168 } else {
169 delete reader_;
170 reader_ = 0;
171 delete header_;
172 header_ = 0;
173 throw(AipsError("Illegal IF selection"));
174 }
175 }
176 beamOffset_ = 0;
177 if (whichBeam>=0) {
178 if (whichBeam>=0 && whichBeam<nBeam_) {
179 beams = False;
180 beams(whichBeam) = True;
181 header_->nbeam = 1;
182 nBeam_ = 1;
183 beamOffset_ = whichBeam;
184 } else {
185 delete reader_;
186 reader_ = 0;
187 delete header_;
188 header_ = 0;
189 throw(AipsError("Illegal Beam selection"));
190 }
191 }
192 Vector<Int> start(nIF_, 1);
193 Vector<Int> end(nIF_, 0);
194 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
195 table_->setHeader(*header_);
196 //For MS, add the location of POINTING of the input MS so one get
197 //pointing data from there, if necessary.
198 //Also find nrow in MS
199 nInDataRow = 0;
200 if (format == "MS2") {
201 Path datapath(inName);
202 String ptTabPath = datapath.absoluteName();
203 Table inMS(ptTabPath);
204 nInDataRow = inMS.nrow();
205 ptTabPath.append("/POINTING");
206 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
207 if ((header_->antennaname).matches("GBT")) {
208 String GOTabPath = datapath.absoluteName();
209 GOTabPath.append("/GBT_GO");
210 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
211 }
212 }
213
214}
215
216void STFiller::close( )
217{
218 delete reader_;reader_=0;
219 delete header_;header_=0;
220 table_ = 0;
221}
222
223int asap::STFiller::read( )
224{
225 int status = 0;
226
227 Int beamNo, IFno, refBeam, scanNo, cycleNo;
228 Float azimuth, elevation, focusAxi, focusRot, focusTan,
229 humidity, parAngle, pressure, temperature, windAz, windSpeed;
230 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
231 String fieldName, srcName, tcalTime, obsType;
232 Vector<Float> calFctr, sigma, tcal, tsys;
233 Matrix<Float> baseLin, baseSub;
234 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
235 Matrix<Float> spectra;
236 Matrix<uChar> flagtra;
237 Complex xCalFctr;
238 Vector<Complex> xPol;
239 Double min = 0.0;
240 Double max = nInDataRow;
241 //ProgressMeter fillpm(min, max, "Data importing progress");
242 int n = 0;
243 while ( status == 0 ) {
244 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
245 srcName, srcDir, srcPM, srcVel, obsType, IFno,
246 refFreq, bandwidth, freqInc, restFreq, tcal, tcalTime,
247 azimuth, elevation, parAngle, focusAxi,
248 focusTan, focusRot, temperature, pressure,
249 humidity, windSpeed, windAz, refBeam,
250 beamNo, direction, scanRate,
251 tsys, sigma, calFctr, baseLin, baseSub,
252 spectra, flagtra, xCalFctr, xPol);
253 if ( status != 0 ) break;
254 n += 1;
255
256 Regex filterrx(".*[SL|PA]$");
257 Regex obsrx("^AT.+");
258 if ( header_->antennaname.matches(obsrx) && obsType.matches(filterrx)) {
259 //cerr << "ignoring paddle scan" << endl;
260 continue;
261 }
262 TableRow row(table_->table());
263 TableRecord& rec = row.record();
264 // fields that don't get used and are just passed through asap
265 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
266 *srateCol = scanRate;
267 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
268 *spmCol = srcPM;
269 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
270 *sdirCol = srcDir;
271 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
272 *svelCol = srcVel;
273 // the real stuff
274 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
275 *fitCol = -1;
276 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
277 *scanoCol = scanNo-1;
278 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
279 *cyclenoCol = cycleNo-1;
280 RecordFieldPtr<Double> mjdCol(rec, "TIME");
281 *mjdCol = mjd;
282 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
283 *intCol = interval;
284 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
285 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
286 RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
287 *fieldnCol = fieldName;
288 // try to auto-identify if it is on or off.
289 Regex rx(".*[e|w|_R]$");
290 Regex rx2("_S$");
291 Int match = srcName.matches(rx);
292 if (match) {
293 *srcnCol = srcName;
294 } else {
295 *srcnCol = srcName.before(rx2);
296 }
297 //*srcnCol = srcName;//.before(rx2);
298 *srctCol = match;
299 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
300 *beamCol = beamNo-beamOffset_-1;
301 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
302 Int rb = -1;
303 if (nBeam_ > 1 ) rb = refBeam-1;
304 *rbCol = rb;
305 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
306 *ifCol = IFno-ifOffset_- 1;
307 uInt id;
308 /// @todo this has to change when nchan isn't global anymore
309 id = table_->frequencies().addEntry(Double(header_->nchan/2),
310 refFreq, freqInc);
311 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
312 *mfreqidCol = id;
313
314 id = table_->molecules().addEntry(restFreq);
315 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
316 *molidCol = id;
317
318 id = table_->tcal().addEntry(tcalTime, tcal);
319 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
320 *mcalidCol = id;
321 id = table_->weather().addEntry(temperature, pressure, humidity,
322 windSpeed, windAz);
323 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
324 *mweatheridCol = id;
325 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
326 id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
327 *mfocusidCol = id;
328 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
329 *dirCol = direction;
330 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
331 *azCol = azimuth;
332 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
333 *elCol = elevation;
334
335 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
336 *parCol = parAngle;
337
338 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
339 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
340 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
341
342 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
343 // Turn the (nchan,npol) matrix and possible complex xPol vector
344 // into 2-4 rows in the scantable
345 Vector<Float> tsysvec(1);
346 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
347 uInt npol = (spectra.ncolumn()==1 ? 1: 2);
348 for ( uInt i=0; i< npol; ++i ) {
349 tsysvec = tsys(i);
350 *tsysCol = tsysvec;
351 *polnoCol = i;
352
353 *specCol = spectra.column(i);
354 *flagCol = flagtra.column(i);
355 table_->table().addRow();
356 row.put(table_->table().nrow()-1, rec);
357 }
358 if ( haveXPol_[0] ) {
359 // no tsys given for xpol, so emulate it
360 tsysvec = sqrt(tsys[0]*tsys[1]);
361 *tsysCol = tsysvec;
362 // add real part of cross pol
363 *polnoCol = 2;
364 Vector<Float> r(real(xPol));
365 *specCol = r;
366 // make up flags from linears
367 /// @fixme this has to be a bitwise or of both pols
368 *flagCol = flagtra.column(0);// | flagtra.column(1);
369 table_->table().addRow();
370 row.put(table_->table().nrow()-1, rec);
371 // ad imaginary part of cross pol
372 *polnoCol = 3;
373 Vector<Float> im(imag(xPol));
374 *specCol = im;
375 table_->table().addRow();
376 row.put(table_->table().nrow()-1, rec);
377 }
378 //fillpm._update(n);
379 }
380 if (status > 0) {
381 close();
382 throw(AipsError("Reading error occured, data possibly corrupted."));
383 }
384 //fillpm.done();
385 return status;
386}
387
388}//namespace asap
Note: See TracBrowser for help on using the repository browser.