source: trunk/src/PKSFiller.cpp@ 1900

Last change on this file since 1900 was 1880, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

An implementation of asapmath.splitant() is changed since
the new filler doesn't support antenna parameter at the
moment.

The default behavior of MS filler also changed.


File size: 9.8 KB
Line 
1//
2// C++ Implementation: PKSFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2010
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#include <casa/Logging/LogIO.h>
23
24#include <measures/Measures/MDirection.h>
25#include <measures/Measures/MeasConvert.h>
26
27#include <atnf/PKSIO/PKSrecord.h>
28#include <atnf/PKSIO/PKSreader.h>
29
30#ifdef HAS_ALMA
31 #include <casa/System/ProgressMeter.h>
32#endif
33#include <casa/Logging/LogIO.h>
34
35#include <time.h>
36
37#include "STDefs.h"
38#include "STAttr.h"
39
40#include "PKSFiller.h"
41#include "STHeader.h"
42
43using namespace casa;
44
45namespace asap {
46
47
48PKSFiller::PKSFiller( CountedPtr< Scantable > stbl ) :
49 FillerBase(stbl),
50 reader_(0)
51{
52 setReferenceRegex(".*(e|w|_R)$");
53}
54
55PKSFiller::~PKSFiller()
56{
57 close();
58}
59
60bool PKSFiller::open( const std::string& filename)
61{
62 Bool haveBase, haveSpectra;
63
64 String inName(filename);
65 Path path(inName);
66 inName = path.expandedName();
67
68 File file(inName);
69 filename_ = inName;
70
71 // Create reader and fill in values for arguments
72 String format;
73 Vector<Bool> beams, ifs;
74 Vector<uInt> nchans,npols;
75
76 String antenna("");
77
78 reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
79 nchans, npols, haveXPol_, haveBase, haveSpectra);
80 if (reader_.null() == True) {
81 return False;
82 }
83
84 if (!haveSpectra) {
85 reader_ = 0;
86 throw(AipsError("No spectral data in file."));
87 }
88
89 LogIO os( LogOrigin( "PKSFiller", "open()", WHERE ) ) ;
90 nBeam_ = beams.nelements();
91 nIF_ = ifs.nelements();
92 // Get basic parameters.
93 if ( anyEQ(haveXPol_, True) ) {
94 os << "Cross polarization present" << LogIO::POST;
95 for (uInt i=0; i< npols.nelements();++i) {
96 if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
97 }
98 }
99 STHeader header;
100 header.nchan = max(nchans);
101 header.npol = max(npols);
102 header.nbeam = nBeam_;
103
104 Int status = reader_->getHeader(header.observer, header.project,
105 header.antennaname, header.antennaposition,
106 header.obstype,
107 header.fluxunit,
108 header.equinox,
109 header.freqref,
110 header.utc, header.reffreq,
111 header.bandwidth);
112
113 if (status) {
114 reader_ = 0;
115 throw(AipsError("Failed to get header."));
116 }
117 os << "Found " << header.antennaname << " data." << LogIO::POST;
118 if ((header.obstype).matches("*SW*")) {
119 // need robust way here - probably read ahead of next timestamp
120 os << "Header indicates frequency switched observation.\n"
121 << "setting # of IFs = 1 " << LogIO::POST;
122
123 nIF_ = 1;
124 header.obstype = String("fswitch");
125 }
126 // Determine Telescope and set brightness unit
127 Bool throwIt = False;
128 Instrument inst = STAttr::convertInstrument(header.antennaname, throwIt);
129
130 if (inst==ATMOPRA || inst==TIDBINBILLA) {
131 header.fluxunit = "K";
132 } else {
133 // downcase for use with Quanta
134 if (header.fluxunit == "JY") {
135 header.fluxunit = "Jy";
136 }
137 }
138 STAttr stattr;
139 header.poltype = stattr.feedPolType(inst);
140 header.nif = nIF_;
141 header.epoch = "UTC";
142 // *** header.frequnit = "Hz"
143 // Apply selection criteria.
144 Vector<Int> ref;
145
146 Vector<Int> start(nIF_, 1);
147 Vector<Int> end(nIF_, 0);
148 Bool getPt = False;
149 reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
150 setHeader(header);
151 //For MS, add the location of POINTING of the input MS so one get
152 //pointing data from there, if necessary.
153 //Also find nrow in MS
154 nInDataRow = 0;
155 if (format == "MS2") {
156 Path datapath(inName);
157 String ptTabPath = datapath.absoluteName();
158 Table inMS(ptTabPath);
159 nInDataRow = inMS.nrow();
160 ptTabPath.append("/POINTING");
161 table_->table().rwKeywordSet().define("POINTING", ptTabPath);
162 if (header.antennaname.matches("GBT")) {
163 String GOTabPath = datapath.absoluteName();
164 GOTabPath.append("/GBT_GO");
165 table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
166 }
167 }
168 String freqFrame = header.freqref;
169 //translate frequency reference frame back to
170 //MS style (as PKSMS2reader converts the original frame
171 //in FITS standard style)
172 if (freqFrame == "TOPOCENT") {
173 freqFrame = "TOPO";
174 } else if (freqFrame == "GEOCENER") {
175 freqFrame = "GEO";
176 } else if (freqFrame == "BARYCENT") {
177 freqFrame = "BARY";
178 } else if (freqFrame == "GALACTOC") {
179 freqFrame = "GALACTO";
180 } else if (freqFrame == "LOCALGRP") {
181 freqFrame = "LGROUP";
182 } else if (freqFrame == "CMBDIPOL") {
183 freqFrame = "CMB";
184 } else if (freqFrame == "SOURCE") {
185 freqFrame = "REST";
186 }
187 // set both "FRAME" and "BASEFRAME"
188 table_->frequencies().setFrame(freqFrame, false);
189 table_->frequencies().setFrame(freqFrame, true);
190
191 return true;
192}
193
194void PKSFiller::close( )
195{
196 if (reader_.null() != False) {
197 reader_ = 0;
198 }
199 table_ = 0;
200}
201
202void PKSFiller::fill( )
203{
204 int status = 0;
205 LogIO os( LogOrigin( "PKSFiller", "fill()", WHERE ) ) ;
206
207/**
208 Int beamNo, IFno, refBeam, scanNo, cycleNo;
209 Float azimuth, elevation, focusAxi, focusRot, focusTan,
210 humidity, parAngle, pressure, temperature, windAz, windSpeed;
211 Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
212 String fieldName, srcName, tcalTime, obsType;
213 Vector<Float> calFctr, sigma, tcal, tsys;
214 Matrix<Float> baseLin, baseSub;
215 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
216 Matrix<Float> spectra;
217 Matrix<uChar> flagtra;
218 Complex xCalFctr;
219 Vector<Complex> xPol;
220**/
221
222 Double min = 0.0;
223 Double max = nInDataRow;
224#ifdef HAS_ALMA
225 ProgressMeter fillpm(min, max, "Data importing progress");
226#endif
227 PKSrecord pksrec;
228 pksrec.srcType=-1;
229 int n = 0;
230 bool isGBTFITS = false ;
231 if ((table_->getAntennaName().find( "GBT" ) != String::npos)
232 && File(filename_).isRegular()) {
233 FILE *fp = fopen( filename_.c_str(), "r" ) ;
234 fseek( fp, 640, SEEK_SET ) ;
235 char buf[81] ;
236 fread( buf, 80, 1, fp ) ;
237 buf[80] = '\0' ;
238 if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
239 isGBTFITS = true ;
240 }
241 fclose( fp ) ;
242 }
243
244 while ( status == 0 ) {
245
246 status = reader_->read(pksrec);
247 if ( status != 0 ) break;
248 n += 1;
249 Regex filterrx(".*[SL|PA]$");
250 Regex obsrx("^AT.+");
251 if ( table_->getAntennaName().matches(obsrx) &&
252 pksrec.obsType.matches(filterrx)) {
253 //cerr << "ignoring paddle scan" << endl;
254 continue;
255 }
256
257 Vector<Double> sratedbl(pksrec.scanRate.nelements());
258 convertArray(sratedbl, pksrec.scanRate);
259 setScanRate(sratedbl);
260
261 Regex rx(getReferenceRegex());
262 Regex rx2("_S$");
263 Int match = pksrec.srcName.matches(rx);
264 std::string srcname;
265 Int srctype = Int(SrcType::NOTYPE);
266 if (match) {
267 srcname = pksrec.srcName;
268 srctype = Int(SrcType::PSOFF);
269 } else {
270 srcname = pksrec.srcName.before(rx2);
271 srctype = Int(SrcType::PSON);
272 }
273 if ( pksrec.srcType != Int(SrcType::NOTYPE)) {
274 srctype = pksrec.srcType ;
275 }
276 setTime(pksrec.mjd, pksrec.interval);
277 setSource(srcname, srctype, pksrec.fieldName,
278 pksrec.srcDir, pksrec.srcPM, pksrec.srcVel);
279
280 if (nBeam_ > 1 ) {
281 setReferenceBeam(pksrec.refBeam-1);
282 }
283
284 setMolecule(pksrec.restFreq);
285 setTcal(pksrec.tcalTime, pksrec.tcal);
286 setWeather(pksrec.temperature, pksrec.pressure,
287 pksrec.humidity, pksrec.windSpeed,
288 pksrec.windAz);
289 setFocus(pksrec.parAngle, pksrec.focusAxi, pksrec.focusTan,
290 pksrec.focusRot);
291 setDirection(pksrec.direction, pksrec.azimuth, pksrec.elevation);
292
293 setFrequency(Double(pksrec.spectra.nrow()/2),
294 pksrec.refFreq, pksrec.freqInc);
295
296 // Note: this is only one value for all polarisations!
297 setFlagrow(pksrec.flagrow);
298 // Turn the (nchan,npol) matrix and possible complex xPol vector
299 // into 2-4 rows in the scantable
300 Vector<Float> tsysvec(1);
301 // Why is spectra.ncolumn() == 3 for haveXPol_ == True
302 uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
303 uInt polno =0;
304 for ( uInt i=0; i< npol; ++i ) {
305 tsysvec = pksrec.tsys(i);
306 if (isGBTFITS) {
307 polno = pksrec.polNo ;
308 } else {
309 polno = i;
310 }
311 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
312 pksrec.beamNo-1);
313 setSpectrum(pksrec.spectra.column(i), pksrec.flagged.column(i), tsysvec);
314 commitRow();
315 }
316 if ( haveXPol_[0] ) {
317 // no tsys given for xpol, so emulate it
318 tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
319 // add real part of cross pol
320 polno = 2;
321 Vector<Float> r(real(pksrec.xPol));
322 // make up flags from linears
323 /// @fixme this has to be a bitwise or of both pols
324 /// pksrec.flagged.column(0) | pksrec.flagged.column(1);
325
326 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
327 pksrec.beamNo-1);
328 setSpectrum(r, pksrec.flagged.column(0), tsysvec);
329 commitRow();
330
331 // ad imaginary part of cross pol
332 polno = 3;
333 Vector<Float> im(imag(pksrec.xPol));
334 setIndex(pksrec.scanNo-1, pksrec.cycleNo-1, pksrec.IFno-1, polno,
335 pksrec.beamNo-1);
336 setSpectrum(im, pksrec.flagged.column(0), tsysvec);
337 commitRow();
338 }
339#ifdef HAS_ALMA
340 fillpm._update(n);
341#endif
342 }
343 if (status > 0) {
344 close();
345 throw(AipsError("Reading error occured, data possibly corrupted."));
346 }
347#ifdef HAS_ALMA
348 fillpm.done();
349#endif
350}
351
352}//namespace asap
Note: See TracBrowser for help on using the repository browser.