source: branches/Release2.1/src/STFiller.cpp@ 2524

Last change on this file since 2524 was 1188, checked in by mar637, 18 years ago

changed feedtype to be string and added detection of feedtype to filler.

  • 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 Int status = reader_->getHeader(header_->observer, header_->project,
115 header_->antennaname, header_->antennaposition,
116 header_->obstype,header_->equinox,
117 header_->freqref,
118 header_->utc, header_->reffreq,
119 header_->bandwidth);
120
121 if (status) {
122 delete reader_;
123 reader_ = 0;
124 delete header_;
125 header_ = 0;
126 throw(AipsError("Failed to get header."));
127 }
128 if ((header_->obstype).matches("*SW*")) {
129 // need robust way here - probably read ahead of next timestamp
130 pushLog("Header indicates frequency switched observation.\n"
131 "setting # of IFs = 1 ");
132 nIF_ = 1;
133 header_->obstype = String("fswitch");
134 }
135 // Determine Telescope and set brightness unit
136
137 Bool throwIt = False;
138 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
139 header_->fluxunit = "Jy";
140 if (inst==ATMOPRA || inst==TIDBINBILLA) {
141 header_->fluxunit = "K";
142 }
143 STAttr stattr;
144 header_->poltype = stattr.feedPolType(inst);
145 header_->nif = nIF_;
146 header_->epoch = "UTC";
147 // *** header_->frequnit = "Hz"
148 // Apply selection criteria.
149 Vector<Int> ref;
150 ifOffset_ = 0;
151 if (whichIF>=0) {
152 if (whichIF>=0 && whichIF<nIF_) {
153 ifs = False;
154 ifs(whichIF) = True;
155 header_->nif = 1;
156 nIF_ = 1;
157 ifOffset_ = whichIF;
158 } else {
159 delete reader_;
160 reader_ = 0;
161 delete header_;
162 header_ = 0;
163 throw(AipsError("Illegal IF selection"));
164 }
165 }
166 beamOffset_ = 0;
167 if (whichBeam>=0) {
168 if (whichBeam>=0 && whichBeam<nBeam_) {
169 beams = False;
170 beams(whichBeam) = True;
171 header_->nbeam = 1;
172 nBeam_ = 1;
173 beamOffset_ = whichBeam;
174 } else {
175 delete reader_;
176 reader_ = 0;
177 delete header_;
178 header_ = 0;
179 throw(AipsError("Illegal Beam selection"));
180 }
181 }
182 Vector<Int> start(nIF_, 1);
183 Vector<Int> end(nIF_, 0);
184 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
185 table_->setHeader(*header_);
186}
187
188void STFiller::close( )
189{
190 delete reader_;reader_=0;
191 delete header_;header_=0;
192 table_ = 0;
193}
194
195int asap::STFiller::read( )
196{
197 int status = 0;
198
199 Int beamNo, IFno, refBeam, scanNo, cycleNo;
200 Float azimuth, elevation, focusAxi, focusRot, focusTan,
201 humidity, parAngle, pressure, temperature, windAz, windSpeed;
202 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
203 String fieldName, srcName, tcalTime, obsType;
204 Vector<Float> calFctr, sigma, tcal, tsys;
205 Matrix<Float> baseLin, baseSub;
206 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
207 Matrix<Float> spectra;
208 Matrix<uChar> flagtra;
209 Complex xCalFctr;
210 Vector<Complex> xPol;
211 while ( status == 0 ) {
212 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
213 srcName, srcDir, srcPM, srcVel, obsType, IFno,
214 refFreq, bandwidth, freqInc, restFreq, tcal, tcalTime,
215 azimuth, elevation, parAngle, focusAxi,
216 focusTan, focusRot, temperature, pressure,
217 humidity, windSpeed, windAz, refBeam,
218 beamNo, direction, scanRate,
219 tsys, sigma, calFctr, baseLin, baseSub,
220 spectra, flagtra, xCalFctr, xPol);
221 if ( status != 0 ) break;
222 Regex filterrx(".*[SL|PA]$");
223 Regex obsrx("^AT.+");
224 if ( header_->antennaname.matches(obsrx) && obsType.matches(filterrx)) {
225 //cerr << "ignoring paddle scan" << endl;
226 continue;
227 }
228 TableRow row(table_->table());
229 TableRecord& rec = row.record();
230 // fields that don't get used and are just passed through asap
231 RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
232 *srateCol = scanRate;
233 RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
234 *spmCol = srcPM;
235 RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
236 *sdirCol = srcDir;
237 RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
238 *svelCol = srcVel;
239 // the real stuff
240 RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
241 *fitCol = -1;
242 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
243 *scanoCol = scanNo-1;
244 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
245 *cyclenoCol = cycleNo-1;
246 RecordFieldPtr<Double> mjdCol(rec, "TIME");
247 *mjdCol = mjd;
248 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
249 *intCol = interval;
250 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
251 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
252 // try to auto-identify if it is on or off.
253 Regex rx(".*[e|w|_R]$");
254 Regex rx2("_S$");
255 Int match = srcName.matches(rx);
256 if (match) {
257 *srcnCol = srcName;
258 } else {
259 *srcnCol = srcName.before(rx2);
260 }
261 //*srcnCol = srcName;//.before(rx2);
262 *srctCol = match;
263 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
264 *beamCol = beamNo-beamOffset_-1;
265 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
266 Int rb = -1;
267 if (nBeam_ > 1 ) rb = refBeam-1;
268 *rbCol = rb;
269 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
270 *ifCol = IFno-ifOffset_- 1;
271 uInt id;
272 /// @todo this has to change when nchan isn't global anymore
273 id = table_->frequencies().addEntry(Double(header_->nchan/2),
274 refFreq, freqInc);
275 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
276 *mfreqidCol = id;
277
278 id = table_->molecules().addEntry(restFreq);
279 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
280 *molidCol = id;
281
282 id = table_->tcal().addEntry(tcalTime, tcal);
283 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
284 *mcalidCol = id;
285 id = table_->weather().addEntry(temperature, pressure, humidity,
286 windSpeed, windAz);
287 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
288 *mweatheridCol = id;
289 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
290 id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
291 *mfocusidCol = id;
292 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
293 *dirCol = direction;
294 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
295 *azCol = azimuth;
296 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
297 *elCol = elevation;
298
299 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
300 *parCol = parAngle;
301
302 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
303 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
304 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
305
306 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
307 // Turn the (nchan,npol) matrix and possible complex xPol vector
308 // into 2-4 rows in the scantable
309 Vector<Float> tsysvec(1);
310 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
311 uInt npol = (spectra.ncolumn()==1 ? 1: 2);
312 for ( uInt i=0; i< npol; ++i ) {
313 tsysvec = tsys(i);
314 *tsysCol = tsysvec;
315 *polnoCol = i;
316
317 *specCol = spectra.column(i);
318 *flagCol = flagtra.column(i);
319 table_->table().addRow();
320 row.put(table_->table().nrow()-1, rec);
321 }
322 if ( haveXPol_[0] ) {
323 // no tsys given for xpol, so emulate it
324 tsysvec = sqrt(tsys[0]*tsys[1]);
325 *tsysCol = tsysvec;
326 // add real part of cross pol
327 *polnoCol = 2;
328 Vector<Float> r(real(xPol));
329 *specCol = r;
330 // make up flags from linears
331 /// @fixme this has to be a bitwise or of both pols
332 *flagCol = flagtra.column(0);// | flagtra.column(1);
333 table_->table().addRow();
334 row.put(table_->table().nrow()-1, rec);
335 // ad imaginary part of cross pol
336 *polnoCol = 3;
337 Vector<Float> im(imag(xPol));
338 *specCol = im;
339 table_->table().addRow();
340 row.put(table_->table().nrow()-1, rec);
341 }
342 }
343 if (status > 0) {
344 close();
345 throw(AipsError("Reading error occured, data possibly corrupted."));
346 }
347 return status;
348}
349
350}//namespace asap
Note: See TracBrowser for help on using the repository browser.