source: trunk/src/STFiller.cpp@ 995

Last change on this file since 995 was 972, checked in by mar637, 19 years ago

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