source: trunk/src/MSFiller.cpp@ 2520

Last change on this file since 2520 was 2303, checked in by Takeshi Nakazato, 13 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...

Properly handle the case when unexisting antenna is specified.


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 // 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 ; }
711private:
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
1319class BaseTcalVisitor: public TableVisitor {
1320 uInt lastRecordNo ;
1321 Int lastAntennaId ;
1322 Int lastFeedId ;
1323 Int lastSpwId ;
1324 Double lastTime ;
1325protected:
1326 const Table &table;
1327 uInt count;
1328public:
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
1435class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
1436public:
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 }
1526private:
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
1547MSFiller::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
1569MSFiller::~MSFiller()
1570{
1571 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
1572}
1573
1574bool 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 else {
1609 delete tmpMS ;
1610 //throw( AipsError( "Antenna " + antennaStr_ + " doesn't exist." ) ) ;
1611 os_ << LogIO::SEVERE << "Antenna " << antennaStr_ << " doesn't exist." << LogIO::POST ;
1612 return False ;
1613 }
1614 }
1615
1616 os_ << "Parsing MS options" << endl ;
1617 os_ << " getPt = " << getPt_ << endl ;
1618 os_ << " antenna = " << antenna_ << endl ;
1619 os_ << " antennaStr = " << antennaStr_ << LogIO::POST ;
1620
1621 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
1622 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
1623
1624 delete tmpMS ;
1625
1626 // check which data column exists
1627 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
1628 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
1629
1630 //double endSec = mathutil::gettimeofday_sec() ;
1631 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1632 return true ;
1633}
1634
1635void MSFiller::fill()
1636{
1637 //double startSec = mathutil::gettimeofday_sec() ;
1638 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
1639
1640 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
1641
1642 // Initialize header
1643 STHeader sdh ;
1644 initHeader( sdh ) ;
1645 table_->setHeader( sdh ) ;
1646
1647 // check if optional table exists
1648 const TableRecord &msrec = mstable_.keywordSet() ;
1649 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
1650 if ( isDoppler_ )
1651 if ( mstable_.doppler().nrow() == 0 )
1652 isDoppler_ = False ;
1653 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
1654 if ( isFlagCmd_ )
1655 if ( mstable_.flagCmd().nrow() == 0 )
1656 isFlagCmd_ = False ;
1657 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
1658 if ( isFreqOffset_ )
1659 if ( mstable_.freqOffset().nrow() == 0 )
1660 isFreqOffset_ = False ;
1661 isHistory_ = msrec.isDefined( "HISTORY" ) ;
1662 if ( isHistory_ )
1663 if ( mstable_.history().nrow() == 0 )
1664 isHistory_ = False ;
1665 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
1666 if ( isProcessor_ )
1667 if ( mstable_.processor().nrow() == 0 )
1668 isProcessor_ = False ;
1669 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
1670 if ( isSysCal_ )
1671 if ( mstable_.sysCal().nrow() == 0 )
1672 isSysCal_ = False ;
1673 isWeather_ = msrec.isDefined( "WEATHER" ) ;
1674 if ( isWeather_ )
1675 if ( mstable_.weather().nrow() == 0 )
1676 isWeather_ = False ;
1677
1678 // column name for Tsys and Tcal
1679 if ( isSysCal_ ) {
1680 const MSSysCal &caltab = mstable_.sysCal() ;
1681 if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
1682 colTcal_ = "TCAL" ;
1683 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1684 colTcal_ = "NONE" ;
1685 }
1686 if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
1687 colTsys_ = "TSYS" ;
1688 if ( !caltab.tableDesc().isColumn( colTcal_ ) )
1689 colTsys_ = "NONE" ;
1690 }
1691 }
1692 else {
1693 colTcal_ = "NONE" ;
1694 colTsys_ = "NONE" ;
1695 }
1696
1697 // Access to MS subtables
1698 //MSField &fieldtab = mstable_.field() ;
1699 //MSPolarization &poltab = mstable_.polarization() ;
1700 //MSDataDescription &ddtab = mstable_.dataDescription() ;
1701 //MSObservation &obstab = mstable_.observation() ;
1702 //MSSource &srctab = mstable_.source() ;
1703 //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
1704 //MSSysCal &caltab = mstable_.sysCal() ;
1705 MSPointing &pointtab = mstable_.pointing() ;
1706 //MSState &stattab = mstable_.state() ;
1707 //MSAntenna &anttab = mstable_.antenna() ;
1708
1709 // SUBTABLES: FREQUENCIES
1710 //string freqFrame = getFrame() ;
1711 string freqFrame = "LSRK" ;
1712 table_->frequencies().setFrame( freqFrame ) ;
1713 table_->frequencies().setFrame( freqFrame, True ) ;
1714
1715 // SUBTABLES: WEATHER
1716 fillWeather() ;
1717
1718 // SUBTABLES: FOCUS
1719 fillFocus() ;
1720
1721 // SUBTABLES: TCAL
1722 fillTcal() ;
1723
1724 // SUBTABLES: FIT
1725 //fillFit() ;
1726
1727 // SUBTABLES: HISTORY
1728 //fillHistory() ;
1729
1730 /***
1731 * Start iteration using TableVisitor
1732 ***/
1733 Table stab = table_->table() ;
1734 {
1735 static const char *cols[] = {
1736 "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
1737 "STATE_ID", "TIME",
1738 NULL
1739 };
1740 static const TypeManagerImpl<Int> tmInt;
1741 static const TypeManagerImpl<Double> tmDouble;
1742 static const TypeManager *const tms[] = {
1743 &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
1744 };
1745 //double t0 = mathutil::gettimeofday_sec() ;
1746 MSFillerVisitor myVisitor(mstable_, *table_ );
1747 //double t1 = mathutil::gettimeofday_sec() ;
1748 //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
1749 myVisitor.setAntenna( antenna_ ) ;
1750 //myVisitor.setHeader( sdh ) ;
1751 if ( getPt_ ) {
1752 Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
1753 myVisitor.setPointingTable( ptsel ) ;
1754 }
1755 if ( isWeather_ )
1756 myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ;
1757 if ( isSysCal_ )
1758 myVisitor.setSysCalRecord( tcalrec_ ) ;
1759
1760 //double t2 = mathutil::gettimeofday_sec() ;
1761 traverseTable(mstable_, cols, tms, &myVisitor);
1762 //double t3 = mathutil::gettimeofday_sec() ;
1763 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
1764
1765 sdh = myVisitor.getHeader() ;
1766 }
1767 /***
1768 * End iteration using TableVisitor
1769 ***/
1770
1771 // set header
1772 //sdh = myVisitor.getHeader() ;
1773 //table_->setHeader( sdh ) ;
1774
1775 // save path to POINTING table
1776 // 2011/07/06 TN
1777 // Path to POINTING table in original MS will not be written
1778 // if getPt_ is True
1779 Path datapath( tablename_ ) ;
1780 if ( !getPt_ ) {
1781 String pTabName = datapath.absoluteName() + "/POINTING" ;
1782 stab.rwKeywordSet().define( "POINTING", pTabName ) ;
1783 }
1784
1785 // for GBT
1786 if ( sdh.antennaname.contains( "GBT" ) ) {
1787 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
1788 stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
1789 }
1790
1791 // for MS created from ASDM
1792 const TableRecord &msKeys = mstable_.keywordSet() ;
1793 uInt nfields = msKeys.nfields() ;
1794 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
1795 String name = msKeys.name( ifield ) ;
1796 //os_ << "name = " << name << LogIO::POST ;
1797 if ( name.find( "ASDM" ) != String::npos ) {
1798 String asdmpath = msKeys.asTable( ifield ).tableName() ;
1799 os_ << "ASDM table: " << asdmpath << LogIO::POST ;
1800 stab.rwKeywordSet().define( name, asdmpath ) ;
1801 }
1802 }
1803
1804 //double endSec = mathutil::gettimeofday_sec() ;
1805 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1806}
1807
1808void MSFiller::close()
1809{
1810 //tablesel_.closeSubTables() ;
1811 mstable_.closeSubTables() ;
1812 //tablesel_.unlock() ;
1813 mstable_.unlock() ;
1814}
1815
1816void MSFiller::fillWeather()
1817{
1818 //double startSec = mathutil::gettimeofday_sec() ;
1819 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1820
1821 if ( !isWeather_ ) {
1822 // add dummy row
1823 table_->weather().table().addRow(1,True) ;
1824 return ;
1825 }
1826
1827 Table mWeather = mstable_.weather() ;
1828 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1829 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1830 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1831 if ( mWeatherSel.nrow() == 0 ) {
1832 os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1833 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1834 if ( mWeatherSel.nrow() == 0 ) {
1835 os_ << "No rows in WEATHER table" << LogIO::POST ;
1836 }
1837 }
1838 uInt wnrow = mWeatherSel.nrow() ;
1839 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1840
1841 if ( wnrow == 0 )
1842 return ;
1843
1844 Table wtab = table_->weather().table() ;
1845 wtab.addRow( wnrow ) ;
1846
1847 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
1848 Int stationId = -1 ;
1849 if ( stationInfoExists ) {
1850 // determine which station is closer
1851 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1852 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
1853 Vector<Int> stationIds = stationCol.getColumn() ;
1854 Vector<Int> stationIdList( 0 ) ;
1855 Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
1856 uInt numStation = 0 ;
1857 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
1858 if ( !anyEQ( stationIdList, stationIds[i] ) ) {
1859 numStation++ ;
1860 stationIdList.resize( numStation, True ) ;
1861 stationIdList[numStation-1] = stationIds[i] ;
1862 stationPosList.resize( numStation, 3, True ) ;
1863 stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
1864 }
1865 }
1866 //os_ << "staionIdList = " << stationIdList << endl ;
1867 Table mAntenna = mstable_.antenna() ;
1868 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
1869 Vector<Double> antpos = antposCol( antenna_ ) ;
1870 Double minDiff = -1.0 ;
1871 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
1872 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
1873 if ( minDiff < 0.0 || minDiff > diff ) {
1874 minDiff = diff ;
1875 stationId = stationIdList[i] ;
1876 }
1877 }
1878 }
1879 //os_ << "stationId = " << stationId << endl ;
1880
1881 ScalarColumn<Float> *fCol ;
1882 ROScalarColumn<Float> *sharedFloatCol ;
1883 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1884 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1885 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1886 fCol->putColumn( *sharedFloatCol ) ;
1887 delete sharedFloatCol ;
1888 delete fCol ;
1889 }
1890 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1891 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1892 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1893 fCol->putColumn( *sharedFloatCol ) ;
1894 delete sharedFloatCol ;
1895 delete fCol ;
1896 }
1897 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1898 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1899 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1900 fCol->putColumn( *sharedFloatCol ) ;
1901 delete sharedFloatCol ;
1902 delete fCol ;
1903 }
1904 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1905 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1906 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1907 fCol->putColumn( *sharedFloatCol ) ;
1908 delete sharedFloatCol ;
1909 delete fCol ;
1910 }
1911 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1912 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1913 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1914 fCol->putColumn( *sharedFloatCol ) ;
1915 delete sharedFloatCol ;
1916 delete fCol ;
1917 }
1918 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1919 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1920 idCol.put( irow, irow ) ;
1921
1922 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1923 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1924 String tUnit = tqCol.getUnits() ;
1925 Vector<Double> mwTime = tCol.getColumn() ;
1926 if ( tUnit == "d" )
1927 mwTime *= 86400.0 ;
1928 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1929 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1930 String iUnit = tqCol.getUnits() ;
1931 Vector<Double> mwInterval = tCol.getColumn() ;
1932 if ( iUnit == "d" )
1933 mwInterval *= 86400.0 ;
1934
1935 if ( stationId > 0 ) {
1936 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
1937 Vector<Int> stationVec = stationCol.getColumn() ;
1938 uInt wsnrow = ntrue( stationVec == stationId ) ;
1939 mwTime_.resize( wsnrow ) ;
1940 mwInterval_.resize( wsnrow ) ;
1941 mwIndex_.resize( wsnrow ) ;
1942 uInt wsidx = 0 ;
1943 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
1944 if ( stationId == stationVec[irow] ) {
1945 mwTime_[wsidx] = mwTime[irow] ;
1946 mwInterval_[wsidx] = mwInterval[irow] ;
1947 mwIndex_[wsidx] = irow ;
1948 wsidx++ ;
1949 }
1950 }
1951 }
1952 else {
1953 mwTime_ = mwTime ;
1954 mwInterval_ = mwInterval ;
1955 mwIndex_.resize( mwTime_.size() ) ;
1956 indgen( mwIndex_ ) ;
1957 }
1958 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1959 //double endSec = mathutil::gettimeofday_sec() ;
1960 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1961}
1962
1963void MSFiller::fillFocus()
1964{
1965 //double startSec = mathutil::gettimeofday_sec() ;
1966 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
1967 // tentative
1968 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
1969 //double endSec = mathutil::gettimeofday_sec() ;
1970 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1971}
1972
1973void MSFiller::fillTcal()
1974{
1975 //double startSec = mathutil::gettimeofday_sec() ;
1976 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
1977
1978 if ( !isSysCal_ ) {
1979 // add dummy row
1980 os_ << "No SYSCAL rows" << LogIO::POST ;
1981 table_->tcal().table().addRow(1,True) ;
1982 Vector<Float> defaultTcal( 1, 1.0 ) ;
1983 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1984 tcalCol.put( 0, defaultTcal ) ;
1985 return ;
1986 }
1987
1988 if ( colTcal_ == "NONE" ) {
1989 // add dummy row
1990 os_ << "No TCAL column" << 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 Table &sctab = mstable_.sysCal() ;
1999 if ( sctab.nrow() == 0 ) {
2000 os_ << "No SYSCAL rows" << LogIO::POST ;
2001 return ;
2002 }
2003 ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
2004 Vector<Int> ant = antCol.getColumn() ;
2005 if ( allNE( ant, antenna_ ) ) {
2006 os_ << "No SYSCAL rows" << LogIO::POST ;
2007 return ;
2008 }
2009 ROTableColumn tcalCol( sctab, colTcal_ ) ;
2010 Bool notDefined = False ;
2011 for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
2012 if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
2013 notDefined = True ;
2014 break ;
2015 }
2016 }
2017 if ( notDefined ) {
2018 os_ << "No TCAL value" << LogIO::POST ;
2019 table_->tcal().table().addRow(1,True) ;
2020 Vector<Float> defaultTcal( 1, 1.0 ) ;
2021 ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
2022 tcalCol.put( 0, defaultTcal ) ;
2023 return ;
2024 }
2025
2026 static const char *cols[] = {
2027 "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
2028 NULL
2029 };
2030 static const TypeManagerImpl<Int> tmInt;
2031 static const TypeManagerImpl<Double> tmDouble;
2032 static const TypeManager *const tms[] = {
2033 &tmInt, &tmInt, &tmInt, &tmDouble, NULL
2034 };
2035 Table tab = table_->tcal().table() ;
2036 TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
2037 visitor.setTcalColumn( colTcal_ ) ;
2038
2039 traverseTable(sctab, cols, tms, &visitor);
2040
2041 //tcalrec_.print( std::cout ) ;
2042 //double endSec = mathutil::gettimeofday_sec() ;
2043 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2044}
2045
2046string MSFiller::getFrame()
2047{
2048 MFrequency::Types frame = MFrequency::DEFAULT ;
2049 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
2050 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
2051 uInt nrow = numChanCol.nrow() ;
2052 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
2053 uInt nref = 0 ;
2054 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
2055 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
2056 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
2057 nref++ ;
2058 }
2059 }
2060 if ( nref > 0 )
2061 frame = (MFrequency::Types)measFreqRef[0] ;
2062
2063 return MFrequency::showType( frame ) ;
2064}
2065
2066void MSFiller::initHeader( STHeader &header )
2067{
2068 header.nchan = 0 ;
2069 header.npol = 0 ;
2070 header.nif = 0 ;
2071 header.nbeam = 0 ;
2072 header.observer = "" ;
2073 header.project = "" ;
2074 header.obstype = "" ;
2075 header.antennaname = "" ;
2076 header.antennaposition.resize( 3 ) ;
2077 header.equinox = 0.0 ;
2078 header.freqref = "" ;
2079 header.reffreq = -1.0 ;
2080 header.bandwidth = 0.0 ;
2081 header.utc = 0.0 ;
2082 header.fluxunit = "" ;
2083 header.epoch = "" ;
2084 header.poltype = "" ;
2085}
2086
2087};
Note: See TracBrowser for help on using the repository browser.