source: tags/asap2alpha/src/STFiller.cpp@ 1141

Last change on this file since 1141 was 942, checked in by mar637, 18 years ago

fixed _S stripping in sourcename; reverted auto-detection of npol==1 == stokes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 9.5 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/Utilities/Regex.h>
21
22#include <casa/Containers/RecordField.h>
23
24#include <tables/Tables/TableRow.h>
25
26#include <atnf/PKSIO/PKSreader.h>
27
28#include "STDefs.h"
29#include "STAttr.h"
30
31#include "STFiller.h"
32#include "STHeader.h"
33
34using namespace casa;
35
36namespace asap {
37
38STFiller::STFiller() :
39 reader_(0),
40 header_(0),
41 table_(0)
42{
43}
44
45STFiller::STFiller( CountedPtr< Scantable > stbl ) :
46 reader_(0),
47 header_(0),
48 table_(stbl)
49{
50}
51
52STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
53 reader_(0),
54 header_(0),
55 table_(0)
56{
57 open(filename, whichIF, whichBeam);
58}
59
60STFiller::~STFiller()
61{
62 close();
63}
64
65void STFiller::open( const std::string& filename, int whichIF, int whichBeam )
66{
67 if (table_.null()) {
68 table_ = new Scantable();
69 }
70 if (reader_) { delete reader_; reader_ = 0; }
71 Bool haveBase, haveSpectra;
72
73 String inName(filename);
74 Path path(inName);
75 inName = path.expandedName();
76
77 File file(inName);
78 if ( !file.exists() ) {
79 throw(AipsError("File does not exist"));
80 }
81 filename_ = inName;
82
83 // Create reader and fill in values for arguments
84 String format;
85 Vector<Bool> beams;
86 if ( (reader_ = getPKSreader(inName, 0, 0, format, beams, nIF_,
87 nChan_, nPol_, haveBase, haveSpectra,
88 haveXPol_)) == 0 ) {
89 throw(AipsError("Creation of PKSreader failed"));
90 }
91 if (!haveSpectra) {
92 delete reader_;
93 reader_ = 0;
94 throw(AipsError("No spectral data in file."));
95 return;
96 }
97
98 nBeam_ = beams.nelements();
99 // Get basic parameters.
100 if ( haveXPol_ ) {
101 pushLog("Cross polarization present");
102 nPol_ += 2; // Convert Complex -> 2 Floats
103 }
104 if (header_) delete header_;
105 header_ = new STHeader();
106 header_->nchan = nChan_;
107 header_->npol = nPol_;
108 header_->nbeam = nBeam_;
109
110 // not the right thing to do?!
111 //if ( nPol_ == 1 ) header_->poltype = "stokes";
112 //else header_->poltype = "linear";
113 header_->poltype = "linear";
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
136 // Determine Telescope and set brightness unit
137
138 Bool throwIt = False;
139 Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
140 header_->fluxunit = "Jy";
141 if (inst==ATMOPRA || inst==TIDBINBILLA) {
142 header_->fluxunit = "K";
143 }
144
145 header_->nif = nIF_;
146 header_->epoch = "UTC";
147 // *** header_->frequnit = "Hz"
148 // Apply selection criteria.
149
150 Vector<Int> ref;
151 Vector<Bool> beamSel(nBeam_,True);
152 Vector<Bool> IFsel(nIF_,True);
153
154 ifOffset_ = 0;
155 if (whichIF>=0) {
156 if (whichIF>=0 && whichIF<nIF_) {
157 IFsel = False;
158 IFsel(whichIF) = True;
159 header_->nif = 1;
160 nIF_ = 1;
161 ifOffset_ = whichIF;
162 } else {
163 delete reader_;
164 reader_ = 0;
165 delete header_;
166 header_ = 0;
167 throw(AipsError("Illegal IF selection"));
168 }
169 }
170
171 beamOffset_ = 0;
172 if (whichBeam>=0) {
173 if (whichBeam>=0 && whichBeam<nBeam_) {
174 beamSel = False;
175 beamSel(whichBeam) = True;
176 header_->nbeam = 1;
177 nBeam_ = 1;
178 beamOffset_ = whichBeam;
179 } else {
180 delete reader_;
181 reader_ = 0;
182 delete header_;
183 header_ = 0;
184 throw(AipsError("Illegal Beam selection"));
185 }
186 }
187 Vector<Int> start(nIF_, 1);
188 Vector<Int> end(nIF_, 0);
189 reader_->select(beamSel, IFsel, start, end, ref, True, haveXPol_);
190 table_->setHeader(*header_);
191}
192
193void STFiller::close( )
194{
195 delete reader_;reader_=0;
196 delete header_;header_=0;
197 table_ = 0;
198}
199
200int asap::STFiller::read( )
201{
202 int status = 0;
203
204 Int beamNo, IFno, refBeam, scanNo, cycleNo;
205 Float azimuth, elevation, focusAxi, focusRot, focusTan,
206 humidity, parAngle, pressure, temperature, windAz, windSpeed;
207 Double bandwidth, freqInc, interval, mjd, refFreq, restFreq, srcVel;
208 String fieldName, srcName, tcalTime;
209 Vector<Float> calFctr, sigma, tcal, tsys;
210 Matrix<Float> baseLin, baseSub;
211 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2);
212 Matrix<Float> spectra;
213 Matrix<uChar> flagtra;
214 Complex xCalFctr;
215 Vector<Complex> xPol;
216 while ( status == 0 ) {
217 status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
218 srcName, srcDir, srcPM, srcVel, IFno, refFreq,
219 bandwidth, freqInc, restFreq, tcal, tcalTime,
220 azimuth, elevation, parAngle, focusAxi,
221 focusTan, focusRot, temperature, pressure,
222 humidity, windSpeed, windAz, refBeam,
223 beamNo, direction, scanRate,
224 tsys, sigma, calFctr, baseLin, baseSub,
225 spectra, flagtra, xCalFctr, xPol);
226 if ( status != 0 ) break;
227 TableRow row(table_->table());
228 TableRecord& rec = row.record();
229 RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
230 *scanoCol = scanNo-1;
231 RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
232 *cyclenoCol = cycleNo-1;
233 RecordFieldPtr<Double> mjdCol(rec, "TIME");
234 *mjdCol = mjd;
235 RecordFieldPtr<Double> intCol(rec, "INTERVAL");
236 *intCol = interval;
237 RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
238 RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
239 // try to auto-identify if it is on or off.
240 Regex rx(".*[e|w|_R]$");
241 Regex rx2("_S$");
242 Int match = srcName.matches(rx);
243 if (match) {
244 *srcnCol = srcName;
245 } else {
246 *srcnCol = srcName.before(rx2);
247 }
248 //*srcnCol = srcName;//.before(rx2);
249 *srctCol = match;
250 RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
251 *beamCol = beamNo-beamOffset_-1;
252 RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
253 Int rb = -1;
254 if (nBeam_ > 1 ) rb = refBeam-1;
255 *rbCol = rb;
256 RecordFieldPtr<uInt> ifCol(rec, "IFNO");
257 *ifCol = IFno-ifOffset_- 1;
258 uInt id;
259 /// @todo this has to change when nchan isn't global anymore
260 id = table_->frequencies().addEntry(Double(header_->nchan/2),
261 refFreq, freqInc);
262 RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
263 *mfreqidCol = id;
264
265 id = table_->molecules().addEntry(restFreq);
266 RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
267 *molidCol = id;
268
269 id = table_->tcal().addEntry(tcalTime, tcal);
270 RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
271 *mcalidCol = id;
272 id = table_->weather().addEntry(temperature, pressure, humidity,
273 windSpeed, windAz);
274 RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
275 *mweatheridCol = id;
276 RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
277 id = table_->focus().addEntry(focusAxi, focusTan, focusRot);
278 *mfocusidCol = id;
279 RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
280 *dirCol = direction;
281 RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
282 *azCol = azimuth;
283 RecordFieldPtr<Float> elCol(rec, "ELEVATION");
284 *elCol = elevation;
285
286 RecordFieldPtr<Float> parCol(rec, "PARANGLE");
287 *parCol = parAngle;
288
289 RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
290 RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
291 RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
292
293 RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
294 // Turn the (nchan,npol) matrix and possible complex xPol vector
295 // into 2-4 rows in the scantable
296 Vector<Float> tsysvec(1);
297 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
298 uInt npol = (spectra.ncolumn()==1 ? 1: 2);
299 for ( uInt i=0; i< npol; ++i ) {
300 tsysvec = tsys(i);
301 *tsysCol = tsysvec;
302 *polnoCol = i;
303
304 *specCol = spectra.column(i);
305 *flagCol = flagtra.column(i);
306 table_->table().addRow();
307 row.put(table_->table().nrow()-1, rec);
308 }
309 if ( haveXPol_ ) {
310 // no tsys given for xpol, so emulate it
311 tsysvec = sqrt(tsys[0]*tsys[1]);
312 *tsysCol = tsysvec;
313 // add real part of cross pol
314 *polnoCol = 2;
315 Vector<Float> r(real(xPol));
316 *specCol = r;
317 // make up flags from linears
318 /// @fixme this has to be a bitwise or of both pols
319 *flagCol = flagtra.column(0);// | flagtra.column(1);
320 table_->table().addRow();
321 row.put(table_->table().nrow()-1, rec);
322 // ad imaginary part of cross pol
323 *polnoCol = 3;
324 Vector<Float> im(imag(xPol));
325 *specCol = im;
326 table_->table().addRow();
327 row.put(table_->table().nrow()-1, rec);
328 }
329 }
330 if (status > 0) {
331 close();
332 throw(AipsError("Reading error occured, data possibly corrupted."));
333 }
334 return status;
335}
336
337}//namespace asap
Note: See TracBrowser for help on using the repository browser.