source: trunk/src/STFiller.cpp@ 904

Last change on this file since 904 was 901, checked in by mar637, 19 years ago

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