source: trunk/src/STFiller.cpp@ 1731

Last change on this file since 1731 was 1725, checked in by Malte Marquarding, 15 years ago

Finishing touches to opacity calculations, docs, plotting and model

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