source: trunk/src/STFiller.cpp@ 1068

Last change on this file since 1068 was 1060, checked in by mar637, 18 years ago

Tracking changes in PKSreader/writer for MOPS support.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 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;
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, IFno, refFreq,
218 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 TableRow row(table_->table());
227 TableRecord& rec = row.record();
228 // fields that don't get used and are just passed through asap
229 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
230 *srateCol = scanRate;
231 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
232 *spmCol = srcPM;
233 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
234 *sdirCol = srcDir;
235 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
236 *svelCol = srcVel;
237 // the real stuff
238 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
239 *fitCol = -1;
240 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
241 *scanoCol = scanNo-1;
242 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
243 *cyclenoCol = cycleNo-1;
244 RecordFieldPtr<Double> mjdCol(rec, "TIME");
245 *mjdCol = mjd;
246 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
247 *intCol = interval;
248 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
249 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
250 // try to auto-identify if it is on or off.
251 Regex rx(".*[e|w|_R]$");
252 Regex rx2("_S$");
253 Int match = srcName.matches(rx);
254 if (match) {
255 *srcnCol = srcName;
256 } else {
257 *srcnCol = srcName.before(rx2);
258 }
259 //*srcnCol = srcName;//.before(rx2);
260 *srctCol = match;
261 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
262 *beamCol = beamNo-beamOffset_-1;
263 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
264 Int rb = -1;
265 if (nBeam_ > 1 ) rb = refBeam-1;
266 *rbCol = rb;
267 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
268 *ifCol = IFno-ifOffset_- 1;
269 uInt id;
270 /// @todo this has to change when nchan isn't global anymore
271 id = table_->frequencies().addEntry(Double(header_->nchan/2),
272 refFreq, freqInc);
273 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
274 *mfreqidCol = id;
275
276 id = table_->molecules().addEntry(restFreq);
277 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
278 *molidCol = id;
279
280 id = table_->tcal().addEntry(tcalTime, tcal);
281 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
282 *mcalidCol = id;
283 id = table_->weather().addEntry(temperature, pressure, humidity,
284 windSpeed, windAz);
285 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
286 *mweatheridCol = id;
287 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
288 id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
289 *mfocusidCol = id;
290 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
291 *dirCol = direction;
292 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
293 *azCol = azimuth;
294 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
295 *elCol = elevation;
296
297 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
298 *parCol = parAngle;
299
300 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
301 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
302 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
303
304 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
305 // Turn the (nchan,npol) matrix and possible complex xPol vector
306 // into 2-4 rows in the scantable
307 Vector<Float> tsysvec(1);
308 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
309 uInt npol = (spectra.ncolumn()==1 ? 1: 2);
310 for ( uInt i=0; i< npol; ++i ) {
311 tsysvec = tsys(i);
312 *tsysCol = tsysvec;
313 *polnoCol = i;
314
315 *specCol = spectra.column(i);
316 *flagCol = flagtra.column(i);
317 table_->table().addRow();
318 row.put(table_->table().nrow()-1, rec);
319 }
320 if ( haveXPol_[0] ) {
321 // no tsys given for xpol, so emulate it
322 tsysvec = sqrt(tsys[0]*tsys[1]);
323 *tsysCol = tsysvec;
324 // add real part of cross pol
325 *polnoCol = 2;
326 Vector<Float> r(real(xPol));
327 *specCol = r;
328 // make up flags from linears
329 /// @fixme this has to be a bitwise or of both pols
330 *flagCol = flagtra.column(0);// | flagtra.column(1);
331 table_->table().addRow();
332 row.put(table_->table().nrow()-1, rec);
333 // ad imaginary part of cross pol
334 *polnoCol = 3;
335 Vector<Float> im(imag(xPol));
336 *specCol = im;
337 table_->table().addRow();
338 row.put(table_->table().nrow()-1, rec);
339 }
340 }
341 if (status > 0) {
342 close();
343 throw(AipsError("Reading error occured, data possibly corrupted."));
344 }
345 return status;
346}
347
348}//namespace asap
Note: See TracBrowser for help on using the repository browser.