source: branches/Release2.0/src/STFiller.cpp@ 1153

Last change on this file since 1153 was 999, checked in by mar637, 19 years ago

added "redundant" fields from reader. The aren't used in asap but should be passed on for consistency.

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