source: trunk/src/MSFiller.cpp@ 2737

Last change on this file since 2737 was 2710, checked in by Takeshi Nakazato, 12 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...

Reference epoch for frequency frame conversion is taken from header
which is set as start time of the observation (OBSERVATION/TIME_RANGE[0]
for MS, ExecBlock/startTime for ASDM).


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