1 |
|
---|
2 | //
|
---|
3 | // C++ Interface: MSFiller
|
---|
4 | //
|
---|
5 | // Description:
|
---|
6 | //
|
---|
7 | // This class is specific filler for MS format
|
---|
8 | // New version that is implemented using TableVisitor instead of TableIterator
|
---|
9 | //
|
---|
10 | // Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
|
---|
11 | //
|
---|
12 | // Copyright: See COPYING file that comes with this distribution
|
---|
13 | //
|
---|
14 | //
|
---|
15 |
|
---|
16 | #include <assert.h>
|
---|
17 | #include <iostream>
|
---|
18 | #include <map>
|
---|
19 |
|
---|
20 | #include <tables/Tables/ExprNode.h>
|
---|
21 | #include <tables/Tables/TableIter.h>
|
---|
22 | #include <tables/Tables/TableColumn.h>
|
---|
23 | #include <tables/Tables/ScalarColumn.h>
|
---|
24 | #include <tables/Tables/ArrayColumn.h>
|
---|
25 | #include <tables/Tables/TableParse.h>
|
---|
26 | #include <tables/Tables/TableRow.h>
|
---|
27 |
|
---|
28 | #include <casa/Containers/RecordField.h>
|
---|
29 | #include <casa/Logging/LogIO.h>
|
---|
30 | #include <casa/Arrays/Slicer.h>
|
---|
31 | #include <casa/Quanta/MVTime.h>
|
---|
32 | #include <casa/OS/Path.h>
|
---|
33 |
|
---|
34 | #include <measures/Measures/Stokes.h>
|
---|
35 | #include <measures/Measures/MEpoch.h>
|
---|
36 | #include <measures/Measures/MCEpoch.h>
|
---|
37 | #include <measures/Measures/MFrequency.h>
|
---|
38 | #include <measures/Measures/MCFrequency.h>
|
---|
39 | #include <measures/Measures/MPosition.h>
|
---|
40 | #include <measures/Measures/MCPosition.h>
|
---|
41 | #include <measures/Measures/MDirection.h>
|
---|
42 | #include <measures/Measures/MCDirection.h>
|
---|
43 | #include <measures/Measures/MeasConvert.h>
|
---|
44 | #include <measures/TableMeasures/ScalarMeasColumn.h>
|
---|
45 | #include <measures/TableMeasures/ArrayMeasColumn.h>
|
---|
46 | #include <measures/TableMeasures/ScalarQuantColumn.h>
|
---|
47 | #include <measures/TableMeasures/ArrayQuantColumn.h>
|
---|
48 |
|
---|
49 | #include <ms/MeasurementSets/MSAntennaIndex.h>
|
---|
50 |
|
---|
51 | #include <atnf/PKSIO/SrcType.h>
|
---|
52 |
|
---|
53 | #include "MSFiller.h"
|
---|
54 | #include "STHeader.h"
|
---|
55 |
|
---|
56 | #include "MathUtils.h"
|
---|
57 |
|
---|
58 | using namespace casa ;
|
---|
59 | using namespace std ;
|
---|
60 |
|
---|
61 | namespace asap {
|
---|
62 |
|
---|
63 | class BaseMSFillerVisitor: public TableVisitor {
|
---|
64 | uInt lastRecordNo ;
|
---|
65 | Int lastObservationId ;
|
---|
66 | Int lastFeedId ;
|
---|
67 | Int lastFieldId ;
|
---|
68 | Int lastDataDescId ;
|
---|
69 | Int lastScanNo ;
|
---|
70 | Int lastStateId ;
|
---|
71 | Double lastTime ;
|
---|
72 | protected:
|
---|
73 | const Table &table;
|
---|
74 | uInt count;
|
---|
75 | public:
|
---|
76 | BaseMSFillerVisitor(const Table &table)
|
---|
77 | : table(table)
|
---|
78 | {
|
---|
79 | count = 0;
|
---|
80 | }
|
---|
81 |
|
---|
82 | virtual void enterObservationId(const uInt recordNo, Int columnValue) { }
|
---|
83 | virtual void leaveObservationId(const uInt recordNo, Int columnValue) { }
|
---|
84 | virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
|
---|
85 | virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
|
---|
86 | virtual void enterFieldId(const uInt recordNo, Int columnValue) { }
|
---|
87 | virtual void leaveFieldId(const uInt recordNo, Int columnValue) { }
|
---|
88 | virtual void enterDataDescId(const uInt recordNo, Int columnValue) { }
|
---|
89 | virtual void leaveDataDescId(const uInt recordNo, Int columnValue) { }
|
---|
90 | virtual void enterScanNo(const uInt recordNo, Int columnValue) { }
|
---|
91 | virtual void leaveScanNo(const uInt recordNo, Int columnValue) { }
|
---|
92 | virtual void enterStateId(const uInt recordNo, Int columnValue) { }
|
---|
93 | virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
|
---|
94 | virtual void enterTime(const uInt recordNo, Double columnValue) { }
|
---|
95 | virtual void leaveTime(const uInt recordNo, Double columnValue) { }
|
---|
96 |
|
---|
97 | virtual Bool visitRecord(const uInt recordNo,
|
---|
98 | const Int ObservationId,
|
---|
99 | const Int feedId,
|
---|
100 | const Int fieldId,
|
---|
101 | const Int dataDescId,
|
---|
102 | const Int scanNo,
|
---|
103 | const Int stateId,
|
---|
104 | const Double time) { return True ; }
|
---|
105 |
|
---|
106 | virtual Bool visit(Bool isFirst, const uInt recordNo,
|
---|
107 | const uInt nCols, void const *const colValues[]) {
|
---|
108 | Int observationId, feedId, fieldId, dataDescId, scanNo, stateId;
|
---|
109 | Double time;
|
---|
110 | { // prologue
|
---|
111 | uInt i = 0;
|
---|
112 | {
|
---|
113 | const Int *col = (const Int *)colValues[i++];
|
---|
114 | observationId = col[recordNo];
|
---|
115 | }
|
---|
116 | {
|
---|
117 | const Int *col = (const Int *)colValues[i++];
|
---|
118 | feedId = col[recordNo];
|
---|
119 | }
|
---|
120 | {
|
---|
121 | const Int *col = (const Int *)colValues[i++];
|
---|
122 | fieldId = col[recordNo];
|
---|
123 | }
|
---|
124 | {
|
---|
125 | const Int *col = (const Int *)colValues[i++];
|
---|
126 | dataDescId = col[recordNo];
|
---|
127 | }
|
---|
128 | {
|
---|
129 | const Int *col = (const Int *)colValues[i++];
|
---|
130 | scanNo = col[recordNo];
|
---|
131 | }
|
---|
132 | {
|
---|
133 | const Int *col = (const Int *)colValues[i++];
|
---|
134 | stateId = col[recordNo];
|
---|
135 | }
|
---|
136 | {
|
---|
137 | const Double *col = (const Double *)colValues[i++];
|
---|
138 | time = col[recordNo];
|
---|
139 | }
|
---|
140 | assert(nCols == i);
|
---|
141 | }
|
---|
142 |
|
---|
143 | if (isFirst) {
|
---|
144 | enterObservationId(recordNo, observationId);
|
---|
145 | enterFeedId(recordNo, feedId);
|
---|
146 | enterFieldId(recordNo, fieldId);
|
---|
147 | enterDataDescId(recordNo, dataDescId);
|
---|
148 | enterScanNo(recordNo, scanNo);
|
---|
149 | enterStateId(recordNo, stateId);
|
---|
150 | enterTime(recordNo, time);
|
---|
151 | } else {
|
---|
152 | if (lastObservationId != observationId) {
|
---|
153 | leaveTime(lastRecordNo, lastTime);
|
---|
154 | leaveStateId(lastRecordNo, lastStateId);
|
---|
155 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
156 | leaveDataDescId(lastRecordNo, lastDataDescId);
|
---|
157 | leaveFieldId(lastRecordNo, lastFieldId);
|
---|
158 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
159 | leaveObservationId(lastRecordNo, lastObservationId);
|
---|
160 |
|
---|
161 | enterObservationId(recordNo, observationId);
|
---|
162 | enterFeedId(recordNo, feedId);
|
---|
163 | enterFieldId(recordNo, fieldId);
|
---|
164 | enterDataDescId(recordNo, dataDescId);
|
---|
165 | enterScanNo(recordNo, scanNo);
|
---|
166 | enterStateId(recordNo, stateId);
|
---|
167 | enterTime(recordNo, time);
|
---|
168 | } else if (lastFeedId != feedId) {
|
---|
169 | leaveTime(lastRecordNo, lastTime);
|
---|
170 | leaveStateId(lastRecordNo, lastStateId);
|
---|
171 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
172 | leaveDataDescId(lastRecordNo, lastDataDescId);
|
---|
173 | leaveFieldId(lastRecordNo, lastFieldId);
|
---|
174 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
175 |
|
---|
176 | enterFeedId(recordNo, feedId);
|
---|
177 | enterFieldId(recordNo, fieldId);
|
---|
178 | enterDataDescId(recordNo, dataDescId);
|
---|
179 | enterScanNo(recordNo, scanNo);
|
---|
180 | enterStateId(recordNo, stateId);
|
---|
181 | enterTime(recordNo, time);
|
---|
182 | } else if (lastFieldId != fieldId) {
|
---|
183 | leaveTime(lastRecordNo, lastTime);
|
---|
184 | leaveStateId(lastRecordNo, lastStateId);
|
---|
185 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
186 | leaveDataDescId(lastRecordNo, lastDataDescId);
|
---|
187 | leaveFieldId(lastRecordNo, lastFieldId);
|
---|
188 |
|
---|
189 | enterFieldId(recordNo, fieldId);
|
---|
190 | enterDataDescId(recordNo, dataDescId);
|
---|
191 | enterScanNo(recordNo, scanNo);
|
---|
192 | enterStateId(recordNo, stateId);
|
---|
193 | enterTime(recordNo, time);
|
---|
194 | } else if (lastDataDescId != dataDescId) {
|
---|
195 | leaveTime(lastRecordNo, lastTime);
|
---|
196 | leaveStateId(lastRecordNo, lastStateId);
|
---|
197 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
198 | leaveDataDescId(lastRecordNo, lastDataDescId);
|
---|
199 |
|
---|
200 | enterDataDescId(recordNo, dataDescId);
|
---|
201 | enterScanNo(recordNo, scanNo);
|
---|
202 | enterStateId(recordNo, stateId);
|
---|
203 | enterTime(recordNo, time);
|
---|
204 | } else if (lastScanNo != scanNo) {
|
---|
205 | leaveTime(lastRecordNo, lastTime);
|
---|
206 | leaveStateId(lastRecordNo, lastStateId);
|
---|
207 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
208 |
|
---|
209 | enterScanNo(recordNo, scanNo);
|
---|
210 | enterStateId(recordNo, stateId);
|
---|
211 | enterTime(recordNo, time);
|
---|
212 | } else if (lastStateId != stateId) {
|
---|
213 | leaveTime(lastRecordNo, lastTime);
|
---|
214 | leaveStateId(lastRecordNo, lastStateId);
|
---|
215 |
|
---|
216 | enterStateId(recordNo, stateId);
|
---|
217 | enterTime(recordNo, time);
|
---|
218 | } else if (lastTime != time) {
|
---|
219 | leaveTime(lastRecordNo, lastTime);
|
---|
220 | enterTime(recordNo, time);
|
---|
221 | }
|
---|
222 | }
|
---|
223 | count++;
|
---|
224 | Bool result = visitRecord(recordNo, observationId, feedId, fieldId, dataDescId,
|
---|
225 | scanNo, stateId, time);
|
---|
226 |
|
---|
227 | { // epilogue
|
---|
228 | lastRecordNo = recordNo;
|
---|
229 |
|
---|
230 | lastObservationId = observationId;
|
---|
231 | lastFeedId = feedId;
|
---|
232 | lastFieldId = fieldId;
|
---|
233 | lastDataDescId = dataDescId;
|
---|
234 | lastScanNo = scanNo;
|
---|
235 | lastStateId = stateId;
|
---|
236 | lastTime = time;
|
---|
237 | }
|
---|
238 | return result ;
|
---|
239 | }
|
---|
240 |
|
---|
241 | virtual void finish() {
|
---|
242 | if (count > 0) {
|
---|
243 | leaveTime(lastRecordNo, lastTime);
|
---|
244 | leaveStateId(lastRecordNo, lastStateId);
|
---|
245 | leaveScanNo(lastRecordNo, lastScanNo);
|
---|
246 | leaveDataDescId(lastRecordNo, lastDataDescId);
|
---|
247 | leaveFieldId(lastRecordNo, lastFieldId);
|
---|
248 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
249 | leaveObservationId(lastRecordNo, lastObservationId);
|
---|
250 | }
|
---|
251 | }
|
---|
252 | };
|
---|
253 |
|
---|
254 | class MSFillerVisitor: public BaseMSFillerVisitor, public MSFillerUtils {
|
---|
255 | public:
|
---|
256 | MSFillerVisitor(const Table &from, Scantable &to)
|
---|
257 | : BaseMSFillerVisitor(from),
|
---|
258 | scantable( to )
|
---|
259 | {
|
---|
260 | antennaId = 0 ;
|
---|
261 | rowidx = 0 ;
|
---|
262 | tablerow = TableRow( scantable.table() ) ;
|
---|
263 | feedEntry = Vector<Int>( 64, -1 ) ;
|
---|
264 | nbeam = 0 ;
|
---|
265 | ifmap.clear() ;
|
---|
266 | const TableDesc &desc = table.tableDesc() ;
|
---|
267 | if ( desc.isColumn( "DATA" ) )
|
---|
268 | dataColumnName = "DATA" ;
|
---|
269 | else if ( desc.isColumn( "FLOAT_DATA" ) )
|
---|
270 | dataColumnName = "FLOAT_DATA" ;
|
---|
271 | getpt = False ;
|
---|
272 | isWeather = False ;
|
---|
273 | isSysCal = False ;
|
---|
274 | isTcal = False ;
|
---|
275 | cycleNo = 0 ;
|
---|
276 | numSysCalRow = 0 ;
|
---|
277 | header = scantable.getHeader() ;
|
---|
278 | fluxUnit( header.fluxunit ) ;
|
---|
279 |
|
---|
280 | // MS subtables
|
---|
281 | const TableRecord &hdr = table.keywordSet();
|
---|
282 | obstab = hdr.asTable( "OBSERVATION" ) ;
|
---|
283 | spwtab = hdr.asTable( "SPECTRAL_WINDOW" ) ;
|
---|
284 | statetab = hdr.asTable( "STATE" ) ;
|
---|
285 | ddtab = hdr.asTable( "DATA_DESCRIPTION" ) ;
|
---|
286 | poltab = hdr.asTable( "POLARIZATION" ) ;
|
---|
287 | fieldtab = hdr.asTable( "FIELD" ) ;
|
---|
288 | anttab = hdr.asTable( "ANTENNA" ) ;
|
---|
289 | if ( hdr.isDefined( "SYSCAL" ) )
|
---|
290 | sctab = hdr.asTable( "SYSCAL" ) ;
|
---|
291 | if ( hdr.isDefined( "SOURCE" ) )
|
---|
292 | srctab = hdr.asTable( "SOURCE" ) ;
|
---|
293 |
|
---|
294 | // attach to columns
|
---|
295 | // MS MAIN
|
---|
296 | intervalCol.attach( table, "INTERVAL" ) ;
|
---|
297 | flagRowCol.attach( table, "FLAG_ROW" ) ;
|
---|
298 | flagCol.attach( table, "FLAG" ) ;
|
---|
299 | if ( dataColumnName.compare( "DATA" ) == 0 )
|
---|
300 | dataCol.attach( table, dataColumnName ) ;
|
---|
301 | else
|
---|
302 | floatDataCol.attach( table, dataColumnName ) ;
|
---|
303 |
|
---|
304 | // set dummy epoch
|
---|
305 | mf.set( currentTime ) ;
|
---|
306 |
|
---|
307 | // initialize dirType
|
---|
308 | dirType = MDirection::N_Types ;
|
---|
309 |
|
---|
310 | //
|
---|
311 | // add rows to scantable
|
---|
312 | //
|
---|
313 | // number of polarization is up to 4
|
---|
314 | uInt addrow = table.nrow() * maxNumPol() ;
|
---|
315 | scantable.table().addRow( addrow ) ;
|
---|
316 |
|
---|
317 | // attach to columns
|
---|
318 | // Scantable MAIN
|
---|
319 | TableRecord &r = tablerow.record() ;
|
---|
320 | timeRF.attachToRecord( r, "TIME" ) ;
|
---|
321 | intervalRF.attachToRecord( r, "INTERVAL" ) ;
|
---|
322 | directionRF.attachToRecord( r, "DIRECTION" ) ;
|
---|
323 | azimuthRF.attachToRecord( r, "AZIMUTH" ) ;
|
---|
324 | elevationRF.attachToRecord( r, "ELEVATION" ) ;
|
---|
325 | scanRateRF.attachToRecord( r, "SCANRATE" ) ;
|
---|
326 | weatherIdRF.attachToRecord( r, "WEATHER_ID" ) ;
|
---|
327 | cycleNoRF.attachToRecord( r, "CYCLENO" ) ;
|
---|
328 | flagRowRF.attachToRecord( r, "FLAGROW" ) ;
|
---|
329 | polNoRF.attachToRecord( r, "POLNO" ) ;
|
---|
330 | tcalIdRF.attachToRecord( r, "TCAL_ID" ) ;
|
---|
331 | spectraRF.attachToRecord( r, "SPECTRA" ) ;
|
---|
332 | flagtraRF.attachToRecord( r, "FLAGTRA" ) ;
|
---|
333 | tsysRF.attachToRecord( r, "TSYS" ) ;
|
---|
334 | beamNoRF.attachToRecord( r, "BEAMNO" ) ;
|
---|
335 | ifNoRF.attachToRecord( r, "IFNO" ) ;
|
---|
336 | freqIdRF.attachToRecord( r, "FREQ_ID" ) ;
|
---|
337 | moleculeIdRF.attachToRecord( r, "MOLECULE_ID" ) ;
|
---|
338 | sourceNameRF.attachToRecord( r, "SRCNAME" ) ;
|
---|
339 | sourceProperMotionRF.attachToRecord( r, "SRCPROPERMOTION" ) ;
|
---|
340 | sourceDirectionRF.attachToRecord( r, "SRCDIRECTION" ) ;
|
---|
341 | sourceVelocityRF.attachToRecord( r, "SRCVELOCITY" ) ;
|
---|
342 | focusIdRF.attachToRecord( r, "FOCUS_ID" ) ;
|
---|
343 | fieldNameRF.attachToRecord( r, "FIELDNAME" ) ;
|
---|
344 | sourceTypeRF.attachToRecord( r, "SRCTYPE" ) ;
|
---|
345 | scanNoRF.attachToRecord( r, "SCANNO" ) ;
|
---|
346 |
|
---|
347 | // put values
|
---|
348 | RecordFieldPtr<Int> refBeamNoRF( r, "REFBEAMNO" ) ;
|
---|
349 | *refBeamNoRF = -1 ;
|
---|
350 | RecordFieldPtr<Int> fitIdRF( r, "FIT_ID" ) ;
|
---|
351 | *fitIdRF = -1 ;
|
---|
352 | RecordFieldPtr<Float> opacityRF( r, "OPACITY" ) ;
|
---|
353 | *opacityRF = 0.0 ;
|
---|
354 | }
|
---|
355 |
|
---|
356 | virtual void enterObservationId(const uInt recordNo, Int columnValue) {
|
---|
357 | //printf("%u: ObservationId: %d\n", recordNo, columnValue);
|
---|
358 | // update header
|
---|
359 | if ( header.observer.empty() )
|
---|
360 | getScalar( String("OBSERVER"), (uInt)columnValue, obstab, header.observer ) ;
|
---|
361 | if ( header.project.empty() )
|
---|
362 | getScalar( "PROJECT", (uInt)columnValue, obstab, header.project ) ;
|
---|
363 | if ( header.utc == 0.0 ) {
|
---|
364 | Vector<MEpoch> amp ;
|
---|
365 | getArrayMeas( "TIME_RANGE", (uInt)columnValue, obstab, amp ) ;
|
---|
366 | header.utc = amp[0].get( "d" ).getValue() ;
|
---|
367 | }
|
---|
368 | if ( header.antennaname.empty() )
|
---|
369 | getScalar( "TELESCOPE_NAME", (uInt)columnValue, obstab, header.antennaname ) ;
|
---|
370 | }
|
---|
371 | virtual void leaveObservationId(const uInt recordNo, Int columnValue) {
|
---|
372 | // update header
|
---|
373 | header.nbeam = max( header.nbeam, (Int)nbeam ) ;
|
---|
374 |
|
---|
375 | nbeam = 0 ;
|
---|
376 | feedEntry = -1 ;
|
---|
377 | }
|
---|
378 | virtual void enterFeedId(const uInt recordNo, Int columnValue) {
|
---|
379 | //printf("%u: FeedId: %d\n", recordNo, columnValue);
|
---|
380 |
|
---|
381 | // update feed entry
|
---|
382 | if ( allNE( feedEntry, columnValue ) ) {
|
---|
383 | feedEntry[nbeam] = columnValue ;
|
---|
384 | nbeam++ ;
|
---|
385 | }
|
---|
386 |
|
---|
387 | // put values
|
---|
388 | *beamNoRF = (uInt)columnValue ;
|
---|
389 | *focusIdRF = (uInt)0 ;
|
---|
390 | }
|
---|
391 | virtual void leaveFeedId(const uInt recordNo, Int columnValue) {
|
---|
392 | Int nelem = (Int)feedEntry.nelements() ;
|
---|
393 | if ( nbeam > nelem ) {
|
---|
394 | feedEntry.resize( nelem+64, True ) ;
|
---|
395 | Slicer slice( IPosition( 1, nelem ), IPosition( 1, feedEntry.nelements()-1 ) ) ;
|
---|
396 | feedEntry( slice ) = -1 ;
|
---|
397 | }
|
---|
398 | }
|
---|
399 | virtual void enterFieldId(const uInt recordNo, Int columnValue) {
|
---|
400 | //printf("%u: FieldId: %d\n", recordNo, columnValue);
|
---|
401 | // update sourceId and fieldName
|
---|
402 | getScalar( "SOURCE_ID", (uInt)columnValue, fieldtab, sourceId ) ;
|
---|
403 | String fieldName ;
|
---|
404 | getScalar( "NAME", (uInt)columnValue, fieldtab, fieldName ) ;
|
---|
405 | fieldName += "__" + String::toString( columnValue ) ;
|
---|
406 |
|
---|
407 | // put values
|
---|
408 | *fieldNameRF = fieldName ;
|
---|
409 | }
|
---|
410 | virtual void leaveFieldId(const uInt recordNo, Int columnValue) {
|
---|
411 | sourceId = -1 ;
|
---|
412 | }
|
---|
413 | virtual void enterDataDescId(const uInt recordNo, Int columnValue) {
|
---|
414 | //printf("%u: DataDescId: %d\n", recordNo, columnValue);
|
---|
415 | // update polarization and spectral window ids
|
---|
416 | getScalar( "POLARIZATION_ID", (uInt)columnValue, ddtab, polId ) ;
|
---|
417 | getScalar( "SPECTRAL_WINDOW_ID", (uInt)columnValue, ddtab, spwId ) ;
|
---|
418 |
|
---|
419 | // polarization setup
|
---|
420 | getScalar( "NUM_CORR", (uInt)polId, poltab, npol ) ;
|
---|
421 | Vector<Int> corrtype ;
|
---|
422 | getArray( "CORR_TYPE", (uInt)polId, poltab, corrtype ) ;
|
---|
423 | polnos = getPolNos( corrtype ) ;
|
---|
424 |
|
---|
425 | // process SOURCE table
|
---|
426 | String sourceName ;
|
---|
427 | Vector<Double> sourcePM, restFreqs, sysVels ;
|
---|
428 | Vector<String> transition ;
|
---|
429 | processSource( sourceId, spwId, sourceName, sourceDir, sourcePM,
|
---|
430 | restFreqs, transition, sysVels ) ;
|
---|
431 |
|
---|
432 | // spectral setup
|
---|
433 | uInt freqId ;
|
---|
434 | Double reffreq, bandwidth ;
|
---|
435 | String freqref ;
|
---|
436 | getScalar( "NUM_CHAN", (uInt)spwId, spwtab, nchan ) ;
|
---|
437 | Bool iswvr = (Bool)(nchan == 4) ;
|
---|
438 | map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
|
---|
439 | if ( iter == ifmap.end() ) {
|
---|
440 | MEpoch me ;
|
---|
441 | getScalarMeas( "TIME", recordNo, table, me ) ;
|
---|
442 | spectralSetup( spwId, me, antpos, sourceDir,
|
---|
443 | freqId, nchan,
|
---|
444 | freqref, reffreq, bandwidth ) ;
|
---|
445 | ifmap.insert( pair<Int,uInt>(spwId,freqId) ) ;
|
---|
446 | }
|
---|
447 | else {
|
---|
448 | freqId = iter->second ;
|
---|
449 | }
|
---|
450 | sp.resize( npol, nchan ) ;
|
---|
451 | fl.resize( npol, nchan ) ;
|
---|
452 |
|
---|
453 |
|
---|
454 | // molecular setup
|
---|
455 | STMolecules mtab = scantable.molecules() ;
|
---|
456 | uInt molId = mtab.addEntry( restFreqs, transition, transition ) ;
|
---|
457 |
|
---|
458 | // process SYSCAL table
|
---|
459 | if ( isSysCal )
|
---|
460 | processSysCal( spwId ) ;
|
---|
461 |
|
---|
462 | // update header
|
---|
463 | if ( !iswvr ) {
|
---|
464 | header.nchan = max( header.nchan, nchan ) ;
|
---|
465 | header.bandwidth = max( header.bandwidth, bandwidth ) ;
|
---|
466 | if ( header.reffreq == -1.0 )
|
---|
467 | header.reffreq = reffreq ;
|
---|
468 | header.npol = max( header.npol, npol ) ;
|
---|
469 | if ( header.poltype.empty() )
|
---|
470 | header.poltype = getPolType( corrtype[0] ) ;
|
---|
471 | if ( header.freqref.empty() )
|
---|
472 | header.freqref = freqref ;
|
---|
473 | }
|
---|
474 |
|
---|
475 | // put values
|
---|
476 | *ifNoRF = (uInt)spwId ;
|
---|
477 | *freqIdRF = freqId ;
|
---|
478 | *moleculeIdRF = molId ;
|
---|
479 | *sourceNameRF = sourceName ;
|
---|
480 | sourceProperMotionRF.define( sourcePM ) ;
|
---|
481 | Vector<Double> srcD = sourceDir.getAngle().getValue( "rad" ) ;
|
---|
482 | sourceDirectionRF.define( srcD ) ;
|
---|
483 | if ( !sysVels.empty() )
|
---|
484 | *sourceVelocityRF = sysVels[0] ;
|
---|
485 | else {
|
---|
486 | *sourceVelocityRF = (Double)0.0 ;
|
---|
487 | }
|
---|
488 | }
|
---|
489 | virtual void leaveDataDescId(const uInt recordNo, Int columnValue) {
|
---|
490 | npol = 0 ;
|
---|
491 | nchan = 0 ;
|
---|
492 | numSysCalRow = 0 ;
|
---|
493 | }
|
---|
494 | virtual void enterScanNo(const uInt recordNo, Int columnValue) {
|
---|
495 | //printf("%u: ScanNo: %d\n", recordNo, columnValue);
|
---|
496 | // put value
|
---|
497 | // scan number is 1-based in MS while 0-based in Scantable
|
---|
498 | *scanNoRF = (uInt)columnValue - 1 ;
|
---|
499 | }
|
---|
500 | virtual void leaveScanNo(const uInt recordNo, Int columnValue) {
|
---|
501 | cycleNo = 0 ;
|
---|
502 | }
|
---|
503 | virtual void enterStateId(const uInt recordNo, Int columnValue) {
|
---|
504 | //printf("%u: StateId: %d\n", recordNo, columnValue);
|
---|
505 | // SRCTYPE
|
---|
506 | Int srcType = getSrcType( columnValue ) ;
|
---|
507 |
|
---|
508 | // update header
|
---|
509 | if ( header.obstype.empty() )
|
---|
510 | getScalar( "OBS_MODE", (uInt)columnValue, statetab, header.obstype ) ;
|
---|
511 |
|
---|
512 | // put value
|
---|
513 | *sourceTypeRF = srcType ;
|
---|
514 | }
|
---|
515 | virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
|
---|
516 | virtual void enterTime(const uInt recordNo, Double columnValue) {
|
---|
517 | //printf("%u: Time: %f\n", recordNo, columnValue);
|
---|
518 | currentTime = MEpoch( Quantity( columnValue, "s" ), MEpoch::UTC ) ;
|
---|
519 |
|
---|
520 | // DIRECTION, AZEL, and SCANRATE
|
---|
521 | Vector<Double> direction, azel ;
|
---|
522 | Vector<Double> scanrate( 2, 0.0 ) ;
|
---|
523 | if ( getpt )
|
---|
524 | getDirection( direction, azel, scanrate ) ;
|
---|
525 | else
|
---|
526 | getSourceDirection( direction, azel, scanrate ) ;
|
---|
527 |
|
---|
528 | // INTERVAL
|
---|
529 | Double interval = intervalCol.asdouble( recordNo ) ;
|
---|
530 |
|
---|
531 | // WEATHER_ID
|
---|
532 | uInt wid = 0 ;
|
---|
533 | if ( isWeather )
|
---|
534 | wid = getWeatherId() ;
|
---|
535 |
|
---|
536 | // put value
|
---|
537 | Double t = currentTime.get( "d" ).getValue() ;
|
---|
538 | *timeRF = t ;
|
---|
539 | *intervalRF = interval ;
|
---|
540 | directionRF.define( direction ) ;
|
---|
541 | *azimuthRF = (Float)azel[0] ;
|
---|
542 | *elevationRF = (Float)azel[1] ;
|
---|
543 | scanRateRF.define( scanrate ) ;
|
---|
544 | *weatherIdRF = wid ;
|
---|
545 | }
|
---|
546 | virtual void leaveTime(const uInt recordNo, Double columnValue) { }
|
---|
547 | virtual Bool visitRecord(const uInt recordNo,
|
---|
548 | const Int observationId,
|
---|
549 | const Int feedId,
|
---|
550 | const Int fieldId,
|
---|
551 | const Int dataDescId,
|
---|
552 | const Int scanNo,
|
---|
553 | const Int stateId,
|
---|
554 | const Double time)
|
---|
555 | {
|
---|
556 | //printf("%u: %d, %d, %d, %d, %d, %d, %f\n", recordNo,
|
---|
557 | //observationId, feedId, fieldId, dataDescId, scanNo, stateId, time);
|
---|
558 |
|
---|
559 | // SPECTRA and FLAGTRA
|
---|
560 | //Matrix<Float> sp;
|
---|
561 | //Matrix<uChar> fl;
|
---|
562 | spectraAndFlagtra( recordNo, sp, fl ) ;
|
---|
563 |
|
---|
564 | // FLAGROW
|
---|
565 | Bool flr = flagRowCol.asBool( recordNo ) ;
|
---|
566 |
|
---|
567 | // TSYS
|
---|
568 | Matrix<Float> tsys ;
|
---|
569 | uInt scIdx = getSysCalIndex() ;
|
---|
570 | if ( numSysCalRow > 0 ) {
|
---|
571 | tsys = sysCalTsysCol( syscalRow[scIdx] ) ;
|
---|
572 | }
|
---|
573 | else {
|
---|
574 | tsys.resize( npol, 1 ) ;
|
---|
575 | tsys = 1.0 ;
|
---|
576 | }
|
---|
577 |
|
---|
578 | // TCAL_ID
|
---|
579 | Block<uInt> tcalids( npol, 0 ) ;
|
---|
580 | if ( numSysCalRow > 0 ) {
|
---|
581 | tcalids = getTcalId( syscalTime[scIdx] ) ;
|
---|
582 | }
|
---|
583 |
|
---|
584 | // put value
|
---|
585 | *cycleNoRF = cycleNo ;
|
---|
586 | *flagRowRF = (uInt)flr ;
|
---|
587 |
|
---|
588 | // for each polarization component
|
---|
589 | for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
|
---|
590 | // put value depending on polarization component
|
---|
591 | *polNoRF = polnos[ipol] ;
|
---|
592 | *tcalIdRF = tcalids[ipol] ;
|
---|
593 | spectraRF.define( sp.row( ipol ) ) ;
|
---|
594 | flagtraRF.define( fl.row( ipol ) ) ;
|
---|
595 | tsysRF.define( tsys.row( ipol ) ) ;
|
---|
596 |
|
---|
597 | // commit row
|
---|
598 | tablerow.put( rowidx ) ;
|
---|
599 | rowidx++ ;
|
---|
600 | }
|
---|
601 |
|
---|
602 | // increment CYCLENO
|
---|
603 | cycleNo++ ;
|
---|
604 |
|
---|
605 | return True ;
|
---|
606 | }
|
---|
607 | virtual void finish()
|
---|
608 | {
|
---|
609 | BaseMSFillerVisitor::finish();
|
---|
610 | //printf("Total: %u\n", count);
|
---|
611 | // remove redundant rows
|
---|
612 | //cout << "filled " << rowidx << " rows out of " << scantable.nrow() << " rows" << endl ;
|
---|
613 | if ( scantable.nrow() > rowidx ) {
|
---|
614 | uInt numRemove = scantable.nrow() - rowidx ;
|
---|
615 | //cout << "numRemove = " << numRemove << endl ;
|
---|
616 | Vector<uInt> rows( numRemove ) ;
|
---|
617 | indgen( rows, rowidx ) ;
|
---|
618 | scantable.table().removeRow( rows ) ;
|
---|
619 | }
|
---|
620 |
|
---|
621 | // antenna name and station name
|
---|
622 | String antennaName ;
|
---|
623 | getScalar( "NAME", (uInt)antennaId, anttab, antennaName ) ;
|
---|
624 | String stationName ;
|
---|
625 | getScalar( "STATION", (uInt)antennaId, anttab, stationName ) ;
|
---|
626 |
|
---|
627 | // update header
|
---|
628 | header.nif = ifmap.size() ;
|
---|
629 | header.antennaposition = antpos.get( "m" ).getValue() ;
|
---|
630 | if ( header.antennaname.empty() || header.antennaname == antennaName )
|
---|
631 | header.antennaname = antennaName ;
|
---|
632 | else
|
---|
633 | header.antennaname += "//" + antennaName ;
|
---|
634 | if ( !stationName.empty() && stationName != antennaName )
|
---|
635 | header.antennaname += "@" + stationName ;
|
---|
636 | if ( header.fluxunit.empty() || header.fluxunit == "CNTS" )
|
---|
637 | header.fluxunit = "K" ;
|
---|
638 | header.epoch = "UTC" ;
|
---|
639 | header.equinox = 2000.0 ;
|
---|
640 | if (header.freqref == "TOPO") {
|
---|
641 | header.freqref = "TOPOCENT";
|
---|
642 | } else if (header.freqref == "GEO") {
|
---|
643 | header.freqref = "GEOCENTR";
|
---|
644 | } else if (header.freqref == "BARY") {
|
---|
645 | header.freqref = "BARYCENT";
|
---|
646 | } else if (header.freqref == "GALACTO") {
|
---|
647 | header.freqref = "GALACTOC";
|
---|
648 | } else if (header.freqref == "LGROUP") {
|
---|
649 | header.freqref = "LOCALGRP";
|
---|
650 | } else if (header.freqref == "CMB") {
|
---|
651 | header.freqref = "CMBDIPOL";
|
---|
652 | } else if (header.freqref == "REST") {
|
---|
653 | header.freqref = "SOURCE";
|
---|
654 | }
|
---|
655 | scantable.setHeader( header ) ;
|
---|
656 | }
|
---|
657 | void setAntenna( Int id )
|
---|
658 | {
|
---|
659 | antennaId = id ;
|
---|
660 |
|
---|
661 | Vector< Quantum<Double> > pos ;
|
---|
662 | getArrayQuant( "POSITION", (uInt)antennaId, anttab, pos ) ;
|
---|
663 | antpos = MPosition( MVPosition( pos ), MPosition::ITRF ) ;
|
---|
664 | mf.set( antpos ) ;
|
---|
665 | }
|
---|
666 | void setPointingTable( const Table &tab, String columnToUse="DIRECTION" )
|
---|
667 | {
|
---|
668 | // input POINTING table must be
|
---|
669 | // 1) selected by antenna
|
---|
670 | // 2) sorted by TIME
|
---|
671 | ROScalarColumn<Double> tcol( tab, "TIME" ) ;
|
---|
672 | ROArrayColumn<Double> dcol( tab, columnToUse ) ;
|
---|
673 | tcol.getColumn( pointingTime ) ;
|
---|
674 | dcol.getColumn( pointingDirection ) ;
|
---|
675 | const TableRecord &rec = dcol.keywordSet() ;
|
---|
676 | String pointingRef = rec.asRecord( "MEASINFO" ).asString( "Ref" ) ;
|
---|
677 | MDirection::getType( dirType, pointingRef ) ;
|
---|
678 | getpt = True ;
|
---|
679 | }
|
---|
680 | void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it )
|
---|
681 | {
|
---|
682 | isWeather = True ;
|
---|
683 | weatherTime = t ;
|
---|
684 | weatherInterval = it ;
|
---|
685 | }
|
---|
686 | void setSysCalRecord( const Record &r )
|
---|
687 | //void setSysCalRecord( const map< String,Vector<uInt> > &r )
|
---|
688 | {
|
---|
689 | isSysCal = True ;
|
---|
690 | isTcal = True ;
|
---|
691 | syscalRecord = r ;
|
---|
692 | if ( syscalRecord.nfields() == 0 )
|
---|
693 | isTcal = False ;
|
---|
694 |
|
---|
695 | const TableDesc &desc = sctab.tableDesc() ;
|
---|
696 | uInt nrow = sctab.nrow() ;
|
---|
697 | syscalRow.resize( nrow ) ;
|
---|
698 | syscalTime.resize( nrow ) ;
|
---|
699 | syscalInterval.resize( nrow ) ;
|
---|
700 | String tsysCol = "NONE" ;
|
---|
701 | Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ;
|
---|
702 | for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) {
|
---|
703 | if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) )
|
---|
704 | tsysCol = tsysCols[i] ;
|
---|
705 | }
|
---|
706 | sysCalTsysCol.attach( sctab, tsysCol ) ;
|
---|
707 | }
|
---|
708 | STHeader getHeader() { return header ; }
|
---|
709 | uInt getNumBeam() { return nbeam ; }
|
---|
710 | uInt getFilledRowNum() { return rowidx ; }
|
---|
711 | private:
|
---|
712 | void fluxUnit( String &u )
|
---|
713 | {
|
---|
714 | ROTableColumn col( table, dataColumnName ) ;
|
---|
715 | const TableRecord &rec = col.keywordSet() ;
|
---|
716 | if ( rec.isDefined( "UNIT" ) )
|
---|
717 | u = rec.asString( "UNIT" ) ;
|
---|
718 | else if ( rec.isDefined( "QuantumUnits" ) )
|
---|
719 | u = rec.asString( "QuantumUnits" ) ;
|
---|
720 | if ( u.empty() )
|
---|
721 | u = "K" ;
|
---|
722 | }
|
---|
723 | void processSource( Int sourceId, Int spwId,
|
---|
724 | String &name, MDirection &dir, Vector<Double> &pm,
|
---|
725 | Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel )
|
---|
726 | {
|
---|
727 | // find row
|
---|
728 | uInt nrow = srctab.nrow() ;
|
---|
729 | Int idx = -1 ;
|
---|
730 | ROTableRow row( srctab ) ;
|
---|
731 | for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
|
---|
732 | const TableRecord &r = row.get( irow ) ;
|
---|
733 | if ( r.asInt( "SOURCE_ID" ) == sourceId ) {
|
---|
734 | Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ;
|
---|
735 | if ( tmpSpwId == spwId || tmpSpwId == -1 ) {
|
---|
736 | idx = (Int)irow ;
|
---|
737 | break ;
|
---|
738 | }
|
---|
739 | }
|
---|
740 | }
|
---|
741 |
|
---|
742 | // fill
|
---|
743 | Int numLines = 0 ;
|
---|
744 | if ( idx != -1 ) {
|
---|
745 | const TableRecord &r = row.get( idx ) ;
|
---|
746 | name = r.asString( "NAME" ) ;
|
---|
747 | getScalarMeas( "DIRECTION", idx, srctab, dir ) ;
|
---|
748 | pm = r.toArrayDouble( "PROPER_MOTION" ) ;
|
---|
749 | numLines = r.asInt( "NUM_LINES" ) ;
|
---|
750 | }
|
---|
751 | else {
|
---|
752 | name = "" ;
|
---|
753 | pm = Vector<Double>( 2, 0.0 ) ;
|
---|
754 | dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
|
---|
755 | }
|
---|
756 | if ( !getpt ) {
|
---|
757 | String ref = dir.getRefString() ;
|
---|
758 | MDirection::getType( dirType, ref ) ;
|
---|
759 | }
|
---|
760 |
|
---|
761 | rf.resize( numLines ) ;
|
---|
762 | trans.resize( numLines ) ;
|
---|
763 | vel.resize( numLines ) ;
|
---|
764 | if ( numLines > 0 ) {
|
---|
765 | Block<Bool> isDefined = row.getDefined() ;
|
---|
766 | Vector<String> colNames = row.columnNames() ;
|
---|
767 | Vector<Int> indexes( 3, -1 ) ;
|
---|
768 | Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ;
|
---|
769 | for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) {
|
---|
770 | if ( anyEQ( indexes, -1 ) ) {
|
---|
771 | for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) {
|
---|
772 | if ( colNames[icol] == cols[jcol] )
|
---|
773 | indexes[jcol] = icol ;
|
---|
774 | }
|
---|
775 | }
|
---|
776 | }
|
---|
777 | if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) {
|
---|
778 | Vector< Quantum<Double> > qrf ;
|
---|
779 | getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ;
|
---|
780 | for ( int i = 0 ; i < numLines ; i++ )
|
---|
781 | rf[i] = qrf[i].getValue( "Hz" ) ;
|
---|
782 | }
|
---|
783 | if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) {
|
---|
784 | getArray( "TRANSITION", idx, srctab, trans ) ;
|
---|
785 | }
|
---|
786 | if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) {
|
---|
787 | Vector< Quantum<Double> > qsv ;
|
---|
788 | getArrayQuant( "SYSVEL", idx, srctab, qsv ) ;
|
---|
789 | for ( int i = 0 ; i < numLines ; i++ )
|
---|
790 | vel[i] = qsv[i].getValue( "m/s" ) ;
|
---|
791 | }
|
---|
792 | }
|
---|
793 | }
|
---|
794 | void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md,
|
---|
795 | uInt &freqId, Int &nchan,
|
---|
796 | String &freqref, Double &reffreq, Double &bandwidth )
|
---|
797 | {
|
---|
798 | // fill
|
---|
799 | Int measFreqRef ;
|
---|
800 | getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ;
|
---|
801 | MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ;
|
---|
802 | //freqref = MFrequency::showType( freqRef ) ;
|
---|
803 | freqref = "LSRK" ;
|
---|
804 | Quantum<Double> q ;
|
---|
805 | getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ;
|
---|
806 | bandwidth = q.getValue( "Hz" ) ;
|
---|
807 | getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ;
|
---|
808 | reffreq = q.getValue( "Hz" ) ;
|
---|
809 | Double refpix = 0.5 * ( (Double)nchan-1.0 ) ;
|
---|
810 | Int refchan = ( nchan - 1 ) / 2 ;
|
---|
811 | Bool even = (Bool)( nchan % 2 == 0 ) ;
|
---|
812 | Vector< Quantum<Double> > qa ;
|
---|
813 | getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ;
|
---|
814 | Double increment = qa[refchan].getValue( "Hz" ) ;
|
---|
815 | getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ;
|
---|
816 | if ( nchan == 1 ) {
|
---|
817 | Int netSideband ;
|
---|
818 | getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ;
|
---|
819 | if ( netSideband == 1 ) increment *= -1.0 ;
|
---|
820 | }
|
---|
821 | else {
|
---|
822 | if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) )
|
---|
823 | increment *= -1.0 ;
|
---|
824 | }
|
---|
825 | Double refval = qa[refchan].getValue( "Hz" ) ;
|
---|
826 | if ( even )
|
---|
827 | refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ;
|
---|
828 | if ( freqRef != MFrequency::LSRK ) {
|
---|
829 | MeasFrame mframe( me, mp, md ) ;
|
---|
830 | MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mframe ) ) ;
|
---|
831 | refval = tolsr( Quantum<Double>( refval, "Hz" ) ).get( "Hz" ).getValue() ;
|
---|
832 | }
|
---|
833 |
|
---|
834 | // add new row to FREQUENCIES
|
---|
835 | Table ftab = scantable.frequencies().table() ;
|
---|
836 | freqId = ftab.nrow() ;
|
---|
837 | ftab.addRow() ;
|
---|
838 | TableRow row( ftab ) ;
|
---|
839 | TableRecord &r = row.record() ;
|
---|
840 | RecordFieldPtr<uInt> idRF( r, "ID" ) ;
|
---|
841 | *idRF = freqId ;
|
---|
842 | RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ;
|
---|
843 | RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ;
|
---|
844 | RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ;
|
---|
845 | *refpixRF = refpix ;
|
---|
846 | *refvalRF = refval ;
|
---|
847 | *incrRF = increment ;
|
---|
848 | row.put( freqId ) ;
|
---|
849 | }
|
---|
850 | void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl )
|
---|
851 | {
|
---|
852 | Matrix<Bool> b = flagCol( recordNo ) ;
|
---|
853 | if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) {
|
---|
854 | sp = floatDataCol( recordNo ) ;
|
---|
855 | convertArray( fl, b ) ;
|
---|
856 | }
|
---|
857 | else {
|
---|
858 | Bool notyet = True ;
|
---|
859 | Matrix<Complex> c = dataCol( recordNo ) ;
|
---|
860 | for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
|
---|
861 | if ( ( header.poltype == "linear" || header.poltype == "circular" )
|
---|
862 | && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) {
|
---|
863 | if ( notyet ) {
|
---|
864 | Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
|
---|
865 | IPosition start( 1, 0 ) ;
|
---|
866 | IPosition end( 1, 2*nchan-1 ) ;
|
---|
867 | IPosition inc( 1, 2 ) ;
|
---|
868 | if ( polnos[ipol] == 2 ) {
|
---|
869 | sp.row( ipol ) = tmp( start, end, inc ) ;
|
---|
870 | Vector<Bool> br = b.row( ipol ) ;
|
---|
871 | Vector<uChar> flr = fl.row( ipol ) ;
|
---|
872 | convertArray( flr, br ) ;
|
---|
873 | start = IPosition( 1, 1 ) ;
|
---|
874 | Int jpol = ipol+1 ;
|
---|
875 | while( polnos[jpol] != 3 && jpol < npol )
|
---|
876 | jpol++ ;
|
---|
877 | sp.row( jpol ) = tmp( start, end, inc ) ;
|
---|
878 | flr.reference( fl.row( jpol ) ) ;
|
---|
879 | convertArray( flr, br ) ;
|
---|
880 | }
|
---|
881 | else if ( polnos[ipol] == 3 ) {
|
---|
882 | sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ;
|
---|
883 | Int jpol = ipol+1 ;
|
---|
884 | while( polnos[jpol] != 2 && jpol < npol )
|
---|
885 | jpol++ ;
|
---|
886 | Vector<Bool> br = b.row( ipol ) ;
|
---|
887 | Vector<uChar> flr = fl.row( jpol ) ;
|
---|
888 | sp.row( jpol ) = tmp( start, end, inc ) ;
|
---|
889 | convertArray( flr, br ) ;
|
---|
890 | start = IPosition( 1, 1 ) ;
|
---|
891 | sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ;
|
---|
892 | flr.reference( fl.row( ipol ) ) ;
|
---|
893 | convertArray( flr, br ) ;
|
---|
894 | }
|
---|
895 | notyet = False ;
|
---|
896 | }
|
---|
897 | }
|
---|
898 | else {
|
---|
899 | Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
|
---|
900 | IPosition start( 1, 0 ) ;
|
---|
901 | IPosition end( 1, 2*nchan-1 ) ;
|
---|
902 | IPosition inc( 1, 2 ) ;
|
---|
903 | sp.row( ipol ) = tmp( start, end, inc ) ;
|
---|
904 | Vector<Bool> br = b.row( ipol ) ;
|
---|
905 | Vector<uChar> flr = fl.row( ipol ) ;
|
---|
906 | convertArray( flr, br ) ;
|
---|
907 | }
|
---|
908 | }
|
---|
909 | }
|
---|
910 | }
|
---|
911 | uInt binarySearch( Vector<Double> &timeList, Double target )
|
---|
912 | {
|
---|
913 | Int low = 0 ;
|
---|
914 | Int high = timeList.nelements() ;
|
---|
915 | uInt idx = 0 ;
|
---|
916 |
|
---|
917 | while ( low <= high ) {
|
---|
918 | idx = (Int)( 0.5 * ( low + high ) ) ;
|
---|
919 | Double t = timeList[idx] ;
|
---|
920 | if ( t < target )
|
---|
921 | low = idx + 1 ;
|
---|
922 | else if ( t > target )
|
---|
923 | high = idx - 1 ;
|
---|
924 | else {
|
---|
925 | return idx ;
|
---|
926 | }
|
---|
927 | }
|
---|
928 |
|
---|
929 | idx = max( 0, min( low, high ) ) ;
|
---|
930 | return idx ;
|
---|
931 | }
|
---|
932 | void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
|
---|
933 | {
|
---|
934 | // @todo At the moment, do binary search every time
|
---|
935 | // if this is bottleneck, frequency of binary search must be reduced
|
---|
936 | Double t = currentTime.get( "s" ).getValue() ;
|
---|
937 | uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ;
|
---|
938 | Matrix<Double> d ;
|
---|
939 | if ( pointingTime[idx] == t )
|
---|
940 | d = pointingDirection.xyPlane( idx ) ;
|
---|
941 | else if ( pointingTime[idx] < t ) {
|
---|
942 | if ( idx == pointingTime.nelements()-1 )
|
---|
943 | d = pointingDirection.xyPlane( idx ) ;
|
---|
944 | else
|
---|
945 | d = interp( pointingTime[idx], pointingTime[idx+1], t,
|
---|
946 | pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ;
|
---|
947 | }
|
---|
948 | else {
|
---|
949 | if ( idx == 0 )
|
---|
950 | d = pointingDirection.xyPlane( idx ) ;
|
---|
951 | else
|
---|
952 | d = interp( pointingTime[idx-1], pointingTime[idx], t,
|
---|
953 | pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
|
---|
954 | }
|
---|
955 | mf.resetEpoch( currentTime ) ;
|
---|
956 | Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
|
---|
957 | if ( dirType != MDirection::J2000 ) {
|
---|
958 | MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
|
---|
959 | dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
|
---|
960 | }
|
---|
961 | else {
|
---|
962 | dir = d.column( 0 ) ;
|
---|
963 | }
|
---|
964 | if ( dirType != MDirection::AZELGEO ) {
|
---|
965 | MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
|
---|
966 | //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
|
---|
967 | azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
|
---|
968 | }
|
---|
969 | else {
|
---|
970 | azel = d.column( 0 ) ;
|
---|
971 | }
|
---|
972 | if ( d.ncolumn() > 1 )
|
---|
973 | srate = d.column( 1 ) ;
|
---|
974 | }
|
---|
975 | void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
|
---|
976 | {
|
---|
977 | dir = sourceDir.getAngle( "rad" ).getValue() ;
|
---|
978 | mf.resetEpoch( currentTime ) ;
|
---|
979 | MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
|
---|
980 | azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
|
---|
981 | if ( dirType != MDirection::J2000 ) {
|
---|
982 | MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
|
---|
983 | dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
|
---|
984 | }
|
---|
985 | }
|
---|
986 | String detectSeparator( String &s )
|
---|
987 | {
|
---|
988 | String tmp = s.substr( 0, s.find_first_of( "," ) ) ;
|
---|
989 | Char *separators[] = { ":", "#", ".", "_" } ;
|
---|
990 | uInt nsep = 4 ;
|
---|
991 | for ( uInt i = 0 ; i < nsep ; i++ ) {
|
---|
992 | if ( tmp.find( separators[i] ) != String::npos )
|
---|
993 | return separators[i] ;
|
---|
994 | }
|
---|
995 | return "" ;
|
---|
996 | }
|
---|
997 | Int getSrcType( Int stateId )
|
---|
998 | {
|
---|
999 | // get values
|
---|
1000 | Bool sig ;
|
---|
1001 | getScalar( "SIG", stateId, statetab, sig ) ;
|
---|
1002 | Bool ref ;
|
---|
1003 | getScalar( "REF", stateId, statetab, ref ) ;
|
---|
1004 | Double cal ;
|
---|
1005 | getScalar( "CAL", stateId, statetab, cal ) ;
|
---|
1006 | String obsmode ;
|
---|
1007 | getScalar( "OBS_MODE", stateId, statetab, obsmode ) ;
|
---|
1008 | String sep = detectSeparator( obsmode ) ;
|
---|
1009 |
|
---|
1010 | Int srcType = SrcType::NOTYPE ;
|
---|
1011 | if ( sep == ":" )
|
---|
1012 | srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ;
|
---|
1013 | else if ( sep == "." || sep == "#" )
|
---|
1014 | srcTypeALMA( srcType, sep, obsmode ) ;
|
---|
1015 | else if ( sep == "_" )
|
---|
1016 | srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ;
|
---|
1017 | else
|
---|
1018 | srcTypeDefault( srcType, sig, ref ) ;
|
---|
1019 |
|
---|
1020 | return srcType ;
|
---|
1021 | }
|
---|
1022 | void srcTypeDefault( Int &st, Bool &sig, Bool &ref )
|
---|
1023 | {
|
---|
1024 | if ( sig ) st = SrcType::SIG ;
|
---|
1025 | else if ( ref ) st = SrcType::REF ;
|
---|
1026 | }
|
---|
1027 | void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal )
|
---|
1028 | {
|
---|
1029 | Int epos = mode.find_first_of( sep ) ;
|
---|
1030 | Int nextpos = mode.find_first_of( sep, epos+1 ) ;
|
---|
1031 | String m1 = mode.substr( 0, epos ) ;
|
---|
1032 | String m2 = mode.substr( epos+1, nextpos-epos-1 ) ;
|
---|
1033 | if ( m1 == "Nod" ) {
|
---|
1034 | st = SrcType::NOD ;
|
---|
1035 | }
|
---|
1036 | else if ( m1 == "OffOn" ) {
|
---|
1037 | if ( m2 == "PSWITCHON" ) st = SrcType::PSON ;
|
---|
1038 | if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ;
|
---|
1039 | }
|
---|
1040 | else {
|
---|
1041 | if ( m2 == "FSWITCH" ) {
|
---|
1042 | if ( sig ) st = SrcType::FSON ;
|
---|
1043 | else if ( ref ) st = SrcType::FSOFF ;
|
---|
1044 | }
|
---|
1045 | }
|
---|
1046 | if ( cal > 0.0 ) {
|
---|
1047 | if ( st == SrcType::NOD )
|
---|
1048 | st = SrcType::NODCAL ;
|
---|
1049 | else if ( st == SrcType::PSON )
|
---|
1050 | st = SrcType::PONCAL ;
|
---|
1051 | else if ( st == SrcType::PSOFF )
|
---|
1052 | st = SrcType::POFFCAL ;
|
---|
1053 | else if ( st == SrcType::FSON )
|
---|
1054 | st = SrcType::FONCAL ;
|
---|
1055 | else if ( st == SrcType::FSOFF )
|
---|
1056 | st = SrcType::FOFFCAL ;
|
---|
1057 | else
|
---|
1058 | st = SrcType::CAL ;
|
---|
1059 | }
|
---|
1060 | }
|
---|
1061 | void srcTypeALMA( Int &st, String &sep, String &mode )
|
---|
1062 | {
|
---|
1063 | Int epos = mode.find_first_of( "," ) ;
|
---|
1064 | String first = mode.substr( 0, epos ) ;
|
---|
1065 | epos = first.find_first_of( sep ) ;
|
---|
1066 | Int nextpos = first.find_first_of( sep, epos+1 ) ;
|
---|
1067 | String m1 = first.substr( 0, epos ) ;
|
---|
1068 | String m2 = first.substr( epos+1, nextpos-epos-1 ) ;
|
---|
1069 | if ( m1.find( "CALIBRATE_" ) == 0 ) {
|
---|
1070 | if ( m2.find( "ON_SOURCE" ) == 0 )
|
---|
1071 | st = SrcType::PONCAL ;
|
---|
1072 | else if ( m2.find( "OFF_SOURCE" ) == 0 )
|
---|
1073 | st = SrcType::POFFCAL ;
|
---|
1074 | }
|
---|
1075 | else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) {
|
---|
1076 | if ( m2.find( "ON_SOURCE" ) == 0 )
|
---|
1077 | st = SrcType::PSON ;
|
---|
1078 | else if ( m2.find( "OFF_SOURCE" ) == 0 )
|
---|
1079 | st = SrcType::PSOFF ;
|
---|
1080 | }
|
---|
1081 | }
|
---|
1082 | void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref )
|
---|
1083 | {
|
---|
1084 | Int epos = mode.find_first_of( "," ) ;
|
---|
1085 | String first = mode.substr( 0, epos ) ;
|
---|
1086 | string substr[4] ;
|
---|
1087 | int numSubstr = split( first, substr, 4, sep ) ;
|
---|
1088 | String m1( substr[0] ) ;
|
---|
1089 | String m2( substr[2] ) ;
|
---|
1090 | if ( numSubstr == 4 ) {
|
---|
1091 | if ( m1.find( "CALIBRATE" ) == 0 ) {
|
---|
1092 | if ( m2.find( "ON" ) == 0 )
|
---|
1093 | st = SrcType::PONCAL ;
|
---|
1094 | else if ( m2.find( "OFF" ) == 0 )
|
---|
1095 | st = SrcType::POFFCAL ;
|
---|
1096 | }
|
---|
1097 | else if ( m1.find( "OBSERVE" ) == 0 ) {
|
---|
1098 | if ( m2.find( "ON" ) == 0 )
|
---|
1099 | st = SrcType::PSON ;
|
---|
1100 | else if ( m2.find( "OFF" ) == 0 )
|
---|
1101 | st = SrcType::PSOFF ;
|
---|
1102 | }
|
---|
1103 | }
|
---|
1104 | else {
|
---|
1105 | if ( sig ) st = SrcType::SIG ;
|
---|
1106 | else if ( ref ) st = SrcType::REF ;
|
---|
1107 | }
|
---|
1108 | }
|
---|
1109 | Block<uInt> getPolNos( Vector<Int> &corr )
|
---|
1110 | {
|
---|
1111 | Block<uInt> polnos( npol ) ;
|
---|
1112 | for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
|
---|
1113 | if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX )
|
---|
1114 | polnos[ipol] = 0 ;
|
---|
1115 | else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY )
|
---|
1116 | polnos[ipol] = 1 ;
|
---|
1117 | else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY )
|
---|
1118 | polnos[ipol] = 2 ;
|
---|
1119 | else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX )
|
---|
1120 | polnos[ipol] = 3 ;
|
---|
1121 | }
|
---|
1122 | return polnos ;
|
---|
1123 | }
|
---|
1124 | String getPolType( Int &corr )
|
---|
1125 | {
|
---|
1126 | String poltype = "" ;
|
---|
1127 | if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V )
|
---|
1128 | poltype = "stokes" ;
|
---|
1129 | else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX )
|
---|
1130 | poltype = "linear" ;
|
---|
1131 | else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR )
|
---|
1132 | poltype = "circular" ;
|
---|
1133 | else if ( corr == Stokes::Plinear || corr == Stokes::Pangle )
|
---|
1134 | poltype = "linpol" ;
|
---|
1135 | return poltype ;
|
---|
1136 | }
|
---|
1137 | uInt getWeatherId()
|
---|
1138 | {
|
---|
1139 | // if only one row, return 0
|
---|
1140 | if ( weatherTime.nelements() == 1 )
|
---|
1141 | return 0 ;
|
---|
1142 |
|
---|
1143 | // @todo At the moment, do binary search every time
|
---|
1144 | // if this is bottleneck, frequency of binary search must be reduced
|
---|
1145 | Double t = currentTime.get( "s" ).getValue() ;
|
---|
1146 | uInt idx = min( binarySearch( weatherTime, t ), weatherTime.nelements()-1 ) ;
|
---|
1147 | if ( weatherTime[idx] < t ) {
|
---|
1148 | if ( idx != weatherTime.nelements()-1 ) {
|
---|
1149 | if ( weatherTime[idx+1] - t < 0.5 * weatherInterval[idx+1] )
|
---|
1150 | idx++ ;
|
---|
1151 | }
|
---|
1152 | }
|
---|
1153 | else if ( weatherTime[idx] > t ) {
|
---|
1154 | if ( idx != 0 ) {
|
---|
1155 | if ( weatherTime[idx] - t > 0.5 * weatherInterval[idx] )
|
---|
1156 | idx-- ;
|
---|
1157 | }
|
---|
1158 | }
|
---|
1159 | return idx ;
|
---|
1160 | }
|
---|
1161 | void processSysCal( Int &spwId )
|
---|
1162 | {
|
---|
1163 | // get feedId from row
|
---|
1164 | Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
|
---|
1165 |
|
---|
1166 | uInt nrow = sctab.nrow() ;
|
---|
1167 | ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ;
|
---|
1168 | Vector<Int> aids = col.getColumn() ;
|
---|
1169 | col.attach( sctab, "FEED_ID" ) ;
|
---|
1170 | Vector<Int> fids = col.getColumn() ;
|
---|
1171 | col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ;
|
---|
1172 | Vector<Int> sids = col.getColumn() ;
|
---|
1173 | ROScalarColumn<Double> timeCol( sctab, "TIME" ) ;
|
---|
1174 | ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ;
|
---|
1175 | for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
|
---|
1176 | if ( aids[irow] == antennaId
|
---|
1177 | && fids[irow] == feedId
|
---|
1178 | && sids[irow] == spwId ) {
|
---|
1179 | syscalRow[numSysCalRow] = irow ;
|
---|
1180 | syscalTime[numSysCalRow] = timeCol( irow ) ;
|
---|
1181 | syscalInterval[numSysCalRow] = intCol( irow ) ;
|
---|
1182 | numSysCalRow++ ;
|
---|
1183 | }
|
---|
1184 | }
|
---|
1185 | }
|
---|
1186 | uInt getSysCalIndex()
|
---|
1187 | {
|
---|
1188 | // if only one row, return 0
|
---|
1189 | if ( numSysCalRow == 1 || !isSysCal )
|
---|
1190 | return 0 ;
|
---|
1191 |
|
---|
1192 | // @todo At the moment, do binary search every time
|
---|
1193 | // if this is bottleneck, frequency of binary search must be reduced
|
---|
1194 | Double t = currentTime.get( "s" ).getValue() ;
|
---|
1195 | Vector<Double> tslice = syscalTime( Slice(0, numSysCalRow) ) ;
|
---|
1196 | uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ;
|
---|
1197 | if ( syscalTime[idx] < t ) {
|
---|
1198 | if ( idx != numSysCalRow-1 ) {
|
---|
1199 | if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] )
|
---|
1200 | idx++ ;
|
---|
1201 | }
|
---|
1202 | }
|
---|
1203 | else if ( syscalTime[idx] > t ) {
|
---|
1204 | if ( idx != 0 ) {
|
---|
1205 | if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] )
|
---|
1206 | idx-- ;
|
---|
1207 | }
|
---|
1208 | }
|
---|
1209 | return idx ;
|
---|
1210 | }
|
---|
1211 | Block<uInt> getTcalId( Double &t )
|
---|
1212 | {
|
---|
1213 | // return 0 if no SysCal table
|
---|
1214 | if ( !isSysCal or !isTcal ) {
|
---|
1215 | return Block<uInt>( 4, 0 ) ;
|
---|
1216 | }
|
---|
1217 |
|
---|
1218 | // get feedId from row
|
---|
1219 | Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
|
---|
1220 |
|
---|
1221 | // key
|
---|
1222 | String key = keyTcal( feedId, spwId, t ) ;
|
---|
1223 |
|
---|
1224 | // retrieve ids
|
---|
1225 | Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ;
|
---|
1226 | //Vector<uInt> ids = syscalRecord[key] ;
|
---|
1227 | uInt np = ids[1] - ids[0] + 1 ;
|
---|
1228 | Block<uInt> tcalids( np ) ;
|
---|
1229 | if ( np > 0 ) {
|
---|
1230 | tcalids[0] = ids[0] ;
|
---|
1231 | if ( np > 1 ) {
|
---|
1232 | tcalids[1] = ids[1] ;
|
---|
1233 | for ( uInt ip = 2 ; ip < np ; ip++ )
|
---|
1234 | tcalids[ip] = ids[0] + ip - 1 ;
|
---|
1235 | }
|
---|
1236 | }
|
---|
1237 | return tcalids ;
|
---|
1238 | }
|
---|
1239 | uInt maxNumPol()
|
---|
1240 | {
|
---|
1241 | ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ;
|
---|
1242 | return max( numCorrCol.getColumn() ) ;
|
---|
1243 | }
|
---|
1244 |
|
---|
1245 | Scantable &scantable;
|
---|
1246 | Int antennaId;
|
---|
1247 | uInt rowidx;
|
---|
1248 | String dataColumnName;
|
---|
1249 | TableRow tablerow;
|
---|
1250 | STHeader header;
|
---|
1251 | Vector<Int> feedEntry;
|
---|
1252 | uInt nbeam;
|
---|
1253 | Int npol;
|
---|
1254 | Int nchan;
|
---|
1255 | Int sourceId;
|
---|
1256 | Int polId;
|
---|
1257 | Int spwId;
|
---|
1258 | uInt cycleNo;
|
---|
1259 | MDirection sourceDir;
|
---|
1260 | MPosition antpos;
|
---|
1261 | MEpoch currentTime;
|
---|
1262 | map<Int,uInt> ifmap;
|
---|
1263 | Block<uInt> polnos;
|
---|
1264 | Bool getpt;
|
---|
1265 | Vector<Double> pointingTime;
|
---|
1266 | Cube<Double> pointingDirection;
|
---|
1267 | MDirection::Types dirType;
|
---|
1268 | Bool isWeather;
|
---|
1269 | Vector<Double> weatherTime;
|
---|
1270 | Vector<Double> weatherInterval;
|
---|
1271 | Bool isSysCal;
|
---|
1272 | Bool isTcal;
|
---|
1273 | Record syscalRecord;
|
---|
1274 | //map< String,Vector<uInt> > syscalRecord;
|
---|
1275 | uInt numSysCalRow ;
|
---|
1276 | Vector<uInt> syscalRow;
|
---|
1277 | Vector<Double> syscalTime;
|
---|
1278 | Vector<Double> syscalInterval;
|
---|
1279 | //String tsysCol;
|
---|
1280 | //String tcalCol;
|
---|
1281 |
|
---|
1282 | // MS subtables
|
---|
1283 | Table obstab;
|
---|
1284 | Table sctab;
|
---|
1285 | Table spwtab;
|
---|
1286 | Table statetab;
|
---|
1287 | Table ddtab;
|
---|
1288 | Table poltab;
|
---|
1289 | Table fieldtab;
|
---|
1290 | Table anttab;
|
---|
1291 | Table srctab;
|
---|
1292 | Matrix<Float> sp;
|
---|
1293 | Matrix<uChar> fl;
|
---|
1294 | MeasFrame mf;
|
---|
1295 |
|
---|
1296 | // MS MAIN columns
|
---|
1297 | ROTableColumn intervalCol;
|
---|
1298 | ROTableColumn flagRowCol;
|
---|
1299 | ROArrayColumn<Float> floatDataCol;
|
---|
1300 | ROArrayColumn<Complex> dataCol;
|
---|
1301 | ROArrayColumn<Bool> flagCol;
|
---|
1302 |
|
---|
1303 | // MS SYSCAL columns
|
---|
1304 | ROArrayColumn<Float> sysCalTsysCol;
|
---|
1305 |
|
---|
1306 | // Scantable MAIN columns
|
---|
1307 | RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF;
|
---|
1308 | RecordFieldPtr< Vector<Double> > directionRF,scanRateRF,
|
---|
1309 | sourceProperMotionRF,sourceDirectionRF;
|
---|
1310 | RecordFieldPtr<Float> azimuthRF,elevationRF;
|
---|
1311 | RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF,
|
---|
1312 | ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF;
|
---|
1313 | RecordFieldPtr< Vector<Float> > spectraRF,tsysRF;
|
---|
1314 | RecordFieldPtr< Vector<uChar> > flagtraRF;
|
---|
1315 | RecordFieldPtr<String> sourceNameRF,fieldNameRF;
|
---|
1316 | RecordFieldPtr<Int> sourceTypeRF;
|
---|
1317 | };
|
---|
1318 |
|
---|
1319 | class BaseTcalVisitor: public TableVisitor {
|
---|
1320 | uInt lastRecordNo ;
|
---|
1321 | Int lastAntennaId ;
|
---|
1322 | Int lastFeedId ;
|
---|
1323 | Int lastSpwId ;
|
---|
1324 | Double lastTime ;
|
---|
1325 | protected:
|
---|
1326 | const Table &table;
|
---|
1327 | uInt count;
|
---|
1328 | public:
|
---|
1329 | BaseTcalVisitor(const Table &table)
|
---|
1330 | : table(table)
|
---|
1331 | {
|
---|
1332 | count = 0;
|
---|
1333 | }
|
---|
1334 |
|
---|
1335 | virtual void enterAntennaId(const uInt recordNo, Int columnValue) { }
|
---|
1336 | virtual void leaveAntennaId(const uInt recordNo, Int columnValue) { }
|
---|
1337 | virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
|
---|
1338 | virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
|
---|
1339 | virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
|
---|
1340 | virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
|
---|
1341 | virtual void enterTime(const uInt recordNo, Double columnValue) { }
|
---|
1342 | virtual void leaveTime(const uInt recordNo, Double columnValue) { }
|
---|
1343 |
|
---|
1344 | virtual Bool visitRecord(const uInt recordNo,
|
---|
1345 | const Int antennaId,
|
---|
1346 | const Int feedId,
|
---|
1347 | const Int spwId,
|
---|
1348 | const Double time) { return True ; }
|
---|
1349 |
|
---|
1350 | virtual Bool visit(Bool isFirst, const uInt recordNo,
|
---|
1351 | const uInt nCols, void const *const colValues[]) {
|
---|
1352 | Int antennaId, feedId, spwId;
|
---|
1353 | Double time;
|
---|
1354 | { // prologue
|
---|
1355 | uInt i = 0;
|
---|
1356 | {
|
---|
1357 | const Int *col = (const Int *)colValues[i++];
|
---|
1358 | antennaId = col[recordNo];
|
---|
1359 | }
|
---|
1360 | {
|
---|
1361 | const Int *col = (const Int *)colValues[i++];
|
---|
1362 | feedId = col[recordNo];
|
---|
1363 | }
|
---|
1364 | {
|
---|
1365 | const Int *col = (const Int *)colValues[i++];
|
---|
1366 | spwId = col[recordNo];
|
---|
1367 | }
|
---|
1368 | {
|
---|
1369 | const Double *col = (const Double *)colValues[i++];
|
---|
1370 | time = col[recordNo];
|
---|
1371 | }
|
---|
1372 | assert(nCols == i);
|
---|
1373 | }
|
---|
1374 |
|
---|
1375 | if (isFirst) {
|
---|
1376 | enterAntennaId(recordNo, antennaId);
|
---|
1377 | enterFeedId(recordNo, feedId);
|
---|
1378 | enterSpwId(recordNo, spwId);
|
---|
1379 | enterTime(recordNo, time);
|
---|
1380 | } else {
|
---|
1381 | if ( lastAntennaId != antennaId ) {
|
---|
1382 | leaveTime(lastRecordNo, lastTime);
|
---|
1383 | leaveSpwId(lastRecordNo, lastSpwId);
|
---|
1384 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
1385 | leaveAntennaId(lastRecordNo, lastAntennaId);
|
---|
1386 |
|
---|
1387 | enterAntennaId(recordNo, antennaId);
|
---|
1388 | enterFeedId(recordNo, feedId);
|
---|
1389 | enterSpwId(recordNo, spwId);
|
---|
1390 | enterTime(recordNo, time);
|
---|
1391 | }
|
---|
1392 | else if (lastFeedId != feedId) {
|
---|
1393 | leaveTime(lastRecordNo, lastTime);
|
---|
1394 | leaveSpwId(lastRecordNo, lastSpwId);
|
---|
1395 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
1396 |
|
---|
1397 | enterFeedId(recordNo, feedId);
|
---|
1398 | enterSpwId(recordNo, spwId);
|
---|
1399 | enterTime(recordNo, time);
|
---|
1400 | } else if (lastSpwId != spwId) {
|
---|
1401 | leaveTime(lastRecordNo, lastTime);
|
---|
1402 | leaveSpwId(lastRecordNo, lastSpwId);
|
---|
1403 |
|
---|
1404 | enterSpwId(recordNo, spwId);
|
---|
1405 | enterTime(recordNo, time);
|
---|
1406 | } else if (lastTime != time) {
|
---|
1407 | leaveTime(lastRecordNo, lastTime);
|
---|
1408 | enterTime(recordNo, time);
|
---|
1409 | }
|
---|
1410 | }
|
---|
1411 | count++;
|
---|
1412 | Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time);
|
---|
1413 |
|
---|
1414 | { // epilogue
|
---|
1415 | lastRecordNo = recordNo;
|
---|
1416 |
|
---|
1417 | lastAntennaId = antennaId;
|
---|
1418 | lastFeedId = feedId;
|
---|
1419 | lastSpwId = spwId;
|
---|
1420 | lastTime = time;
|
---|
1421 | }
|
---|
1422 | return result ;
|
---|
1423 | }
|
---|
1424 |
|
---|
1425 | virtual void finish() {
|
---|
1426 | if (count > 0) {
|
---|
1427 | leaveTime(lastRecordNo, lastTime);
|
---|
1428 | leaveSpwId(lastRecordNo, lastSpwId);
|
---|
1429 | leaveFeedId(lastRecordNo, lastFeedId);
|
---|
1430 | leaveAntennaId(lastRecordNo, lastAntennaId);
|
---|
1431 | }
|
---|
1432 | }
|
---|
1433 | };
|
---|
1434 |
|
---|
1435 | class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
|
---|
1436 | public:
|
---|
1437 | TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid )
|
---|
1438 | //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid )
|
---|
1439 | : BaseTcalVisitor( table ),
|
---|
1440 | tcal(tcaltab),
|
---|
1441 | rec(r),
|
---|
1442 | antenna(aid)
|
---|
1443 | {
|
---|
1444 | process = False ;
|
---|
1445 | rowidx = 0 ;
|
---|
1446 |
|
---|
1447 | // attach to SYSCAL columns
|
---|
1448 | timeCol.attach( table, "TIME" ) ;
|
---|
1449 |
|
---|
1450 | // add rows
|
---|
1451 | uInt addrow = table.nrow() * 4 ;
|
---|
1452 | tcal.addRow( addrow ) ;
|
---|
1453 |
|
---|
1454 | // attach to TCAL columns
|
---|
1455 | row = TableRow( tcal ) ;
|
---|
1456 | TableRecord &trec = row.record() ;
|
---|
1457 | idRF.attachToRecord( trec, "ID" ) ;
|
---|
1458 | timeRF.attachToRecord( trec, "TIME" ) ;
|
---|
1459 | tcalRF.attachToRecord( trec, "TCAL" ) ;
|
---|
1460 | }
|
---|
1461 |
|
---|
1462 | virtual void enterAntennaId(const uInt recordNo, Int columnValue) {
|
---|
1463 | if ( columnValue == antenna )
|
---|
1464 | process = True ;
|
---|
1465 | }
|
---|
1466 | virtual void leaveAntennaId(const uInt recordNo, Int columnValue) {
|
---|
1467 | process = False ;
|
---|
1468 | }
|
---|
1469 | virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
|
---|
1470 | virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
|
---|
1471 | virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
|
---|
1472 | virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
|
---|
1473 | virtual void enterTime(const uInt recordNo, Double columnValue) {
|
---|
1474 | qtime = timeCol( recordNo ) ;
|
---|
1475 | }
|
---|
1476 | virtual void leaveTime(const uInt recordNo, Double columnValue) { }
|
---|
1477 | virtual Bool visitRecord(const uInt recordNo,
|
---|
1478 | const Int antennaId,
|
---|
1479 | const Int feedId,
|
---|
1480 | const Int spwId,
|
---|
1481 | const Double time)
|
---|
1482 | {
|
---|
1483 | //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ;
|
---|
1484 | if ( process ) {
|
---|
1485 | String sTime = MVTime( qtime ).string( MVTime::YMD ) ;
|
---|
1486 | *timeRF = sTime ;
|
---|
1487 | uInt oldidx = rowidx ;
|
---|
1488 | Matrix<Float> subtcal = tcalCol( recordNo ) ;
|
---|
1489 | Vector<uInt> idminmax( 2 ) ;
|
---|
1490 | for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) {
|
---|
1491 | *idRF = rowidx ;
|
---|
1492 | tcalRF.define( subtcal.row( ipol ) ) ;
|
---|
1493 |
|
---|
1494 | // commit row
|
---|
1495 | row.put( rowidx ) ;
|
---|
1496 | rowidx++ ;
|
---|
1497 | }
|
---|
1498 |
|
---|
1499 | idminmax[0] = oldidx ;
|
---|
1500 | idminmax[1] = rowidx - 1 ;
|
---|
1501 |
|
---|
1502 | String key = keyTcal( feedId, spwId, sTime ) ;
|
---|
1503 | rec.define( key, idminmax ) ;
|
---|
1504 | //rec[key] = idminmax ;
|
---|
1505 | }
|
---|
1506 | return True ;
|
---|
1507 | }
|
---|
1508 | virtual void finish()
|
---|
1509 | {
|
---|
1510 | BaseTcalVisitor::finish() ;
|
---|
1511 |
|
---|
1512 | if ( tcal.nrow() > rowidx ) {
|
---|
1513 | uInt numRemove = tcal.nrow() - rowidx ;
|
---|
1514 | //cout << "numRemove = " << numRemove << endl ;
|
---|
1515 | Vector<uInt> rows( numRemove ) ;
|
---|
1516 | indgen( rows, rowidx ) ;
|
---|
1517 | tcal.removeRow( rows ) ;
|
---|
1518 | }
|
---|
1519 |
|
---|
1520 | }
|
---|
1521 | void setTcalColumn( String &col )
|
---|
1522 | {
|
---|
1523 | //colName = col ;
|
---|
1524 | tcalCol.attach( table, col ) ;
|
---|
1525 | }
|
---|
1526 | private:
|
---|
1527 | Table &tcal;
|
---|
1528 | Record &rec;
|
---|
1529 | //map< String,Vector<uInt> > &rec;
|
---|
1530 | Int antenna;
|
---|
1531 | uInt rowidx;
|
---|
1532 | Bool process;
|
---|
1533 | Quantum<Double> qtime;
|
---|
1534 | TableRow row;
|
---|
1535 | String colName;
|
---|
1536 |
|
---|
1537 | // MS SYSCAL columns
|
---|
1538 | ROScalarQuantColumn<Double> timeCol;
|
---|
1539 | ROArrayColumn<Float> tcalCol;
|
---|
1540 |
|
---|
1541 | // TCAL columns
|
---|
1542 | RecordFieldPtr<uInt> idRF;
|
---|
1543 | RecordFieldPtr<String> timeRF;
|
---|
1544 | RecordFieldPtr< Vector<Float> > tcalRF;
|
---|
1545 | };
|
---|
1546 |
|
---|
1547 | MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
|
---|
1548 | : table_( stable ),
|
---|
1549 | tablename_( "" ),
|
---|
1550 | antenna_( -1 ),
|
---|
1551 | antennaStr_(""),
|
---|
1552 | getPt_( True ),
|
---|
1553 | isFloatData_( False ),
|
---|
1554 | isData_( False ),
|
---|
1555 | isDoppler_( False ),
|
---|
1556 | isFlagCmd_( False ),
|
---|
1557 | isFreqOffset_( False ),
|
---|
1558 | isHistory_( False ),
|
---|
1559 | isProcessor_( False ),
|
---|
1560 | isSysCal_( False ),
|
---|
1561 | isWeather_( False ),
|
---|
1562 | colTsys_( "TSYS_SPECTRUM" ),
|
---|
1563 | colTcal_( "TCAL_SPECTRUM" )
|
---|
1564 | {
|
---|
1565 | os_ = LogIO() ;
|
---|
1566 | os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
|
---|
1567 | }
|
---|
1568 |
|
---|
1569 | MSFiller::~MSFiller()
|
---|
1570 | {
|
---|
1571 | os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
|
---|
1572 | }
|
---|
1573 |
|
---|
1574 | bool MSFiller::open( const std::string &filename, const casa::Record &rec )
|
---|
1575 | {
|
---|
1576 | os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
|
---|
1577 | //double startSec = mathutil::gettimeofday_sec() ;
|
---|
1578 | //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
|
---|
1579 | //os_ << " filename = " << filename << endl ;
|
---|
1580 |
|
---|
1581 | // parsing MS options
|
---|
1582 | if ( rec.isDefined( "ms" ) ) {
|
---|
1583 | Record msrec = rec.asRecord( "ms" ) ;
|
---|
1584 | if ( msrec.isDefined( "getpt" ) ) {
|
---|
1585 | getPt_ = msrec.asBool( "getpt" ) ;
|
---|
1586 | }
|
---|
1587 | if ( msrec.isDefined( "antenna" ) ) {
|
---|
1588 | if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
|
---|
1589 | antenna_ = msrec.asInt( "antenna" ) ;
|
---|
1590 | }
|
---|
1591 | else {
|
---|
1592 | //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
|
---|
1593 | antennaStr_ = msrec.asString( "antenna" ) ;
|
---|
1594 | }
|
---|
1595 | }
|
---|
1596 | else {
|
---|
1597 | antenna_ = 0 ;
|
---|
1598 | }
|
---|
1599 | }
|
---|
1600 |
|
---|
1601 | MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
|
---|
1602 | tablename_ = tmpMS->tableName() ;
|
---|
1603 | if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
|
---|
1604 | MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
|
---|
1605 | Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
|
---|
1606 | if ( id.size() > 0 )
|
---|
1607 | antenna_ = id[0] ;
|
---|
1608 | }
|
---|
1609 |
|
---|
1610 | os_ << "Parsing MS options" << endl ;
|
---|
1611 | os_ << " getPt = " << getPt_ << endl ;
|
---|
1612 | os_ << " antenna = " << antenna_ << endl ;
|
---|
1613 | os_ << " antennaStr = " << antennaStr_ << LogIO::POST ;
|
---|
1614 |
|
---|
1615 | mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
|
---|
1616 | && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
|
---|
1617 |
|
---|
1618 | delete tmpMS ;
|
---|
1619 |
|
---|
1620 | // check which data column exists
|
---|
1621 | isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
|
---|
1622 | isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
|
---|
1623 |
|
---|
1624 | //double endSec = mathutil::gettimeofday_sec() ;
|
---|
1625 | //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
|
---|
1626 | return true ;
|
---|
1627 | }
|
---|
1628 |
|
---|
1629 | void MSFiller::fill()
|
---|
1630 | {
|
---|
1631 | //double startSec = mathutil::gettimeofday_sec() ;
|
---|
1632 | //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
|
---|
1633 |
|
---|
1634 | os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
|
---|
1635 |
|
---|
1636 | // Initialize header
|
---|
1637 | STHeader sdh ;
|
---|
1638 | initHeader( sdh ) ;
|
---|
1639 | table_->setHeader( sdh ) ;
|
---|
1640 |
|
---|
1641 | // check if optional table exists
|
---|
1642 | const TableRecord &msrec = mstable_.keywordSet() ;
|
---|
1643 | isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
|
---|
1644 | if ( isDoppler_ )
|
---|
1645 | if ( mstable_.doppler().nrow() == 0 )
|
---|
1646 | isDoppler_ = False ;
|
---|
1647 | isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
|
---|
1648 | if ( isFlagCmd_ )
|
---|
1649 | if ( mstable_.flagCmd().nrow() == 0 )
|
---|
1650 | isFlagCmd_ = False ;
|
---|
1651 | isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
|
---|
1652 | if ( isFreqOffset_ )
|
---|
1653 | if ( mstable_.freqOffset().nrow() == 0 )
|
---|
1654 | isFreqOffset_ = False ;
|
---|
1655 | isHistory_ = msrec.isDefined( "HISTORY" ) ;
|
---|
1656 | if ( isHistory_ )
|
---|
1657 | if ( mstable_.history().nrow() == 0 )
|
---|
1658 | isHistory_ = False ;
|
---|
1659 | isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
|
---|
1660 | if ( isProcessor_ )
|
---|
1661 | if ( mstable_.processor().nrow() == 0 )
|
---|
1662 | isProcessor_ = False ;
|
---|
1663 | isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
|
---|
1664 | if ( isSysCal_ )
|
---|
1665 | if ( mstable_.sysCal().nrow() == 0 )
|
---|
1666 | isSysCal_ = False ;
|
---|
1667 | isWeather_ = msrec.isDefined( "WEATHER" ) ;
|
---|
1668 | if ( isWeather_ )
|
---|
1669 | if ( mstable_.weather().nrow() == 0 )
|
---|
1670 | isWeather_ = False ;
|
---|
1671 |
|
---|
1672 | // column name for Tsys and Tcal
|
---|
1673 | if ( isSysCal_ ) {
|
---|
1674 | const MSSysCal &caltab = mstable_.sysCal() ;
|
---|
1675 | if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
|
---|
1676 | colTcal_ = "TCAL" ;
|
---|
1677 | if ( !caltab.tableDesc().isColumn( colTcal_ ) )
|
---|
1678 | colTcal_ = "NONE" ;
|
---|
1679 | }
|
---|
1680 | if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
|
---|
1681 | colTsys_ = "TSYS" ;
|
---|
1682 | if ( !caltab.tableDesc().isColumn( colTcal_ ) )
|
---|
1683 | colTsys_ = "NONE" ;
|
---|
1684 | }
|
---|
1685 | }
|
---|
1686 | else {
|
---|
1687 | colTcal_ = "NONE" ;
|
---|
1688 | colTsys_ = "NONE" ;
|
---|
1689 | }
|
---|
1690 |
|
---|
1691 | // Access to MS subtables
|
---|
1692 | //MSField &fieldtab = mstable_.field() ;
|
---|
1693 | //MSPolarization &poltab = mstable_.polarization() ;
|
---|
1694 | //MSDataDescription &ddtab = mstable_.dataDescription() ;
|
---|
1695 | //MSObservation &obstab = mstable_.observation() ;
|
---|
1696 | //MSSource &srctab = mstable_.source() ;
|
---|
1697 | //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
|
---|
1698 | //MSSysCal &caltab = mstable_.sysCal() ;
|
---|
1699 | MSPointing &pointtab = mstable_.pointing() ;
|
---|
1700 | //MSState &stattab = mstable_.state() ;
|
---|
1701 | //MSAntenna &anttab = mstable_.antenna() ;
|
---|
1702 |
|
---|
1703 | // SUBTABLES: FREQUENCIES
|
---|
1704 | //string freqFrame = getFrame() ;
|
---|
1705 | string freqFrame = "LSRK" ;
|
---|
1706 | table_->frequencies().setFrame( freqFrame ) ;
|
---|
1707 | table_->frequencies().setFrame( freqFrame, True ) ;
|
---|
1708 |
|
---|
1709 | // SUBTABLES: WEATHER
|
---|
1710 | fillWeather() ;
|
---|
1711 |
|
---|
1712 | // SUBTABLES: FOCUS
|
---|
1713 | fillFocus() ;
|
---|
1714 |
|
---|
1715 | // SUBTABLES: TCAL
|
---|
1716 | fillTcal() ;
|
---|
1717 |
|
---|
1718 | // SUBTABLES: FIT
|
---|
1719 | //fillFit() ;
|
---|
1720 |
|
---|
1721 | // SUBTABLES: HISTORY
|
---|
1722 | //fillHistory() ;
|
---|
1723 |
|
---|
1724 | /***
|
---|
1725 | * Start iteration using TableVisitor
|
---|
1726 | ***/
|
---|
1727 | Table stab = table_->table() ;
|
---|
1728 | {
|
---|
1729 | static const char *cols[] = {
|
---|
1730 | "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
|
---|
1731 | "STATE_ID", "TIME",
|
---|
1732 | NULL
|
---|
1733 | };
|
---|
1734 | static const TypeManagerImpl<Int> tmInt;
|
---|
1735 | static const TypeManagerImpl<Double> tmDouble;
|
---|
1736 | static const TypeManager *const tms[] = {
|
---|
1737 | &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
|
---|
1738 | };
|
---|
1739 | //double t0 = mathutil::gettimeofday_sec() ;
|
---|
1740 | MSFillerVisitor myVisitor(mstable_, *table_ );
|
---|
1741 | //double t1 = mathutil::gettimeofday_sec() ;
|
---|
1742 | //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
|
---|
1743 | myVisitor.setAntenna( antenna_ ) ;
|
---|
1744 | //myVisitor.setHeader( sdh ) ;
|
---|
1745 | if ( getPt_ ) {
|
---|
1746 | Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
|
---|
1747 | myVisitor.setPointingTable( ptsel ) ;
|
---|
1748 | }
|
---|
1749 | if ( isWeather_ )
|
---|
1750 | myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ;
|
---|
1751 | if ( isSysCal_ )
|
---|
1752 | myVisitor.setSysCalRecord( tcalrec_ ) ;
|
---|
1753 |
|
---|
1754 | //double t2 = mathutil::gettimeofday_sec() ;
|
---|
1755 | traverseTable(mstable_, cols, tms, &myVisitor);
|
---|
1756 | //double t3 = mathutil::gettimeofday_sec() ;
|
---|
1757 | //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
|
---|
1758 |
|
---|
1759 | sdh = myVisitor.getHeader() ;
|
---|
1760 | }
|
---|
1761 | /***
|
---|
1762 | * End iteration using TableVisitor
|
---|
1763 | ***/
|
---|
1764 |
|
---|
1765 | // set header
|
---|
1766 | //sdh = myVisitor.getHeader() ;
|
---|
1767 | //table_->setHeader( sdh ) ;
|
---|
1768 |
|
---|
1769 | // save path to POINTING table
|
---|
1770 | // 2011/07/06 TN
|
---|
1771 | // Path to POINTING table in original MS will not be written
|
---|
1772 | // if getPt_ is True
|
---|
1773 | Path datapath( tablename_ ) ;
|
---|
1774 | if ( !getPt_ ) {
|
---|
1775 | String pTabName = datapath.absoluteName() + "/POINTING" ;
|
---|
1776 | stab.rwKeywordSet().define( "POINTING", pTabName ) ;
|
---|
1777 | }
|
---|
1778 |
|
---|
1779 | // for GBT
|
---|
1780 | if ( sdh.antennaname.contains( "GBT" ) ) {
|
---|
1781 | String goTabName = datapath.absoluteName() + "/GBT_GO" ;
|
---|
1782 | stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
|
---|
1783 | }
|
---|
1784 |
|
---|
1785 | // for MS created from ASDM
|
---|
1786 | const TableRecord &msKeys = mstable_.keywordSet() ;
|
---|
1787 | uInt nfields = msKeys.nfields() ;
|
---|
1788 | for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
|
---|
1789 | String name = msKeys.name( ifield ) ;
|
---|
1790 | //os_ << "name = " << name << LogIO::POST ;
|
---|
1791 | if ( name.find( "ASDM" ) != String::npos ) {
|
---|
1792 | String asdmpath = msKeys.asTable( ifield ).tableName() ;
|
---|
1793 | os_ << "ASDM table: " << asdmpath << LogIO::POST ;
|
---|
1794 | stab.rwKeywordSet().define( name, asdmpath ) ;
|
---|
1795 | }
|
---|
1796 | }
|
---|
1797 |
|
---|
1798 | //double endSec = mathutil::gettimeofday_sec() ;
|
---|
1799 | //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
|
---|
1800 | }
|
---|
1801 |
|
---|
1802 | void MSFiller::close()
|
---|
1803 | {
|
---|
1804 | //tablesel_.closeSubTables() ;
|
---|
1805 | mstable_.closeSubTables() ;
|
---|
1806 | //tablesel_.unlock() ;
|
---|
1807 | mstable_.unlock() ;
|
---|
1808 | }
|
---|
1809 |
|
---|
1810 | void MSFiller::fillWeather()
|
---|
1811 | {
|
---|
1812 | //double startSec = mathutil::gettimeofday_sec() ;
|
---|
1813 | //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
|
---|
1814 |
|
---|
1815 | if ( !isWeather_ ) {
|
---|
1816 | // add dummy row
|
---|
1817 | table_->weather().table().addRow(1,True) ;
|
---|
1818 | return ;
|
---|
1819 | }
|
---|
1820 |
|
---|
1821 | Table mWeather = mstable_.weather() ;
|
---|
1822 | //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
|
---|
1823 | Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
|
---|
1824 | //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
|
---|
1825 | if ( mWeatherSel.nrow() == 0 ) {
|
---|
1826 | os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
|
---|
1827 | mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
|
---|
1828 | if ( mWeatherSel.nrow() == 0 ) {
|
---|
1829 | os_ << "No rows in WEATHER table" << LogIO::POST ;
|
---|
1830 | }
|
---|
1831 | }
|
---|
1832 | uInt wnrow = mWeatherSel.nrow() ;
|
---|
1833 | //os_ << "wnrow = " << wnrow << LogIO::POST ;
|
---|
1834 |
|
---|
1835 | if ( wnrow == 0 )
|
---|
1836 | return ;
|
---|
1837 |
|
---|
1838 | Table wtab = table_->weather().table() ;
|
---|
1839 | wtab.addRow( wnrow ) ;
|
---|
1840 |
|
---|
1841 | Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
|
---|
1842 | Int stationId = -1 ;
|
---|
1843 | if ( stationInfoExists ) {
|
---|
1844 | // determine which station is closer
|
---|
1845 | ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
|
---|
1846 | ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
|
---|
1847 | Vector<Int> stationIds = stationCol.getColumn() ;
|
---|
1848 | Vector<Int> stationIdList( 0 ) ;
|
---|
1849 | Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
|
---|
1850 | uInt numStation = 0 ;
|
---|
1851 | for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
|
---|
1852 | if ( !anyEQ( stationIdList, stationIds[i] ) ) {
|
---|
1853 | numStation++ ;
|
---|
1854 | stationIdList.resize( numStation, True ) ;
|
---|
1855 | stationIdList[numStation-1] = stationIds[i] ;
|
---|
1856 | stationPosList.resize( numStation, 3, True ) ;
|
---|
1857 | stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
|
---|
1858 | }
|
---|
1859 | }
|
---|
1860 | //os_ << "staionIdList = " << stationIdList << endl ;
|
---|
1861 | Table mAntenna = mstable_.antenna() ;
|
---|
1862 | ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
|
---|
1863 | Vector<Double> antpos = antposCol( antenna_ ) ;
|
---|
1864 | Double minDiff = -1.0 ;
|
---|
1865 | for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
|
---|
1866 | Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
|
---|
1867 | if ( minDiff < 0.0 || minDiff > diff ) {
|
---|
1868 | minDiff = diff ;
|
---|
1869 | stationId = stationIdList[i] ;
|
---|
1870 | }
|
---|
1871 | }
|
---|
1872 | }
|
---|
1873 | //os_ << "stationId = " << stationId << endl ;
|
---|
1874 |
|
---|
1875 | ScalarColumn<Float> *fCol ;
|
---|
1876 | ROScalarColumn<Float> *sharedFloatCol ;
|
---|
1877 | if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
|
---|
1878 | fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
|
---|
1879 | sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
|
---|
1880 | fCol->putColumn( *sharedFloatCol ) ;
|
---|
1881 | delete sharedFloatCol ;
|
---|
1882 | delete fCol ;
|
---|
1883 | }
|
---|
1884 | if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
|
---|
1885 | fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
|
---|
1886 | sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
|
---|
1887 | fCol->putColumn( *sharedFloatCol ) ;
|
---|
1888 | delete sharedFloatCol ;
|
---|
1889 | delete fCol ;
|
---|
1890 | }
|
---|
1891 | if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
|
---|
1892 | fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
|
---|
1893 | sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
|
---|
1894 | fCol->putColumn( *sharedFloatCol ) ;
|
---|
1895 | delete sharedFloatCol ;
|
---|
1896 | delete fCol ;
|
---|
1897 | }
|
---|
1898 | if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
|
---|
1899 | fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
|
---|
1900 | sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
|
---|
1901 | fCol->putColumn( *sharedFloatCol ) ;
|
---|
1902 | delete sharedFloatCol ;
|
---|
1903 | delete fCol ;
|
---|
1904 | }
|
---|
1905 | if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
|
---|
1906 | fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
|
---|
1907 | sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
|
---|
1908 | fCol->putColumn( *sharedFloatCol ) ;
|
---|
1909 | delete sharedFloatCol ;
|
---|
1910 | delete fCol ;
|
---|
1911 | }
|
---|
1912 | ScalarColumn<uInt> idCol( wtab, "ID" ) ;
|
---|
1913 | for ( uInt irow = 0 ; irow < wnrow ; irow++ )
|
---|
1914 | idCol.put( irow, irow ) ;
|
---|
1915 |
|
---|
1916 | ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
|
---|
1917 | ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
|
---|
1918 | String tUnit = tqCol.getUnits() ;
|
---|
1919 | Vector<Double> mwTime = tCol.getColumn() ;
|
---|
1920 | if ( tUnit == "d" )
|
---|
1921 | mwTime *= 86400.0 ;
|
---|
1922 | tqCol.attach( mWeatherSel, "INTERVAL" ) ;
|
---|
1923 | tCol.attach( mWeatherSel, "INTERVAL" ) ;
|
---|
1924 | String iUnit = tqCol.getUnits() ;
|
---|
1925 | Vector<Double> mwInterval = tCol.getColumn() ;
|
---|
1926 | if ( iUnit == "d" )
|
---|
1927 | mwInterval *= 86400.0 ;
|
---|
1928 |
|
---|
1929 | if ( stationId > 0 ) {
|
---|
1930 | ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
|
---|
1931 | Vector<Int> stationVec = stationCol.getColumn() ;
|
---|
1932 | uInt wsnrow = ntrue( stationVec == stationId ) ;
|
---|
1933 | mwTime_.resize( wsnrow ) ;
|
---|
1934 | mwInterval_.resize( wsnrow ) ;
|
---|
1935 | mwIndex_.resize( wsnrow ) ;
|
---|
1936 | uInt wsidx = 0 ;
|
---|
1937 | for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
|
---|
1938 | if ( stationId == stationVec[irow] ) {
|
---|
1939 | mwTime_[wsidx] = mwTime[irow] ;
|
---|
1940 | mwInterval_[wsidx] = mwInterval[irow] ;
|
---|
1941 | mwIndex_[wsidx] = irow ;
|
---|
1942 | wsidx++ ;
|
---|
1943 | }
|
---|
1944 | }
|
---|
1945 | }
|
---|
1946 | else {
|
---|
1947 | mwTime_ = mwTime ;
|
---|
1948 | mwInterval_ = mwInterval ;
|
---|
1949 | mwIndex_.resize( mwTime_.size() ) ;
|
---|
1950 | indgen( mwIndex_ ) ;
|
---|
1951 | }
|
---|
1952 | //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
|
---|
1953 | //double endSec = mathutil::gettimeofday_sec() ;
|
---|
1954 | //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
|
---|
1955 | }
|
---|
1956 |
|
---|
1957 | void MSFiller::fillFocus()
|
---|
1958 | {
|
---|
1959 | //double startSec = mathutil::gettimeofday_sec() ;
|
---|
1960 | //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
|
---|
1961 | // tentative
|
---|
1962 | table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
|
---|
1963 | //double endSec = mathutil::gettimeofday_sec() ;
|
---|
1964 | //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
|
---|
1965 | }
|
---|
1966 |
|
---|
1967 | void MSFiller::fillTcal()
|
---|
1968 | {
|
---|
1969 | //double startSec = mathutil::gettimeofday_sec() ;
|
---|
1970 | //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
|
---|
1971 |
|
---|
1972 | if ( !isSysCal_ ) {
|
---|
1973 | // add dummy row
|
---|
1974 | os_ << "No SYSCAL rows" << LogIO::POST ;
|
---|
1975 | table_->tcal().table().addRow(1,True) ;
|
---|
1976 | Vector<Float> defaultTcal( 1, 1.0 ) ;
|
---|
1977 | ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
|
---|
1978 | tcalCol.put( 0, defaultTcal ) ;
|
---|
1979 | return ;
|
---|
1980 | }
|
---|
1981 |
|
---|
1982 | if ( colTcal_ == "NONE" ) {
|
---|
1983 | // add dummy row
|
---|
1984 | os_ << "No TCAL column" << LogIO::POST ;
|
---|
1985 | table_->tcal().table().addRow(1,True) ;
|
---|
1986 | Vector<Float> defaultTcal( 1, 1.0 ) ;
|
---|
1987 | ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
|
---|
1988 | tcalCol.put( 0, defaultTcal ) ;
|
---|
1989 | return ;
|
---|
1990 | }
|
---|
1991 |
|
---|
1992 | Table &sctab = mstable_.sysCal() ;
|
---|
1993 | if ( sctab.nrow() == 0 ) {
|
---|
1994 | os_ << "No SYSCAL rows" << LogIO::POST ;
|
---|
1995 | return ;
|
---|
1996 | }
|
---|
1997 | ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
|
---|
1998 | Vector<Int> ant = antCol.getColumn() ;
|
---|
1999 | if ( allNE( ant, antenna_ ) ) {
|
---|
2000 | os_ << "No SYSCAL rows" << LogIO::POST ;
|
---|
2001 | return ;
|
---|
2002 | }
|
---|
2003 | ROTableColumn tcalCol( sctab, colTcal_ ) ;
|
---|
2004 | Bool notDefined = False ;
|
---|
2005 | for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
|
---|
2006 | if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
|
---|
2007 | notDefined = True ;
|
---|
2008 | break ;
|
---|
2009 | }
|
---|
2010 | }
|
---|
2011 | if ( notDefined ) {
|
---|
2012 | os_ << "No TCAL value" << LogIO::POST ;
|
---|
2013 | table_->tcal().table().addRow(1,True) ;
|
---|
2014 | Vector<Float> defaultTcal( 1, 1.0 ) ;
|
---|
2015 | ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
|
---|
2016 | tcalCol.put( 0, defaultTcal ) ;
|
---|
2017 | return ;
|
---|
2018 | }
|
---|
2019 |
|
---|
2020 | static const char *cols[] = {
|
---|
2021 | "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
|
---|
2022 | NULL
|
---|
2023 | };
|
---|
2024 | static const TypeManagerImpl<Int> tmInt;
|
---|
2025 | static const TypeManagerImpl<Double> tmDouble;
|
---|
2026 | static const TypeManager *const tms[] = {
|
---|
2027 | &tmInt, &tmInt, &tmInt, &tmDouble, NULL
|
---|
2028 | };
|
---|
2029 | Table tab = table_->tcal().table() ;
|
---|
2030 | TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
|
---|
2031 | visitor.setTcalColumn( colTcal_ ) ;
|
---|
2032 |
|
---|
2033 | traverseTable(sctab, cols, tms, &visitor);
|
---|
2034 |
|
---|
2035 | //tcalrec_.print( std::cout ) ;
|
---|
2036 | //double endSec = mathutil::gettimeofday_sec() ;
|
---|
2037 | //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
|
---|
2038 | }
|
---|
2039 |
|
---|
2040 | string MSFiller::getFrame()
|
---|
2041 | {
|
---|
2042 | MFrequency::Types frame = MFrequency::DEFAULT ;
|
---|
2043 | ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
|
---|
2044 | ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
|
---|
2045 | uInt nrow = numChanCol.nrow() ;
|
---|
2046 | Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
|
---|
2047 | uInt nref = 0 ;
|
---|
2048 | for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
|
---|
2049 | if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
|
---|
2050 | measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
|
---|
2051 | nref++ ;
|
---|
2052 | }
|
---|
2053 | }
|
---|
2054 | if ( nref > 0 )
|
---|
2055 | frame = (MFrequency::Types)measFreqRef[0] ;
|
---|
2056 |
|
---|
2057 | return MFrequency::showType( frame ) ;
|
---|
2058 | }
|
---|
2059 |
|
---|
2060 | void MSFiller::initHeader( STHeader &header )
|
---|
2061 | {
|
---|
2062 | header.nchan = 0 ;
|
---|
2063 | header.npol = 0 ;
|
---|
2064 | header.nif = 0 ;
|
---|
2065 | header.nbeam = 0 ;
|
---|
2066 | header.observer = "" ;
|
---|
2067 | header.project = "" ;
|
---|
2068 | header.obstype = "" ;
|
---|
2069 | header.antennaname = "" ;
|
---|
2070 | header.antennaposition.resize( 3 ) ;
|
---|
2071 | header.equinox = 0.0 ;
|
---|
2072 | header.freqref = "" ;
|
---|
2073 | header.reffreq = -1.0 ;
|
---|
2074 | header.bandwidth = 0.0 ;
|
---|
2075 | header.utc = 0.0 ;
|
---|
2076 | header.fluxunit = "" ;
|
---|
2077 | header.epoch = "" ;
|
---|
2078 | header.poltype = "" ;
|
---|
2079 | }
|
---|
2080 |
|
---|
2081 | };
|
---|