source: trunk/src/MSWriter.cpp@ 2859

Last change on this file since 2859 was 2842, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5545

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave[sdsave_flagging]

Put in Release Notes: No

Module(s): sd

Description: Describe your changes here...

Changed a behavior of MSWriter according to CAS-5545.

  • All the FLAG_ROW's are always set to False.
  • If FLAGROW is non-zero in a certain Scantable row, all the FLAG values in the corresponding location in MS/MAIN (FLAG in the corresponding polarization component) have to be set to True.


File size: 79.1 KB
Line 
1//
2// C++ Interface: MSWriter
3//
4// Description:
5//
6// This class is specific writer for MS format
7//
8// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
9//
10// Copyright: See COPYING file that comes with this distribution
11//
12//
13#include <assert.h>
14
15#include <set>
16
17#include <casa/OS/File.h>
18#include <casa/OS/RegularFile.h>
19#include <casa/OS/Directory.h>
20#include <casa/OS/SymLink.h>
21#include <casa/BasicSL/String.h>
22#include <casa/Arrays/Cube.h>
23#include <casa/Containers/RecordField.h>
24
25#include <tables/Tables/ExprNode.h>
26#include <tables/Tables/TableDesc.h>
27#include <tables/Tables/SetupNewTab.h>
28#include <tables/Tables/TableIter.h>
29#include <tables/Tables/RefRows.h>
30#include <tables/Tables/TableRow.h>
31
32#include <ms/MeasurementSets/MeasurementSet.h>
33#include <ms/MeasurementSets/MSColumns.h>
34#include <ms/MeasurementSets/MSPolIndex.h>
35#include <ms/MeasurementSets/MSDataDescIndex.h>
36#include <ms/MeasurementSets/MSSourceIndex.h>
37
38#include "MSWriter.h"
39#include "STHeader.h"
40#include "STFrequencies.h"
41#include "STMolecules.h"
42#include "STTcal.h"
43#include "MathUtils.h"
44#include "TableTraverse.h"
45
46using namespace casa ;
47using namespace std ;
48
49namespace asap {
50
51class CorrTypeHandler {
52public:
53 CorrTypeHandler()
54 {}
55 virtual ~CorrTypeHandler() {}
56 virtual Vector<Stokes::StokesTypes> corrType() = 0 ;
57 virtual void reset()
58 {
59 npol = 0 ;
60 }
61 void append( uInt polno )
62 {
63 polnos[npol] = polno ;
64 npol++ ;
65 }
66 uInt nPol() { return npol ; }
67protected:
68 Vector<Stokes::StokesTypes> polmap ;
69 uInt polnos[4] ;
70 uInt npol ;
71};
72
73class LinearHandler : public CorrTypeHandler {
74public:
75 LinearHandler()
76 : CorrTypeHandler()
77 {
78 initMap() ;
79 }
80 virtual ~LinearHandler() {}
81 virtual Vector<Stokes::StokesTypes> corrType()
82 {
83 Vector<Stokes::StokesTypes> ret( npol, Stokes::Undefined ) ;
84 if ( npol < 4 ) {
85 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
86 ret[ipol] = polmap[polnos[ipol]] ;
87 }
88 else if ( npol == 4 ) {
89 ret[0] = polmap[0] ;
90 ret[1] = polmap[2] ;
91 ret[2] = polmap[3] ;
92 ret[3] = polmap[1] ;
93 }
94 else {
95 throw( AipsError("npol > 4") ) ;
96 }
97 return ret ;
98 }
99protected:
100 void initMap()
101 {
102 polmap.resize( 4 ) ;
103 polmap[0] = Stokes::XX ;
104 polmap[1] = Stokes::YY ;
105 polmap[2] = Stokes::XY ;
106 polmap[3] = Stokes::YX ;
107 }
108};
109class CircularHandler : public CorrTypeHandler {
110public:
111 CircularHandler()
112 : CorrTypeHandler()
113 {
114 initMap() ;
115 }
116 virtual ~CircularHandler() {}
117 virtual Vector<Stokes::StokesTypes> corrType()
118 {
119 Vector<Stokes::StokesTypes> ret( npol, Stokes::Undefined ) ;
120 if ( npol < 4 ) {
121 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
122 ret[ipol] = polmap[polnos[ipol]] ;
123 }
124 else if ( npol == 4 ) {
125 ret[0] = polmap[0] ;
126 ret[1] = polmap[2] ;
127 ret[2] = polmap[3] ;
128 ret[3] = polmap[1] ;
129 }
130 else {
131 throw( AipsError("npol > 4") ) ;
132 }
133 return ret ;
134 }
135private:
136 void initMap()
137 {
138 polmap.resize( 4 ) ;
139 polmap[0] = Stokes::RR ;
140 polmap[1] = Stokes::LL ;
141 polmap[2] = Stokes::RL ;
142 polmap[3] = Stokes::LR ;
143 }
144};
145class StokesHandler : public CorrTypeHandler {
146public:
147 StokesHandler()
148 : CorrTypeHandler()
149 {
150 initMap() ;
151 }
152 virtual ~StokesHandler() {}
153 virtual Vector<Stokes::StokesTypes> corrType()
154 {
155 Vector<Stokes::StokesTypes> ret( npol, Stokes::Undefined ) ;
156 if ( npol <= 4 ) {
157 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
158 ret[ipol] = polmap[polnos[ipol]] ;
159 }
160 else {
161 throw( AipsError("npol > 4") ) ;
162 }
163 return ret ;
164 }
165private:
166 void initMap()
167 {
168 polmap.resize( 4 ) ;
169 polmap[0] = Stokes::I ;
170 polmap[1] = Stokes::Q ;
171 polmap[2] = Stokes::U ;
172 polmap[3] = Stokes::V ;
173 }
174};
175class LinPolHandler : public CorrTypeHandler {
176public:
177 LinPolHandler()
178 : CorrTypeHandler()
179 {
180 initMap() ;
181 }
182 virtual ~LinPolHandler() {}
183 virtual Vector<Stokes::StokesTypes> corrType()
184 {
185 Vector<Stokes::StokesTypes> ret( npol, Stokes::Undefined ) ;
186 if ( npol <= 2 ) {
187 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
188 ret[ipol] = polmap[polnos[ipol]] ;
189 }
190 else {
191 throw( AipsError("npol > 4") ) ;
192 }
193 return ret ;
194 }
195private:
196 void initMap()
197 {
198 polmap.resize( 2 ) ;
199 polmap[0] = Stokes::Plinear ;
200 polmap[1] = Stokes::Pangle ;
201 }
202};
203
204class DataHolder {
205public:
206 DataHolder( TableRow &tableRow, String polType )
207 : row( tableRow )
208 {
209 nchan = 0 ;
210 npol = 0 ;
211 makeCorrTypeHandler( polType ) ;
212 attach() ;
213 flagRow.resize( 4 ) ;
214 reset() ;
215 sigmaTemplate.resize( 4 ) ;
216 sigmaTemplate = 1.0 ;
217 }
218 virtual ~DataHolder() {}
219 void post()
220 {
221 postData() ;
222 postFlag() ;
223 postFlagRow() ;
224 postAuxiliary() ;
225 }
226 virtual void reset()
227 {
228 corr->reset() ;
229 flagRow = False ;
230 npol = 0 ;
231 for ( uInt i = 0 ; i < 4 ; i++ )
232 isFilled[i] = False ;
233 }
234 virtual void accumulate( uInt id, Vector<Float> &sp, Vector<Bool> &fl, Bool &flr )
235 {
236 accumulateCorrType( id ) ;
237 accumulateData( id, sp ) ;
238 accumulateFlag( id, fl ) ;
239 accumulateFlagRow( id, flr ) ;
240 isFilled[id] = True;
241 }
242 uInt nPol() { return npol ; }
243 uInt nChan() { return nchan ; }
244 Vector<Int> corrTypeInt()
245 {
246 Vector<Int> v( npol ) ;
247 convertArray( v, corr->corrType() ) ;
248 return v ;
249 }
250 Vector<Stokes::StokesTypes> corrType() { return corr->corrType() ; }
251 void setNchan( uInt num )
252 {
253 nchan = num ;
254 resize() ;
255 }
256protected:
257 void postAuxiliary()
258 {
259 Vector<Float> w = sigmaTemplate( IPosition(1,0), IPosition(1,npol-1) ) ;
260 sigmaRF.define( w ) ;
261 weightRF.define( w ) ;
262 Cube<Bool> c( npol, nchan, 1, False ) ;
263 flagCategoryRF.define( c ) ;
264 }
265 inline void accumulateCorrType( uInt &id )
266 {
267 corr->append( id ) ;
268 npol = corr->nPol() ;
269 }
270 inline void accumulateFlagRow( uInt &id, Bool &flr )
271 {
272 flagRow[id] = flr ;
273 }
274 void postFlagRow()
275 {
276 // CAS-5545 FLAG_ROW must always be set False
277 //*flagRowRF = anyEQ( flagRow, True ) ;
278 *flagRowRF = False;
279 }
280 inline void accumulateFlag( uInt &id, Vector<Bool> &fl )
281 {
282 flag.row( id ) = fl ;
283 }
284 virtual void postFlag() = 0 ;
285 inline void accumulateData( uInt &id, Vector<Float> &sp )
286 {
287 data.row( id ) = sp ;
288 }
289 uInt filledIndex()
290 {
291 uInt idx = 0;
292 while( !isFilled[idx] && idx < 4 ) ++idx;
293 return idx;
294 }
295 virtual void postData() = 0 ;
296 TableRow &row ;
297 uInt nchan ;
298 uInt npol ;
299 CountedPtr<CorrTypeHandler> corr;
300 RecordFieldPtr< Vector<Float> > sigmaRF ;
301 RecordFieldPtr< Vector<Float> > weightRF ;
302 RecordFieldPtr< Array<Bool> > flagRF ;
303 RecordFieldPtr<Bool> flagRowRF ;
304 RecordFieldPtr< Cube<Bool> > flagCategoryRF ;
305 Vector<Bool> flagRow ;
306 Matrix<Bool> flag ;
307 Matrix<Float> data ;
308 Vector<Float> sigmaTemplate ;
309 Bool isFilled[4] ;
310private:
311 void makeCorrTypeHandler( String &polType )
312 {
313 if ( polType == "linear" )
314 corr = new LinearHandler() ;
315 else if ( polType == "circular" )
316 corr = new CircularHandler() ;
317 else if ( polType == "stokes" )
318 corr = new StokesHandler() ;
319 else if ( polType == "linpol" )
320 corr = new LinPolHandler() ;
321 else
322 throw( AipsError("Invalid polarization type") ) ;
323 }
324 void attach()
325 {
326 TableRecord &rec = row.record() ;
327 sigmaRF.attachToRecord( rec, "SIGMA" ) ;
328 weightRF.attachToRecord( rec, "WEIGHT" ) ;
329 flagRF.attachToRecord( rec, "FLAG" ) ;
330 flagRowRF.attachToRecord( rec, "FLAG_ROW" ) ;
331 flagCategoryRF.attachToRecord( rec, "FLAG_CATEGORY" ) ;
332 }
333 void resize()
334 {
335 flag.resize( 4, nchan ) ;
336 data.resize( 4, nchan ) ;
337 }
338};
339
340class FloatDataHolder : public DataHolder {
341public:
342 FloatDataHolder( TableRow &tableRow, String polType )
343 : DataHolder( tableRow, polType )
344 {
345 attachData() ;
346 }
347 virtual ~FloatDataHolder() {}
348protected:
349 virtual void postFlag()
350 {
351 // CAS-5545 FLAG must all be set True if FLAG_ROW is True
352 for (uInt ipol = 0; ipol < npol; ++ipol) {
353 if (flagRow[ipol]) {
354 flag(IPosition(2, ipol, 0), IPosition(2, ipol, nchan-1)) = True;
355 }
356 }
357
358 if ( npol == 2 ) {
359 flagRF.define( flag( IPosition( 2, 0, 0 ), IPosition( 2, npol-1, nchan-1 ) ) ) ;
360 }
361 else {
362 // should be npol == 1
363 uInt idx = filledIndex() ;
364 flagRF.define( flag( IPosition( 2, idx, 0 ), IPosition( 2, idx, nchan-1 ) ) ) ;
365 }
366 }
367 virtual void postData()
368 {
369 if ( npol == 2 ) {
370 dataRF.define( data( IPosition( 2, 0, 0 ), IPosition( 2, npol-1, nchan-1 ) ) ) ;
371 }
372 else {
373 // should be npol == 1
374 uInt idx = filledIndex() ;
375 dataRF.define( data( IPosition( 2, idx, 0 ), IPosition( 2, idx, nchan-1 ) ) ) ;
376 }
377 }
378private:
379 void attachData()
380 {
381 dataRF.attachToRecord( row.record(), "FLOAT_DATA" ) ;
382 }
383 RecordFieldPtr< Matrix<Float> > dataRF;
384};
385
386class ComplexDataHolder : public DataHolder {
387public:
388 ComplexDataHolder( TableRow &tableRow, String polType )
389 : DataHolder( tableRow, polType )
390 {
391 attachData() ;
392 }
393 virtual ~ComplexDataHolder() {}
394protected:
395 virtual void postFlag()
396 {
397 // CAS-5545 FLAG must all be set True if FLAG_ROW is True
398 for (uInt ipol = 0; ipol < npol; ++ipol) {
399 if (flagRow[ipol]) {
400 flag(IPosition(2, ipol, 0), IPosition(2, ipol, nchan-1)) = True;
401 }
402 }
403
404 if ( npol == 4 ) {
405 Vector<Bool> tmp = flag.row( 3 ) ;
406 flag.row( 3 ) = flag.row( 1 ) ;
407 flag.row( 2 ) = flag.row( 2 ) || tmp ;
408 flag.row( 1 ) = flag.row( 2 ) ;
409 flagRF.define( flag ) ;
410 }
411 else if ( npol == 2 ) {
412 flagRF.define( flag( IPosition( 2, 0, 0 ), IPosition( 2, npol-1, nchan-1 ) ) ) ;
413 }
414 else {
415 // should be npol == 1
416 uInt idx = filledIndex() ;
417 flagRF.define( flag( IPosition( 2, idx, 0 ), IPosition( 2, idx, nchan-1 ) ) ) ;
418 }
419 }
420 virtual void postData()
421 {
422 Matrix<Float> tmp( 2, nchan, 0.0 ) ;
423 Matrix<Complex> v( npol, nchan ) ;
424 if ( isFilled[0] ) {
425 tmp.row( 0 ) = data.row( 0 ) ;
426 v.row( 0 ) = RealToComplex( tmp ) ;
427 }
428 if ( isFilled[1] ) {
429 tmp.row( 0 ) = data.row( 1 ) ;
430 v.row( npol-1 ) = RealToComplex( tmp ) ;
431 }
432 if ( isFilled[2] && isFilled[3] ) {
433 tmp.row( 0 ) = data.row( 2 ) ;
434 tmp.row( 1 ) = data.row( 3 ) ;
435 v.row( 1 ) = RealToComplex( tmp ) ;
436 v.row( 2 ) = conj( v.row( 1 ) ) ;
437 }
438 dataRF.define( v ) ;
439 }
440private:
441 void attachData()
442 {
443 dataRF.attachToRecord( row.record(), "DATA" ) ;
444 }
445 RecordFieldPtr< Matrix<Complex> > dataRF;
446};
447
448class BaseMSWriterVisitor: public TableVisitor {
449 const String *lastFieldName;
450 uInt lastRecordNo;
451 uInt lastBeamNo, lastScanNo, lastIfNo, lastPolNo;
452 Int lastSrcType;
453 uInt lastCycleNo;
454 Double lastTime;
455protected:
456 const Table &table;
457 uInt count;
458public:
459 BaseMSWriterVisitor(const Table &table)
460 : table(table)
461 {
462 static const String dummy;
463 lastFieldName = &dummy;
464 count = 0;
465 }
466
467 virtual void enterFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
468 }
469 virtual void leaveFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
470 }
471 virtual void enterBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
472 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
473 virtual void enterScanNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
474 virtual void leaveScanNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
475 virtual void enterIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
476 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
477 virtual void enterSrcType(const uInt /*recordNo*/, Int /*columnValue*/) { }
478 virtual void leaveSrcType(const uInt /*recordNo*/, Int /*columnValue*/) { }
479 virtual void enterCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
480 virtual void leaveCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
481 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
482 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
483 virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
484 virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
485
486 virtual Bool visitRecord(const uInt /*recordNo*/,
487 const String &/*fieldName*/,
488 const uInt /*beamNo*/,
489 const uInt /*scanNo*/,
490 const uInt /*ifNo*/,
491 const Int /*srcType*/,
492 const uInt /*cycleNo*/,
493 const Double /*time*/,
494 const uInt /*polNo*/) { return True ;}
495
496 virtual Bool visit(Bool isFirst, const uInt recordNo,
497 const uInt nCols, void const *const colValues[]) {
498 const String *fieldName = NULL;
499 uInt beamNo, scanNo, ifNo;
500 Int srcType;
501 uInt cycleNo;
502 Double time;
503 uInt polNo;
504 { // prologue
505 uInt i = 0;
506 {
507 const String *col = (const String*)colValues[i++];
508 fieldName = &col[recordNo];
509 }
510 {
511 const uInt *col = (const uInt *)colValues[i++];
512 beamNo = col[recordNo];
513 }
514 {
515 const uInt *col = (const uInt *)colValues[i++];
516 scanNo = col[recordNo];
517 }
518 {
519 const uInt *col = (const uInt *)colValues[i++];
520 ifNo = col[recordNo];
521 }
522 {
523 const Int *col = (const Int *)colValues[i++];
524 srcType = col[recordNo];
525 }
526 {
527 const uInt *col = (const uInt *)colValues[i++];
528 cycleNo = col[recordNo];
529 }
530 {
531 const Double *col = (const Double *)colValues[i++];
532 time = col[recordNo];
533 }
534 {
535 const Int *col = (const Int *)colValues[i++];
536 polNo = col[recordNo];
537 }
538 assert(nCols == i);
539 }
540
541 if (isFirst) {
542 enterFieldName(recordNo, *fieldName);
543 enterBeamNo(recordNo, beamNo);
544 enterScanNo(recordNo, scanNo);
545 enterIfNo(recordNo, ifNo);
546 enterSrcType(recordNo, srcType);
547 enterCycleNo(recordNo, cycleNo);
548 enterTime(recordNo, time);
549 enterPolNo(recordNo, polNo);
550 } else {
551 if (lastFieldName->compare(*fieldName) != 0) {
552 leavePolNo(lastRecordNo, lastPolNo);
553 leaveTime(lastRecordNo, lastTime);
554 leaveCycleNo(lastRecordNo, lastCycleNo);
555 leaveSrcType(lastRecordNo, lastSrcType);
556 leaveIfNo(lastRecordNo, lastIfNo);
557 leaveScanNo(lastRecordNo, lastScanNo);
558 leaveBeamNo(lastRecordNo, lastBeamNo);
559 leaveFieldName(lastRecordNo, *lastFieldName);
560
561 enterFieldName(recordNo, *fieldName);
562 enterBeamNo(recordNo, beamNo);
563 enterScanNo(recordNo, scanNo);
564 enterIfNo(recordNo, ifNo);
565 enterSrcType(recordNo, srcType);
566 enterCycleNo(recordNo, cycleNo);
567 enterTime(recordNo, time);
568 enterPolNo(recordNo, polNo);
569 } else if (lastBeamNo != beamNo) {
570 leavePolNo(lastRecordNo, lastPolNo);
571 leaveTime(lastRecordNo, lastTime);
572 leaveCycleNo(lastRecordNo, lastCycleNo);
573 leaveSrcType(lastRecordNo, lastSrcType);
574 leaveIfNo(lastRecordNo, lastIfNo);
575 leaveScanNo(lastRecordNo, lastScanNo);
576 leaveBeamNo(lastRecordNo, lastBeamNo);
577
578 enterBeamNo(recordNo, beamNo);
579 enterScanNo(recordNo, scanNo);
580 enterIfNo(recordNo, ifNo);
581 enterSrcType(recordNo, srcType);
582 enterCycleNo(recordNo, cycleNo);
583 enterTime(recordNo, time);
584 enterPolNo(recordNo, polNo);
585 } else if (lastScanNo != scanNo) {
586 leavePolNo(lastRecordNo, lastPolNo);
587 leaveTime(lastRecordNo, lastTime);
588 leaveCycleNo(lastRecordNo, lastCycleNo);
589 leaveSrcType(lastRecordNo, lastSrcType);
590 leaveIfNo(lastRecordNo, lastIfNo);
591 leaveScanNo(lastRecordNo, lastScanNo);
592
593 enterScanNo(recordNo, scanNo);
594 enterIfNo(recordNo, ifNo);
595 enterSrcType(recordNo, srcType);
596 enterCycleNo(recordNo, cycleNo);
597 enterTime(recordNo, time);
598 enterPolNo(recordNo, polNo);
599 } else if (lastIfNo != ifNo) {
600 leavePolNo(lastRecordNo, lastPolNo);
601 leaveTime(lastRecordNo, lastTime);
602 leaveCycleNo(lastRecordNo, lastCycleNo);
603 leaveSrcType(lastRecordNo, lastSrcType);
604 leaveIfNo(lastRecordNo, lastIfNo);
605
606 enterIfNo(recordNo, ifNo);
607 enterSrcType(recordNo, srcType);
608 enterCycleNo(recordNo, cycleNo);
609 enterTime(recordNo, time);
610 enterPolNo(recordNo, polNo);
611 } else if (lastSrcType != srcType) {
612 leavePolNo(lastRecordNo, lastPolNo);
613 leaveTime(lastRecordNo, lastTime);
614 leaveCycleNo(lastRecordNo, lastCycleNo);
615 leaveSrcType(lastRecordNo, lastSrcType);
616
617 enterSrcType(recordNo, srcType);
618 enterCycleNo(recordNo, cycleNo);
619 enterTime(recordNo, time);
620 enterPolNo(recordNo, polNo);
621 } else if (lastCycleNo != cycleNo) {
622 leavePolNo(lastRecordNo, lastPolNo);
623 leaveTime(lastRecordNo, lastTime);
624 leaveCycleNo(lastRecordNo, lastCycleNo);
625
626 enterCycleNo(recordNo, cycleNo);
627 enterTime(recordNo, time);
628 enterPolNo(recordNo, polNo);
629 } else if (lastTime != time) {
630 leavePolNo(lastRecordNo, lastPolNo);
631 leaveTime(lastRecordNo, lastTime);
632
633 enterTime(recordNo, time);
634 enterPolNo(recordNo, polNo);
635 } else if (lastPolNo != polNo) {
636 leavePolNo(lastRecordNo, lastPolNo);
637 enterPolNo(recordNo, polNo);
638 }
639 }
640 count++;
641 Bool result = visitRecord(recordNo, *fieldName, beamNo, scanNo, ifNo, srcType,
642 cycleNo, time, polNo);
643
644 { // epilogue
645 lastRecordNo = recordNo;
646
647 lastFieldName = fieldName;
648 lastBeamNo = beamNo;
649 lastScanNo = scanNo;
650 lastIfNo = ifNo;
651 lastSrcType = srcType;
652 lastCycleNo = cycleNo;
653 lastTime = time;
654 lastPolNo = polNo;
655 }
656 return result ;
657 }
658
659 virtual void finish() {
660 if (count > 0) {
661 leavePolNo(lastRecordNo, lastPolNo);
662 leaveTime(lastRecordNo, lastTime);
663 leaveCycleNo(lastRecordNo, lastCycleNo);
664 leaveSrcType(lastRecordNo, lastSrcType);
665 leaveIfNo(lastRecordNo, lastIfNo);
666 leaveScanNo(lastRecordNo, lastScanNo);
667 leaveBeamNo(lastRecordNo, lastBeamNo);
668 leaveFieldName(lastRecordNo, *lastFieldName);
669 }
670 }
671};
672
673class MSWriterVisitor: public BaseMSWriterVisitor, public MSWriterUtils {
674public:
675 MSWriterVisitor(const Table &table, Table &mstable)
676 : BaseMSWriterVisitor(table),
677 ms(mstable)
678 {
679 rowidx = 0 ;
680 fieldName = "" ;
681 defaultFieldId = 0 ;
682 spwId = -1 ;
683 subscan = 1 ;
684 ptName = "" ;
685 srcId = 0 ;
686
687 row = TableRow( ms ) ;
688
689 initPolarization() ;
690 initFrequencies() ;
691
692 //
693 // add rows to MS
694 //
695 uInt addrow = table.nrow() ;
696 ms.addRow( addrow ) ;
697
698 // attach to Scantable columns
699 spectraCol.attach( table, "SPECTRA" ) ;
700 flagtraCol.attach( table, "FLAGTRA" ) ;
701 flagRowCol.attach( table, "FLAGROW" ) ;
702 tcalIdCol.attach( table, "TCAL_ID" ) ;
703 intervalCol.attach( table, "INTERVAL" ) ;
704 directionCol.attach( table, "DIRECTION" ) ;
705 scanRateCol.attach( table, "SCANRATE" ) ;
706 timeCol.attach( table, "TIME" ) ;
707 freqIdCol.attach( table, "FREQ_ID" ) ;
708 sourceNameCol.attach( table, "SRCNAME" ) ;
709 sourceDirectionCol.attach( table, "SRCDIRECTION" ) ;
710 fieldNameCol.attach( table, "FIELDNAME" ) ;
711
712 // MS subtables
713 attachSubtables() ;
714
715 // attach to MS columns
716 attachMain() ;
717 attachPointing() ;
718 }
719
720 virtual void enterFieldName(const uInt recordNo, const String &/*columnValue*/) {
721 //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str());
722 fieldName = fieldNameCol.asString( recordNo ) ;
723 String::size_type pos = fieldName.find( "__" ) ;
724 if ( pos != String::npos ) {
725 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
726 fieldName = fieldName.substr( 0, pos ) ;
727 }
728 else {
729 fieldId = defaultFieldId ;
730 defaultFieldId++ ;
731 }
732 Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ;
733 Vector<Double> srcDir = sourceDirectionCol( recordNo ) ;
734 Vector<Double> srate = scanRateCol( recordNo ) ;
735 String srcName = sourceNameCol.asString( recordNo ) ;
736
737 addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ;
738
739 // put value
740 *fieldIdRF = fieldId ;
741 }
742 virtual void leaveFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
743 }
744 virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue) {
745 //printf("%u: BeamNo: %u\n", recordNo, columnValue);
746
747 feedId = (Int)columnValue ;
748
749 // put value
750 *feed1RF = feedId ;
751 *feed2RF = feedId ;
752 }
753 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
754 }
755 virtual void enterScanNo(const uInt /*recordNo*/, uInt columnValue) {
756 //printf("%u: ScanNo: %u\n", recordNo, columnValue);
757
758 // put value
759 // SCAN_NUMBER is 0-based in Scantable while 1-based in MS
760 *scanNumberRF = (Int)columnValue + 1 ;
761 }
762 virtual void leaveScanNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
763 subscan = 1 ;
764 }
765 virtual void enterIfNo(const uInt recordNo, uInt columnValue) {
766 //printf("%u: IfNo: %u\n", recordNo, columnValue);
767
768 spwId = (Int)columnValue ;
769 uInt freqId = freqIdCol.asuInt( recordNo ) ;
770
771 Vector<Float> sp = spectraCol( recordNo ) ;
772 uInt nchan = sp.nelements() ;
773 holder->setNchan( nchan ) ;
774
775 addSpectralWindow( spwId, freqId ) ;
776
777 addFeed( feedId, spwId ) ;
778 }
779 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
780 }
781 virtual void enterSrcType(const uInt /*recordNo*/, Int columnValue) {
782 //printf("%u: SrcType: %d\n", recordNo, columnValue);
783
784 Int stateId = addState( columnValue ) ;
785
786 // put value
787 *stateIdRF = stateId ;
788 }
789 virtual void leaveSrcType(const uInt /*recordNo*/, Int /*columnValue*/) {
790 }
791 virtual void enterCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
792 //printf("%u: CycleNo: %u\n", recordNo, columnValue);
793 }
794 virtual void leaveCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
795 }
796 virtual void enterTime(const uInt recordNo, Double columnValue) {
797 //printf("%u: Time: %f\n", recordNo, columnValue);
798
799 Double timeSec = columnValue * 86400.0 ;
800 Double interval = intervalCol.asdouble( recordNo ) ;
801
802 if ( ptName.empty() ) {
803 Vector<Double> dir = directionCol( recordNo ) ;
804 Vector<Double> rate = scanRateCol( recordNo ) ;
805 if ( anyNE( rate, 0.0 ) ) {
806 Matrix<Double> msdir( 2, 2 ) ;
807 msdir.column( 0 ) = dir ;
808 msdir.column( 1 ) = rate ;
809 addPointing( timeSec, interval, msdir ) ;
810 }
811 else {
812 Matrix<Double> msdir( 2, 1 ) ;
813 msdir.column( 0 ) = dir ;
814 addPointing( timeSec, interval, msdir ) ;
815 }
816 }
817
818 // put value
819 *timeRF = timeSec ;
820 *timeCentroidRF = timeSec ;
821 *intervalRF = interval ;
822 *exposureRF = interval ;
823 }
824 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) {
825 if ( holder->nPol() > 0 ) {
826 Int polId = addPolarization() ;
827 Int ddId = addDataDescription( polId, spwId ) ;
828
829 // put field
830 *dataDescIdRF = ddId ;
831 holder->post() ;
832
833 // commit row
834 row.put( rowidx ) ;
835 rowidx++ ;
836
837 // reset holder
838 holder->reset() ;
839 }
840 }
841 virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
842 //printf("%u: PolNo: %d\n", recordNo, columnValue);
843 }
844 virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
845 }
846
847 virtual Bool visitRecord(const uInt recordNo,
848 const String &/*fieldName*/,
849 const uInt /*beamNo*/,
850 const uInt /*scanNo*/,
851 const uInt /*ifNo*/,
852 const Int /*srcType*/,
853 const uInt /*cycleNo*/,
854 const Double /*time*/,
855 const uInt polNo) {
856 //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo,
857 // fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo);
858
859 Vector<Float> sp = spectraCol( recordNo ) ;
860 Vector<uChar> tmp = flagtraCol( recordNo ) ;
861 Vector<Bool> fl( tmp.shape() ) ;
862 convertArray( fl, tmp ) ;
863 Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ;
864 holder->accumulate( polNo, sp, fl, flr ) ;
865
866 return True ;
867 }
868
869 virtual void finish() {
870 BaseMSWriterVisitor::finish();
871 //printf("Total: %u\n", count);
872
873 // remove rows
874 if ( ms.nrow() > rowidx ) {
875 uInt numRemove = ms.nrow() - rowidx ;
876 //cout << "numRemove = " << numRemove << endl ;
877 Vector<uInt> rows( numRemove ) ;
878 indgen( rows, rowidx ) ;
879 ms.removeRow( rows ) ;
880 }
881
882 // fill empty SPECTRAL_WINDOW rows
883 infillSpectralWindow() ;
884
885 // fill empty FIELD rows
886 infillField() ;
887 }
888
889 void dataColumnName( String name )
890 {
891 if ( name == "DATA" )
892 holder = new ComplexDataHolder( row, poltype ) ;
893 else if ( name == "FLOAT_DATA" )
894 holder = new FloatDataHolder( row, poltype ) ;
895 }
896 void pointingTableName( String name ) {
897 ptName = name ;
898 }
899 void setSourceRecord( Record &r ) {
900 srcRec = r ;
901 }
902private:
903 void addField( Int &fid, String &fname, String &srcName,
904 Vector<Double> &sdir, Vector<Double> &srate,
905 Double &tSec )
906 {
907 uInt nrow = fieldtab.nrow() ;
908 while( (Int)nrow <= fid ) {
909 fieldtab.addRow( 1, True ) ;
910 nrow++ ;
911 }
912
913 Matrix<Double> dir ;
914 Int numPoly = 0 ;
915 if ( anyNE( srate, 0.0 ) ) {
916 dir.resize( 2, 2 ) ;
917 dir.column( 0 ) = sdir ;
918 dir.column( 1 ) = srate ;
919 numPoly = 1 ;
920 }
921 else {
922 dir.resize( 2, 1 ) ;
923 dir.column( 0 ) = sdir ;
924 }
925 srcId = srcRec.asInt( srcName ) ;
926
927 TableRow tr( fieldtab ) ;
928 TableRecord &r = tr.record() ;
929 putField( "NAME", r, fname ) ;
930 putField( "NUM_POLY", r, numPoly ) ;
931 putField( "TIME", r, tSec ) ;
932 putField( "SOURCE_ID", r, srcId ) ;
933 defineField( "DELAY_DIR", r, dir ) ;
934 defineField( "REFERENCE_DIR", r, dir ) ;
935 defineField( "PHASE_DIR", r, dir ) ;
936 tr.put( fid ) ;
937
938 // for POINTING table
939 *poNameRF = fname ;
940 }
941 Int addState( Int &id )
942 {
943 String obsMode ;
944 Bool isSignal ;
945 Double tnoise ;
946 Double tload ;
947 queryType( id, obsMode, isSignal, tnoise, tload ) ;
948
949 String key = obsMode+"_"+String::toString( subscan ) ;
950 Int idx = -1 ;
951 uInt nEntry = stateEntry.nelements() ;
952 for ( uInt i = 0 ; i < nEntry ; i++ ) {
953 if ( stateEntry[i] == key ) {
954 idx = i ;
955 break ;
956 }
957 }
958 if ( idx == -1 ) {
959 uInt nrow = statetab.nrow() ;
960 statetab.addRow( 1, True ) ;
961 TableRow tr( statetab ) ;
962 TableRecord &r = tr.record() ;
963 putField( "OBS_MODE", r, obsMode ) ;
964 putField( "SIG", r, isSignal ) ;
965 isSignal = !isSignal ;
966 putField( "REF", r, isSignal ) ;
967 putField( "CAL", r, tnoise ) ;
968 putField( "LOAD", r, tload ) ;
969 tr.put( nrow ) ;
970 idx = nrow ;
971
972 stateEntry.resize( nEntry+1, True ) ;
973 stateEntry[nEntry] = key ;
974 }
975 subscan++ ;
976
977 return idx ;
978 }
979 void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir )
980 {
981 uInt nrow = potab.nrow() ;
982 potab.addRow() ;
983
984 *poNumPolyRF = dir.ncolumn() - 1 ;
985 *poTimeRF = tSec ;
986 *poTimeOriginRF = tSec ;
987 *poIntervalRF = interval ;
988 poDirectionRF.define( dir ) ;
989 poTargetRF.define( dir ) ;
990 porow.put( nrow ) ;
991 }
992 Int addPolarization()
993 {
994 Int idx = -1 ;
995 Vector<Int> corrType = holder->corrTypeInt() ;
996 uInt nEntry = polEntry.size() ;
997 for ( uInt i = 0 ; i < nEntry ; i++ ) {
998 if ( polEntry[i].conform( corrType ) && allEQ( polEntry[i], corrType ) ) {
999 idx = i ;
1000 break ;
1001 }
1002 }
1003
1004 Int numCorr = holder->nPol() ;
1005 Matrix<Int> corrProduct = corrProductTemplate[numCorr].copy() ;
1006 if ( numCorr == 1 && (corrType[0] == Stokes::YY || corrType[0] == Stokes::LL ) ) {
1007 corrProduct = 1;
1008 }
1009
1010 if ( idx == -1 ) {
1011 uInt nrow = poltab.nrow() ;
1012 poltab.addRow( 1, True ) ;
1013 TableRow tr( poltab ) ;
1014 TableRecord &r = tr.record() ;
1015 putField( "NUM_CORR", r, numCorr ) ;
1016 defineField( "CORR_TYPE", r, corrType ) ;
1017 defineField( "CORR_PRODUCT", r, corrProduct ) ;
1018 tr.put( nrow ) ;
1019 idx = nrow ;
1020
1021 polEntry.resize( nEntry+1 ) ;
1022 polEntry[nEntry] = corrType ;
1023 }
1024
1025 return idx ;
1026 }
1027 Int addDataDescription( Int pid, Int sid )
1028 {
1029 Int idx = -1 ;
1030 uInt nItem = 2 ;
1031 uInt len = ddEntry.nelements() ;
1032 uInt nEntry = len / nItem ;
1033 const Int *dd_p = ddEntry.storage() ;
1034 for ( uInt i = 0 ; i < nEntry ; i++ ) {
1035 Int pol = *dd_p ;
1036 dd_p++ ;
1037 Int spw = *dd_p ;
1038 dd_p++ ;
1039 if ( pid == pol && sid == spw ) {
1040 idx = i ;
1041 break ;
1042 }
1043 }
1044
1045 if ( idx == -1 ) {
1046 uInt nrow = ddtab.nrow() ;
1047 ddtab.addRow( 1, True ) ;
1048 TableRow tr( ddtab ) ;
1049 TableRecord &r = tr.record() ;
1050 putField( "POLARIZATION_ID", r, pid ) ;
1051 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
1052 tr.put( nrow ) ;
1053 idx = nrow ;
1054
1055 ddEntry.resize( len+nItem ) ;
1056 ddEntry[len] = pid ;
1057 ddEntry[len+1] = sid ;
1058 }
1059
1060 return idx ;
1061 }
1062 void infillSpectralWindow()
1063 {
1064 ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ;
1065 Vector<Int> nchan = nchanCol.getColumn() ;
1066 TableRow tr( spwtab ) ;
1067 TableRecord &r = tr.record() ;
1068 Int mfr = freqframe ;
1069 Int oneChan = 1 ;
1070 Vector<Double> dummy( 1, 0.0 ) ;
1071 putField( "MEAS_FREQ_REF", r, mfr ) ;
1072 defineField( "CHAN_FREQ", r, dummy ) ;
1073 defineField( "CHAN_WIDTH", r, dummy ) ;
1074 defineField( "EFFECTIVE_BW", r, dummy ) ;
1075 defineField( "RESOLUTION", r, dummy ) ;
1076 putField( "NUM_CHAN", r, oneChan ) ;
1077
1078 for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
1079 if ( nchan[i] == 0 )
1080 tr.put( i ) ;
1081 }
1082 }
1083 void infillField()
1084 {
1085 ScalarColumn<Int> sourceIdCol(fieldtab, "SOURCE_ID");
1086 ArrayColumn<Double> delayDirCol(fieldtab, "DELAY_DIR");
1087 ArrayColumn<Double> phaseDirCol(fieldtab, "PHASE_DIR");
1088 ArrayColumn<Double> referenceDirCol(fieldtab, "REFERENCE_DIR");
1089 uInt nrow = fieldtab.nrow();
1090 Matrix<Double> dummy(IPosition(2, 2, 1), 0.0);
1091 for (uInt irow = 0; irow < nrow; ++irow) {
1092 if (!phaseDirCol.isDefined(irow)) {
1093 delayDirCol.put(irow, dummy);
1094 phaseDirCol.put(irow, dummy);
1095 referenceDirCol.put(irow, dummy);
1096 sourceIdCol.put(irow, -1);
1097 }
1098 }
1099 }
1100 void addSpectralWindow( Int sid, uInt fid )
1101 {
1102 if (processedFreqId.find((uInt)fid) == processedFreqId.end()
1103 || processedIFNO.find((uInt)sid) == processedIFNO.end() ) {
1104 uInt nrow = spwtab.nrow() ;
1105 while( (Int)nrow <= sid ) {
1106 spwtab.addRow( 1, True ) ;
1107 nrow++ ;
1108 }
1109 processedFreqId.insert((uInt)fid);
1110 processedIFNO.insert((uInt)sid);
1111 }
1112 else {
1113 return ;
1114 }
1115
1116
1117 Double rp = refpix[fid] ;
1118 Double rv = refval[fid] ;
1119 Double ic = increment[fid] ;
1120
1121 Int mfrInt = (Int)freqframe ;
1122 Int nchan = holder->nChan() ;
1123 Double bw = nchan * abs( ic ) ;
1124 Double reffreq = rv - rp * ic ;
1125 Int netsb = 0 ; // USB->0, LSB->1
1126 if ( ic < 0 )
1127 netsb = 1 ;
1128 Vector<Double> res( nchan, abs(ic) ) ;
1129 Vector<Double> cw( nchan, ic ) ;
1130 Vector<Double> chanf( nchan ) ;
1131 indgen( chanf, reffreq, ic ) ;
1132
1133 TableRow tr( spwtab ) ;
1134 TableRecord &r = tr.record() ;
1135 putField( "MEAS_FREQ_REF", r, mfrInt ) ;
1136 putField( "NUM_CHAN", r, nchan ) ;
1137 putField( "TOTAL_BANDWIDTH", r, bw ) ;
1138 putField( "REF_FREQUENCY", r, reffreq ) ;
1139 putField( "NET_SIDEBAND", r, netsb ) ;
1140 defineField( "RESOLUTION", r, res ) ;
1141// defineField( "CHAN_WIDTH", r, res ) ;
1142 defineField( "CHAN_WIDTH", r, cw ) ;
1143 defineField( "EFFECTIVE_BW", r, res ) ;
1144 defineField( "CHAN_FREQ", r, chanf ) ;
1145 tr.put( sid ) ;
1146 }
1147 void addFeed( Int fid, Int sid )
1148 {
1149 Int idx = -1 ;
1150 uInt nItem = 2 ;
1151 uInt len = feedEntry.nelements() ;
1152 uInt nEntry = len / nItem ;
1153 const Int *fe_p = feedEntry.storage() ;
1154 for ( uInt i = 0 ; i < nEntry ; i++ ) {
1155 Int feed = *fe_p ;
1156 fe_p++ ;
1157 Int spw = *fe_p ;
1158 fe_p++ ;
1159 if ( fid == feed && sid == spw ) {
1160 idx = i ;
1161 break ;
1162 }
1163 }
1164
1165
1166 if ( idx == -1 ) {
1167 uInt nrow = feedtab.nrow() ;
1168 feedtab.addRow( 1, True ) ;
1169 Int numReceptors = 2 ;
1170 Vector<String> polType( numReceptors ) ;
1171 Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ;
1172 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
1173 if ( poltype == "linear" ) {
1174 polType[0] = "X" ;
1175 polType[1] = "Y" ;
1176 }
1177 else if ( poltype == "circular" ) {
1178 polType[0] = "R" ;
1179 polType[1] = "L" ;
1180 }
1181 else {
1182 polType[0] = "X" ;
1183 polType[1] = "Y" ;
1184 }
1185 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
1186
1187 TableRow tr( feedtab ) ;
1188 TableRecord &r = tr.record() ;
1189 putField( "FEED_ID", r, fid ) ;
1190 putField( "BEAM_ID", r, fid ) ;
1191 Int tmp = 0 ;
1192 putField( "ANTENNA_ID", r, tmp ) ;
1193 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
1194 putField( "NUM_RECEPTORS", r, numReceptors ) ;
1195 defineField( "POLARIZATION_TYPE", r, polType ) ;
1196 defineField( "BEAM_OFFSET", r, beamOffset ) ;
1197 defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ;
1198 defineField( "POL_RESPONSE", r, polResponse ) ;
1199 tr.put( nrow ) ;
1200
1201 feedEntry.resize( len+nItem ) ;
1202 feedEntry[len] = fid ;
1203 feedEntry[len+1] = sid ;
1204 }
1205 }
1206 void initPolarization()
1207 {
1208 const TableRecord &keys = table.keywordSet() ;
1209 poltype = keys.asString( "POLTYPE" ) ;
1210
1211 initCorrProductTemplate() ;
1212 }
1213 void initFrequencies()
1214 {
1215 const TableRecord &keys = table.keywordSet() ;
1216 Table tab = keys.asTable( "FREQUENCIES" ) ;
1217 ROScalarColumn<uInt> idcol( tab, "ID" ) ;
1218 ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ;
1219 ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ;
1220 ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ;
1221 Vector<uInt> id = idcol.getColumn() ;
1222 Vector<Double> rp = rpcol.getColumn() ;
1223 Vector<Double> rv = rvcol.getColumn() ;
1224 Vector<Double> ic = iccol.getColumn() ;
1225 for ( uInt i = 0 ; i < id.nelements() ; i++ ) {
1226 refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ;
1227 refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ;
1228 increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ;
1229 }
1230 String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ;
1231 MFrequency::getType( freqframe, frameStr ) ;
1232 }
1233 void attachSubtables()
1234 {
1235 //const TableRecord &keys = table.keywordSet() ;
1236 TableRecord &mskeys = ms.rwKeywordSet() ;
1237
1238 // FIELD table
1239 fieldtab = mskeys.asTable( "FIELD" ) ;
1240
1241 // SPECTRAL_WINDOW table
1242 spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ;
1243
1244 // POINTING table
1245 potab = mskeys.asTable( "POINTING" ) ;
1246
1247 // POLARIZATION table
1248 poltab = mskeys.asTable( "POLARIZATION" ) ;
1249
1250 // DATA_DESCRIPTION table
1251 ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ;
1252
1253 // STATE table
1254 statetab = mskeys.asTable( "STATE" ) ;
1255
1256 // FEED table
1257 feedtab = mskeys.asTable( "FEED" ) ;
1258 }
1259 void attachMain()
1260 {
1261 TableRecord &r = row.record() ;
1262 dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ;
1263 timeRF.attachToRecord( r, "TIME" ) ;
1264 timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ;
1265 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1266 exposureRF.attachToRecord( r, "EXPOSURE" ) ;
1267 fieldIdRF.attachToRecord( r, "FIELD_ID" ) ;
1268 feed1RF.attachToRecord( r, "FEED1" ) ;
1269 feed2RF.attachToRecord( r, "FEED2" ) ;
1270 scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ;
1271 stateIdRF.attachToRecord( r, "STATE_ID" ) ;
1272
1273 // constant values
1274 //Int id = 0 ;
1275 RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ;
1276 *intRF = 0 ;
1277 intRF.attachToRecord( r, "ANTENNA1" ) ;
1278 *intRF = 0 ;
1279 intRF.attachToRecord( r, "ANTENNA2" ) ;
1280 *intRF = 0 ;
1281 intRF.attachToRecord( r, "ARRAY_ID" ) ;
1282 *intRF = 0 ;
1283 intRF.attachToRecord( r, "PROCESSOR_ID" ) ;
1284 *intRF = 0 ;
1285 RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ;
1286 arrayRF.define( Vector<Double>( 3, 0.0 ) ) ;
1287 }
1288 void attachPointing()
1289 {
1290 porow = TableRow( potab ) ;
1291 TableRecord &r = porow.record() ;
1292 poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ;
1293 poTimeRF.attachToRecord( r, "TIME" ) ;
1294 poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ;
1295 poIntervalRF.attachToRecord( r, "INTERVAL" ) ;
1296 poNameRF.attachToRecord( r, "NAME" ) ;
1297 poDirectionRF.attachToRecord( r, "DIRECTION" ) ;
1298 poTargetRF.attachToRecord( r, "TARGET" ) ;
1299
1300 // constant values
1301 RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ;
1302 *antIdRF = 0 ;
1303 RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ;
1304 *trackingRF = True ;
1305 }
1306 void queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1307 {
1308 t = 0.0 ;
1309 l = 0.0 ;
1310
1311 String sep1="#" ;
1312 String sep2="," ;
1313 String target="OBSERVE_TARGET" ;
1314 String atmcal="CALIBRATE_TEMPERATURE" ;
1315 String onstr="ON_SOURCE" ;
1316 String offstr="OFF_SOURCE" ;
1317 String pswitch="POSITION_SWITCH" ;
1318 String nod="NOD" ;
1319 String fswitch="FREQUENCY_SWITCH" ;
1320 String sigstr="SIG" ;
1321 String refstr="REF" ;
1322 String unspecified="UNSPECIFIED" ;
1323 String ftlow="LOWER" ;
1324 String fthigh="HIGHER" ;
1325 switch ( type ) {
1326 case SrcType::PSON:
1327 stype = target+sep1+onstr+sep2+pswitch ;
1328 b = True ;
1329 break ;
1330 case SrcType::PSOFF:
1331 stype = target+sep1+offstr+sep2+pswitch ;
1332 b = False ;
1333 break ;
1334 case SrcType::NOD:
1335 stype = target+sep1+onstr+sep2+nod ;
1336 b = True ;
1337 break ;
1338 case SrcType::FSON:
1339 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1340 b = True ;
1341 break ;
1342 case SrcType::FSOFF:
1343 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
1344 b = False ;
1345 break ;
1346 case SrcType::SKY:
1347 stype = atmcal+sep1+offstr+sep2+unspecified ;
1348 b = False ;
1349 break ;
1350 case SrcType::HOT:
1351 stype = atmcal+sep1+offstr+sep2+unspecified ;
1352 b = False ;
1353 break ;
1354 case SrcType::WARM:
1355 stype = atmcal+sep1+offstr+sep2+unspecified ;
1356 b = False ;
1357 break ;
1358 case SrcType::COLD:
1359 stype = atmcal+sep1+offstr+sep2+unspecified ;
1360 b = False ;
1361 break ;
1362 case SrcType::PONCAL:
1363 stype = atmcal+sep1+onstr+sep2+pswitch ;
1364 b = True ;
1365 break ;
1366 case SrcType::POFFCAL:
1367 stype = atmcal+sep1+offstr+sep2+pswitch ;
1368 b = False ;
1369 break ;
1370 case SrcType::NODCAL:
1371 stype = atmcal+sep1+onstr+sep2+nod ;
1372 b = True ;
1373 break ;
1374 case SrcType::FONCAL:
1375 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1376 b = True ;
1377 break ;
1378 case SrcType::FOFFCAL:
1379 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
1380 b = False ;
1381 break ;
1382 case SrcType::FSLO:
1383 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
1384 b = True ;
1385 break ;
1386 case SrcType::FLOOFF:
1387 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1388 b = False ;
1389 break ;
1390 case SrcType::FLOSKY:
1391 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1392 b = False ;
1393 break ;
1394 case SrcType::FLOHOT:
1395 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1396 b = False ;
1397 break ;
1398 case SrcType::FLOWARM:
1399 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1400 b = False ;
1401 break ;
1402 case SrcType::FLOCOLD:
1403 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1404 b = False ;
1405 break ;
1406 case SrcType::FSHI:
1407 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
1408 b = True ;
1409 break ;
1410 case SrcType::FHIOFF:
1411 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1412 b = False ;
1413 break ;
1414 case SrcType::FHISKY:
1415 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1416 b = False ;
1417 break ;
1418 case SrcType::FHIHOT:
1419 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1420 b = False ;
1421 break ;
1422 case SrcType::FHIWARM:
1423 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1424 b = False ;
1425 break ;
1426 case SrcType::FHICOLD:
1427 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1428 b = False ;
1429 break ;
1430 case SrcType::SIG:
1431 stype = target+sep1+onstr+sep2+unspecified ;
1432 b = True ;
1433 break ;
1434 case SrcType::REF:
1435 stype = target+sep1+offstr+sep2+unspecified ;
1436 b = False ;
1437 break ;
1438 default:
1439 stype = unspecified ;
1440 b = True ;
1441 break ;
1442 }
1443 }
1444 void initCorrProductTemplate()
1445 {
1446 Int n = 1 ;
1447 {
1448 Matrix<Int> c( 2, n, 0 ) ;
1449 corrProductTemplate[n] = c ;
1450 }
1451 n = 2 ;
1452 {
1453 Matrix<Int> c( 2, n, 0 ) ;
1454 c.column( 1 ) = 1 ;
1455 corrProductTemplate[n] = c ;
1456 }
1457 n = 4 ;
1458 {
1459 Matrix<Int> c( 2, n, 0 ) ;
1460 c( 0, 2 ) = 1 ;
1461 c( 0, 3 ) = 1 ;
1462 c( 1, 1 ) = 1 ;
1463 c( 1, 3 ) = 1 ;
1464 corrProductTemplate[n] = c ;
1465 }
1466 }
1467
1468 Table &ms;
1469 TableRow row;
1470 uInt rowidx;
1471 String fieldName;
1472 Int fieldId;
1473 Int srcId;
1474 Int defaultFieldId;
1475 Int spwId;
1476 Int feedId;
1477 Int subscan;
1478 CountedPtr<DataHolder> holder;
1479 String ptName;
1480 Bool useFloat;
1481 String poltype;
1482
1483 // MS subtables
1484 Table spwtab;
1485 Table statetab;
1486 Table ddtab;
1487 Table poltab;
1488 Table fieldtab;
1489 Table feedtab;
1490 Table potab;
1491
1492 // Scantable MAIN columns
1493 ROArrayColumn<Float> spectraCol;
1494 ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol;
1495 ROArrayColumn<uChar> flagtraCol;
1496 ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol,
1497 sourceNameCol,fieldNameCol;
1498
1499 // MS MAIN columns
1500 RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF,
1501 scanNumberRF,stateIdRF;
1502 RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF;
1503
1504 // MS POINTING columns
1505 TableRow porow;
1506 RecordFieldPtr<Int> poNumPolyRF ;
1507 RecordFieldPtr<Double> poTimeRF,
1508 poTimeOriginRF,
1509 poIntervalRF ;
1510 RecordFieldPtr<String> poNameRF ;
1511 RecordFieldPtr< Matrix<Double> > poDirectionRF,
1512 poTargetRF ;
1513
1514 Vector<String> stateEntry;
1515 Block<Int> ddEntry;
1516 Block<Int> feedEntry;
1517 vector< Vector<Int> > polEntry;
1518 set<uInt> processedFreqId;
1519 set<uInt> processedIFNO;
1520 map<uInt,Double> refpix;
1521 map<uInt,Double> refval;
1522 map<uInt,Double> increment;
1523 MFrequency::Types freqframe;
1524 Record srcRec;
1525 map< Int, Matrix<Int> > corrProductTemplate;
1526};
1527
1528class BaseMSSysCalVisitor: public TableVisitor {
1529 uInt lastRecordNo;
1530 uInt lastBeamNo, lastIfNo, lastPolNo;
1531 Double lastTime;
1532protected:
1533 const Table &table;
1534 uInt count;
1535public:
1536 BaseMSSysCalVisitor(const Table &table)
1537 : table(table)
1538 {
1539 count = 0;
1540 }
1541
1542 virtual void enterBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1543 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1544 virtual void enterIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1545 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1546 virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1547 virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1548 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1549 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1550
1551 virtual Bool visitRecord(const uInt /*recordNo*/,
1552 const uInt /*beamNo*/,
1553 const uInt /*ifNo*/,
1554 const uInt /*polNo*/,
1555 const Double /*time*/) { return True ;}
1556
1557 virtual Bool visit(Bool isFirst, const uInt recordNo,
1558 const uInt nCols, void const *const colValues[]) {
1559 uInt beamNo, ifNo, polNo;
1560 Double time;
1561 { // prologue
1562 uInt i = 0;
1563 {
1564 const uInt *col = (const uInt *)colValues[i++];
1565 beamNo = col[recordNo];
1566 }
1567 {
1568 const uInt *col = (const uInt *)colValues[i++];
1569 ifNo = col[recordNo];
1570 }
1571 {
1572 const Double *col = (const Double *)colValues[i++];
1573 time = col[recordNo];
1574 }
1575 {
1576 const uInt *col = (const uInt *)colValues[i++];
1577 polNo = col[recordNo];
1578 }
1579 assert(nCols == i);
1580 }
1581
1582 if (isFirst) {
1583 enterBeamNo(recordNo, beamNo);
1584 enterIfNo(recordNo, ifNo);
1585 enterTime(recordNo, time);
1586 enterPolNo(recordNo, polNo);
1587 } else {
1588 if (lastBeamNo != beamNo) {
1589 leavePolNo(lastRecordNo, lastPolNo);
1590 leaveTime(lastRecordNo, lastTime);
1591 leaveIfNo(lastRecordNo, lastIfNo);
1592 leaveBeamNo(lastRecordNo, lastBeamNo);
1593
1594 enterBeamNo(recordNo, beamNo);
1595 enterIfNo(recordNo, ifNo);
1596 enterTime(recordNo, time);
1597 enterPolNo(recordNo, polNo);
1598 } else if (lastIfNo != ifNo) {
1599 leavePolNo(lastRecordNo, lastPolNo);
1600 leaveTime(lastRecordNo, lastTime);
1601 leaveIfNo(lastRecordNo, lastIfNo);
1602
1603 enterIfNo(recordNo, ifNo);
1604 enterTime(recordNo, time);
1605 enterPolNo(recordNo, polNo);
1606 } else if (lastTime != time) {
1607 leavePolNo(lastRecordNo, lastPolNo);
1608 leaveTime(lastRecordNo, lastTime);
1609
1610 enterTime(recordNo, time);
1611 enterPolNo(recordNo, polNo);
1612 } else if (lastPolNo != polNo) {
1613 leavePolNo(lastRecordNo, lastPolNo);
1614 enterPolNo(recordNo, polNo);
1615 }
1616 }
1617 count++;
1618 Bool result = visitRecord(recordNo, beamNo, ifNo, polNo, time);
1619
1620 { // epilogue
1621 lastRecordNo = recordNo;
1622
1623 lastBeamNo = beamNo;
1624 lastIfNo = ifNo;
1625 lastPolNo = polNo;
1626 lastTime = time;
1627 }
1628 return result ;
1629 }
1630
1631 virtual void finish() {
1632 if (count > 0) {
1633 leavePolNo(lastRecordNo, lastPolNo);
1634 leaveTime(lastRecordNo, lastTime);
1635 leaveIfNo(lastRecordNo, lastIfNo);
1636 leaveBeamNo(lastRecordNo, lastBeamNo);
1637 }
1638 }
1639};
1640
1641class BaseTsysHolder
1642{
1643public:
1644 BaseTsysHolder( ROArrayColumn<Float> &tsysCol )
1645 : col( tsysCol ),
1646 nchan(0)
1647 {
1648 reset() ;
1649 }
1650 virtual ~BaseTsysHolder() {}
1651 virtual Array<Float> getTsys() = 0 ;
1652 void setNchan( uInt n ) { nchan = n ; }
1653 void appendTsys( uInt row )
1654 {
1655 Vector<Float> v = col( row ) ;
1656 uInt len = tsys.nrow() ;
1657 tsys.resize( len+1, nchan, True ) ;
1658 if ( v.nelements() == nchan )
1659 tsys.row( len ) = v ;
1660 else
1661 tsys.row( len ) = v[0] ;
1662 }
1663 void setTsys( uInt row, uInt idx )
1664 {
1665 if ( idx >= nrow() )
1666 appendTsys( row ) ;
1667 else {
1668 Vector<Float> v = col( row ) ;
1669 if ( v.nelements() == nchan )
1670 tsys.row( idx ) = v ;
1671 else
1672 tsys.row( idx ) = v[0] ;
1673 }
1674 }
1675 void reset()
1676 {
1677 tsys.resize() ;
1678 }
1679 uInt nrow() { return tsys.nrow() ; }
1680 Bool isEffective()
1681 {
1682 return ( !(tsys.empty()) && anyNE( tsys, (Float)1.0 ) ) ;
1683 }
1684 BaseTsysHolder &operator= ( const BaseTsysHolder &v )
1685 {
1686 if ( this != &v )
1687 tsys.assign( v.tsys ) ;
1688 return *this ;
1689 }
1690protected:
1691 ROArrayColumn<Float> col ;
1692 Matrix<Float> tsys ;
1693 uInt nchan ;
1694};
1695
1696class TsysHolder : public BaseTsysHolder
1697{
1698public:
1699 TsysHolder( ROArrayColumn<Float> &tsysCol )
1700 : BaseTsysHolder( tsysCol )
1701 {}
1702 virtual ~TsysHolder() {}
1703 virtual Array<Float> getTsys()
1704 {
1705 return tsys.column( 0 ) ;
1706 }
1707};
1708
1709class TsysSpectrumHolder : public BaseTsysHolder
1710{
1711public:
1712 TsysSpectrumHolder( ROArrayColumn<Float> &tsysCol )
1713 : BaseTsysHolder( tsysCol )
1714 {}
1715 virtual ~TsysSpectrumHolder() {}
1716 virtual Array<Float> getTsys()
1717 {
1718 return tsys ;
1719 }
1720};
1721
1722class BaseTcalProcessor
1723{
1724public:
1725 BaseTcalProcessor( ROArrayColumn<Float> &tcalCol )
1726 : col_( tcalCol )
1727 {}
1728 virtual ~BaseTcalProcessor() {}
1729 void setTcalId( Vector<uInt> &tcalId ) { id_.assign( tcalId ) ; }
1730 virtual Array<Float> getTcal() = 0 ;
1731protected:
1732 ROArrayColumn<Float> col_ ;
1733 Vector<uInt> id_ ;
1734};
1735
1736class TcalProcessor : public BaseTcalProcessor
1737{
1738public:
1739 TcalProcessor( ROArrayColumn<Float> &tcalCol )
1740 : BaseTcalProcessor( tcalCol )
1741 {}
1742 virtual ~TcalProcessor() {}
1743 virtual Array<Float> getTcal()
1744 {
1745 uInt npol = id_.nelements() ;
1746 Vector<Float> tcal( npol ) ;
1747 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
1748 tcal[ipol] = col_( id_[ipol] ).data()[0] ;
1749 //cout << "TcalProcessor: tcal = " << tcal << endl ;
1750 return tcal ;
1751 }
1752};
1753
1754class TcalSpectrumProcessor : public BaseTcalProcessor
1755{
1756public:
1757 TcalSpectrumProcessor( ROArrayColumn<Float> &tcalCol )
1758 : BaseTcalProcessor( tcalCol )
1759 {}
1760 virtual ~TcalSpectrumProcessor() {}
1761 virtual Array<Float> getTcal()
1762 {
1763 uInt npol = id_.nelements() ;
1764 //Vector<Float> tcal0 = col_( 0 ) ;
1765 Vector<Float> tcal0 = col_( id_[0] ) ;
1766 uInt nchan = tcal0.nelements() ;
1767 Matrix<Float> tcal( npol, nchan ) ;
1768 tcal.row( 0 ) = tcal0 ;
1769 for ( uInt ipol = 1 ; ipol < npol ; ipol++ )
1770 tcal.row( ipol ) = col_( id_[ipol] ) ;
1771 return tcal ;
1772 }
1773};
1774
1775class MSSysCalVisitor : public BaseMSSysCalVisitor
1776{
1777public:
1778 MSSysCalVisitor( const Table &from, Table &to )
1779 : BaseMSSysCalVisitor( from ),
1780 sctab( to ),
1781 rowidx( 0 ),
1782 polno_()
1783 {
1784 scrow = TableRow( sctab ) ;
1785
1786 lastTcalId.resize() ;
1787 theTcalId.resize() ;
1788 startTime = 0.0 ;
1789 endTime = 0.0 ;
1790
1791 const TableRecord &keys = table.keywordSet() ;
1792 Table tcalTable = keys.asTable( "TCAL" ) ;
1793 tcalCol.attach( tcalTable, "TCAL" ) ;
1794 tsysCol.attach( table, "TSYS" ) ;
1795 tcalIdCol.attach( table, "TCAL_ID" ) ;
1796 intervalCol.attach( table, "INTERVAL" ) ;
1797 effectiveTcal.resize( tcalTable.nrow() ) ;
1798 for ( uInt irow = 0 ; irow < tcalTable.nrow() ; irow++ ) {
1799 if ( allEQ( tcalCol( irow ), (Float)1.0 ) )
1800 effectiveTcal[irow] = False ;
1801 else
1802 effectiveTcal[irow] = True ;
1803 }
1804
1805 TableRecord &r = scrow.record() ;
1806 RecordFieldPtr<Int> antennaIdRF( r, "ANTENNA_ID" ) ;
1807 *antennaIdRF = 0 ;
1808 feedIdRF.attachToRecord( r, "FEED_ID" ) ;
1809 specWinIdRF.attachToRecord( r, "SPECTRAL_WINDOW_ID" ) ;
1810 timeRF.attachToRecord( r, "TIME" ) ;
1811 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1812 if ( r.isDefined( "TCAL" ) ) {
1813 tcalRF.attachToRecord( r, "TCAL" ) ;
1814 tcalProcessor = new TcalProcessor( tcalCol ) ;
1815 }
1816 else if ( r.isDefined( "TCAL_SPECTRUM" ) ) {
1817 tcalRF.attachToRecord( r, "TCAL_SPECTRUM" ) ;
1818 tcalProcessor = new TcalSpectrumProcessor( tcalCol ) ;
1819 }
1820 if ( r.isDefined( "TSYS" ) ) {
1821 tsysRF.attachToRecord( r, "TSYS" ) ;
1822 theTsys = new TsysHolder( tsysCol ) ;
1823 lastTsys = new TsysHolder( tsysCol ) ;
1824 }
1825 else {
1826 tsysRF.attachToRecord( r, "TSYS_SPECTRUM" ) ;
1827 theTsys = new TsysSpectrumHolder( tsysCol ) ;
1828 lastTsys = new TsysSpectrumHolder( tsysCol ) ;
1829 }
1830
1831 }
1832
1833 virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue)
1834 {
1835 *feedIdRF = (Int)columnValue ;
1836 }
1837 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/)
1838 {
1839 }
1840 virtual void enterIfNo(const uInt recordNo, uInt columnValue)
1841 {
1842 //cout << "enterIfNo" << endl ;
1843 ROArrayColumn<Float> sp( table, "SPECTRA" ) ;
1844 uInt nchan = sp( recordNo ).nelements() ;
1845 theTsys->setNchan( nchan ) ;
1846 lastTsys->setNchan( nchan ) ;
1847
1848 *specWinIdRF = (Int)columnValue ;
1849 }
1850 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/)
1851 {
1852 //cout << "leaveIfNo" << endl ;
1853 post() ;
1854 reset(true);
1855 startTime = 0.0 ;
1856 endTime = 0.0 ;
1857 }
1858 virtual void enterTime(const uInt recordNo, Double columnValue)
1859 {
1860 //cout << "enterTime" << endl ;
1861 interval = intervalCol.asdouble( recordNo ) ;
1862 // start time and end time
1863 if ( startTime == 0.0 ) {
1864 startTime = columnValue * 86400.0 - 0.5 * interval ;
1865 endTime = columnValue * 86400.0 + 0.5 * interval ;
1866 }
1867 }
1868 virtual void leaveTime(const uInt /*recordNo*/, Double columnValue)
1869 {
1870 //cout << "leaveTime" << endl ;
1871 if ( isUpdated() ) {
1872 post() ;
1873 *lastTsys = *theTsys ;
1874 lastTcalId = theTcalId ;
1875 reset(false);
1876 startTime = columnValue * 86400.0 - 0.5 * interval ;
1877 endTime = columnValue * 86400.0 + 0.5 * interval ;
1878 }
1879 else {
1880 endTime = columnValue * 86400.0 + 0.5 * interval ;
1881 }
1882 }
1883 virtual void enterPolNo(const uInt recordNo, uInt columnValue)
1884 {
1885 //cout << "enterPolNo" << endl ;
1886 Vector<Float> tsys = tsysCol( recordNo ) ;
1887 uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
1888 polno_.insert( columnValue ) ;
1889 uInt numPol = polno_.size() ;
1890 if ( lastTsys->nrow() < numPol )
1891 lastTsys->appendTsys( recordNo ) ;
1892 if ( lastTcalId.nelements() <= numPol )
1893 appendTcalId( lastTcalId, tcalId, numPol-1 ) ;
1894 if ( theTsys->nrow() < numPol )
1895 theTsys->appendTsys( recordNo ) ;
1896 else {
1897 theTsys->setTsys( recordNo, numPol-1 ) ;
1898 }
1899 if ( theTcalId.nelements() < numPol )
1900 appendTcalId( theTcalId, tcalId, numPol-1 ) ;
1901 else
1902 setTcalId( theTcalId, tcalId, numPol-1 ) ;
1903 }
1904 virtual void leavePolNo( const uInt /*recordNo*/, uInt /*columnValue*/ )
1905 {
1906 }
1907
1908private:
1909 void reset(bool completely)
1910 {
1911 if (completely) {
1912 lastTsys->reset() ;
1913 lastTcalId.resize() ;
1914 }
1915 theTsys->reset() ;
1916 theTcalId.resize() ;
1917 polno_.clear();
1918 }
1919 void appendTcalId( Vector<uInt> &v, uInt &elem, uInt polId )
1920 {
1921 v.resize( polId+1, True ) ;
1922 v[polId] = elem ;
1923 }
1924 void setTcalId( Vector<uInt> &v, uInt &elem, uInt polId )
1925 {
1926 v[polId] = elem ;
1927 }
1928 void post()
1929 {
1930 // check if given Tcal and Tsys is effective
1931 Bool isEffective = False ;
1932 for ( uInt ipol = 0 ; ipol < lastTcalId.nelements() ; ipol++ ) {
1933 if ( effectiveTcal[lastTcalId[ipol]] ) {
1934 isEffective = True ;
1935 break ;
1936 }
1937 }
1938 if ( !isEffective ) {
1939 if ( !(lastTsys->isEffective()) )
1940 return ;
1941 }
1942
1943 //cout << " interval: " << (endTime-startTime) << " lastTcalId = " << lastTcalId << endl ;
1944 Double midTime = 0.5 * ( startTime + endTime ) ;
1945 Double interval = endTime - startTime ;
1946 *timeRF = midTime ;
1947 *intervalRF = interval ;
1948 tcalProcessor->setTcalId( lastTcalId ) ;
1949 Array<Float> tcal = tcalProcessor->getTcal() ;
1950 Array<Float> tsys = lastTsys->getTsys() ;
1951 tcalRF.define( tcal ) ;
1952 tsysRF.define( tsys ) ;
1953 sctab.addRow( 1, True ) ;
1954 scrow.put( rowidx ) ;
1955 rowidx++ ;
1956 }
1957
1958 Bool isUpdated()
1959 {
1960 Bool ret = (anyNE( theTcalId, lastTcalId ) || anyNE( theTsys->getTsys(), lastTsys->getTsys() )) ;
1961 return ret ;
1962 }
1963
1964 Table &sctab;
1965 TableRow scrow;
1966 uInt rowidx;
1967
1968 Double startTime,endTime,interval;
1969
1970 CountedPtr<BaseTsysHolder> lastTsys,theTsys;
1971 Vector<uInt> lastTcalId,theTcalId;
1972 set<uInt> polno_;
1973 CountedPtr<BaseTcalProcessor> tcalProcessor ;
1974 Vector<Bool> effectiveTcal;
1975
1976 RecordFieldPtr<Int> feedIdRF,specWinIdRF;
1977 RecordFieldPtr<Double> timeRF,intervalRF;
1978 RecordFieldPtr< Array<Float> > tcalRF,tsysRF;
1979
1980 ROArrayColumn<Float> tsysCol,tcalCol;
1981 ROTableColumn tcalIdCol,intervalCol;
1982};
1983
1984MSWriter::MSWriter(CountedPtr<Scantable> stable)
1985 : table_(stable),
1986 mstable_(NULL),
1987 isWeather_(False),
1988 tcalSpec_(False),
1989 tsysSpec_(False),
1990 ptTabName_("")
1991{
1992 os_ = LogIO() ;
1993 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
1994// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
1995
1996 // initialize writer
1997 init() ;
1998}
1999
2000MSWriter::~MSWriter()
2001{
2002 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
2003// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
2004
2005 if ( mstable_ != 0 )
2006 delete mstable_ ;
2007}
2008
2009bool MSWriter::write(const string& filename, const Record& rec)
2010{
2011 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
2012 //double startSec = mathutil::gettimeofday_sec() ;
2013 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
2014
2015 filename_ = filename ;
2016
2017 // parsing MS options
2018 Bool overwrite = False ;
2019 if ( rec.isDefined( "ms" ) ) {
2020 Record msrec = rec.asRecord( "ms" ) ;
2021 if ( msrec.isDefined( "overwrite" ) ) {
2022 overwrite = msrec.asBool( "overwrite" ) ;
2023 }
2024 }
2025
2026 os_ << "Parsing MS options" << endl ;
2027 os_ << " overwrite = " << overwrite << LogIO::POST ;
2028
2029 File file( filename_ ) ;
2030 if ( file.exists() ) {
2031 if ( overwrite ) {
2032 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
2033 if ( file.isRegular() ) RegularFile(file).remove() ;
2034 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
2035 else SymLink(file).remove() ;
2036 }
2037 else {
2038 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
2039 return False ;
2040 }
2041 }
2042
2043 // set up MS
2044 setupMS() ;
2045
2046 // subtables
2047 // OBSERVATION
2048 fillObservation() ;
2049
2050 // ANTENNA
2051 fillAntenna() ;
2052
2053 // PROCESSOR
2054 fillProcessor() ;
2055
2056 // SOURCE
2057 fillSource() ;
2058
2059 // WEATHER
2060 if ( isWeather_ )
2061 fillWeather() ;
2062
2063 // SYSCAL
2064 fillSysCal() ;
2065
2066 /***
2067 * Start iteration using TableVisitor
2068 ***/
2069 {
2070 static const char *cols[] = {
2071 "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME",
2072 "POLNO",
2073 NULL
2074 };
2075 static const TypeManagerImpl<uInt> tmUInt;
2076 static const TypeManagerImpl<Int> tmInt;
2077 static const TypeManagerImpl<Double> tmDouble;
2078 static const TypeManagerImpl<String> tmString;
2079 static const TypeManager *const tms[] = {
2080 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL
2081 };
2082 //double t0 = mathutil::gettimeofday_sec() ;
2083 MSWriterVisitor myVisitor(table_->table(),*mstable_);
2084 //double t1 = mathutil::gettimeofday_sec() ;
2085 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
2086 String dataColName = "FLOAT_DATA" ;
2087 if ( useData_ )
2088 dataColName = "DATA" ;
2089 myVisitor.dataColumnName( dataColName ) ;
2090 myVisitor.pointingTableName( ptTabName_ ) ;
2091 myVisitor.setSourceRecord( srcRec_ ) ;
2092 //double t2 = mathutil::gettimeofday_sec() ;
2093 traverseTable(table_->table(), cols, tms, &myVisitor);
2094 //double t3 = mathutil::gettimeofday_sec() ;
2095 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
2096 }
2097 /***
2098 * End iteration using TableVisitor
2099 ***/
2100
2101 // ASDM tables
2102 const TableRecord &stKeys = table_->table().keywordSet() ;
2103 TableRecord &msKeys = mstable_->rwKeywordSet() ;
2104 uInt nfields = stKeys.nfields() ;
2105 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
2106 String kname = stKeys.name( ifield ) ;
2107 if ( kname.find( "ASDM" ) != String::npos ) {
2108 String asdmpath = stKeys.asString( ifield ) ;
2109 os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
2110 if ( Table::isReadable( asdmpath ) ) {
2111 Table newAsdmTab( asdmpath, Table::Old ) ;
2112 newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
2113 os_ << "add subtable: " << kname << LogIO::POST ;
2114 msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
2115 }
2116 }
2117 }
2118
2119 // replace POINTING table with original one if exists
2120 if ( ptTabName_ != "" ) {
2121 delete mstable_ ;
2122 mstable_ = 0 ;
2123 Table newPtTab( ptTabName_, Table::Old ) ;
2124 newPtTab.copy( filename_+"/POINTING", Table::New ) ;
2125 }
2126
2127 //double endSec = mathutil::gettimeofday_sec() ;
2128 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2129
2130 os_ << "Exported data as MS" << LogIO::POST ;
2131
2132 return True ;
2133}
2134
2135void MSWriter::init()
2136{
2137// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
2138// double startSec = mathutil::gettimeofday_sec() ;
2139// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
2140
2141 // access to scantable
2142 header_ = table_->getHeader() ;
2143
2144 // FLOAT_DATA? or DATA?
2145 if ( header_.npol > 2 ) {
2146 useFloatData_ = False ;
2147 useData_ = True ;
2148 }
2149 else {
2150 useFloatData_ = True ;
2151 useData_ = False ;
2152 }
2153
2154 // polarization type
2155 polType_ = header_.poltype ;
2156 if ( polType_ == "" )
2157 polType_ = "stokes" ;
2158 else if ( polType_.find( "linear" ) != String::npos )
2159 polType_ = "linear" ;
2160 else if ( polType_.find( "circular" ) != String::npos )
2161 polType_ = "circular" ;
2162 else if ( polType_.find( "stokes" ) != String::npos )
2163 polType_ = "stokes" ;
2164 else if ( polType_.find( "linpol" ) != String::npos )
2165 polType_ = "linpol" ;
2166 else
2167 polType_ = "notype" ;
2168
2169 // Check if some subtables are exists
2170 Bool isTcal = False ;
2171 if ( table_->tcal().table().nrow() != 0 ) {
2172 ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
2173 if ( col.isDefined( 0 ) ) {
2174 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
2175 isTcal = True ;
2176 }
2177 else {
2178 os_ << "No TCAL rows" << LogIO::POST ;
2179 }
2180 }
2181 else {
2182 os_ << "No TCAL rows" << LogIO::POST ;
2183 }
2184 if ( table_->weather().table().nrow() != 0 ) {
2185 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
2186 if ( col.isDefined( 0 ) ) {
2187 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
2188 isWeather_ =True ;
2189 }
2190 else {
2191 os_ << "No WEATHER rows" << LogIO::POST ;
2192 }
2193 }
2194 else {
2195 os_ << "No WEATHER rows" << LogIO::POST ;
2196 }
2197
2198 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
2199 if ( header_.nchan != 1 ) {
2200 if ( isTcal ) {
2201 // examine TCAL subtable
2202 Table tcaltab = table_->tcal().table() ;
2203 ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
2204 for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
2205 if ( tcalCol( irow ).size() != 1 )
2206 tcalSpec_ = True ;
2207 }
2208 }
2209 // examine spectral data
2210 TableIterator iter0( table_->table(), "IFNO" ) ;
2211 while( !iter0.pastEnd() ) {
2212 Table t0( iter0.table() ) ;
2213 ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
2214 uInt len = sharedFloatArrCol( 0 ).size() ;
2215 if ( len != 1 ) {
2216 sharedFloatArrCol.attach( t0, "TSYS" ) ;
2217 if ( sharedFloatArrCol( 0 ).size() != 1 )
2218 tsysSpec_ = True ;
2219 }
2220 iter0.next() ;
2221 }
2222 }
2223
2224 // check if reference for POINTING table exists
2225 const TableRecord &rec = table_->table().keywordSet() ;
2226 if ( rec.isDefined( "POINTING" ) ) {
2227 ptTabName_ = rec.asString( "POINTING" ) ;
2228 if ( !Table::isReadable( ptTabName_ ) ) {
2229 ptTabName_ = "" ;
2230 }
2231 }
2232
2233// double endSec = mathutil::gettimeofday_sec() ;
2234// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2235}
2236
2237void MSWriter::setupMS()
2238{
2239// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
2240// double startSec = mathutil::gettimeofday_sec() ;
2241// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
2242
2243 String dunit = table_->getHeader().fluxunit ;
2244
2245 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
2246 if ( useFloatData_ )
2247 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
2248 else if ( useData_ )
2249 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
2250
2251 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
2252
2253 mstable_ = new MeasurementSet( newtab ) ;
2254
2255 TableColumn col ;
2256 if ( useFloatData_ )
2257 col.attach( *mstable_, "FLOAT_DATA" ) ;
2258 else if ( useData_ )
2259 col.attach( *mstable_, "DATA" ) ;
2260 col.rwKeywordSet().define( "UNIT", dunit ) ;
2261
2262 // create subtables
2263 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
2264 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
2265 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
2266
2267 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
2268 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
2269 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
2270
2271 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
2272 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
2273 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
2274
2275 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
2276 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
2277 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
2278
2279 TableDesc fieldDesc = MSField::requiredTableDesc() ;
2280 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
2281 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
2282
2283 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
2284 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
2285 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
2286
2287 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
2288 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
2289 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
2290
2291 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
2292 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
2293 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
2294
2295 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
2296 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
2297 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
2298
2299 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
2300 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
2301 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
2302
2303 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
2304 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
2305 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
2306
2307 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
2308 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
2309 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
2310
2311 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
2312 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
2313 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
2314 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
2315 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
2316 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
2317
2318 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
2319 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
2320 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
2321
2322 TableDesc stateDesc = MSState::requiredTableDesc() ;
2323 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
2324 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
2325
2326 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
2327 if ( tcalSpec_ )
2328 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
2329 else
2330 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ;
2331 if ( tsysSpec_ )
2332 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
2333 else
2334 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ;
2335 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
2336 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
2337
2338 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
2339 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
2340 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
2341 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
2342 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
2343 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
2344 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
2345 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
2346
2347 mstable_->initRefs() ;
2348
2349// double endSec = mathutil::gettimeofday_sec() ;
2350// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2351}
2352
2353void MSWriter::fillObservation()
2354{
2355 //double startSec = mathutil::gettimeofday_sec() ;
2356 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
2357
2358 // only 1 row
2359 mstable_->observation().addRow( 1, True ) ;
2360 MSObservationColumns msObsCols( mstable_->observation() ) ;
2361 msObsCols.observer().put( 0, header_.observer ) ;
2362 // tentatively put antennaname (from ANTENNA subtable)
2363 String hAntennaName = header_.antennaname ;
2364 String::size_type pos = hAntennaName.find( "//" ) ;
2365 String telescopeName ;
2366 if ( pos != String::npos ) {
2367 telescopeName = hAntennaName.substr( 0, pos ) ;
2368 }
2369 else {
2370 pos = hAntennaName.find( "@" ) ;
2371 telescopeName = hAntennaName.substr( 0, pos ) ;
2372 }
2373// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
2374 msObsCols.telescopeName().put( 0, telescopeName ) ;
2375 msObsCols.project().put( 0, header_.project ) ;
2376 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
2377 Table sortedtable = table_->table().sort("TIME") ;
2378 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
2379 Vector<MEpoch> trange( 2 ) ;
2380 trange[0] = timeCol( 0 ) ;
2381 trange[1] = timeCol( table_->nrow()-1 ) ;
2382 msObsCols.timeRangeMeas().put( 0, trange ) ;
2383
2384 //double endSec = mathutil::gettimeofday_sec() ;
2385 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2386}
2387
2388void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d )
2389{
2390 name.upcase() ;
2391
2392 m = "ALT-AZ" ;
2393 t = "GROUND-BASED" ;
2394 if ( name.matches( Regex( "DV[0-9]+$" ) )
2395 || name.matches( Regex( "DA[0-9]+$" ) )
2396 || name.matches( Regex( "PM[0-9]+$" ) ) )
2397 d = 12.0 ;
2398 else if ( name.matches( Regex( "CM[0-9]+$" ) ) )
2399 d = 7.0 ;
2400 else if ( name.contains( "GBT" ) )
2401 d = 104.9 ;
2402 else if ( name.contains( "MOPRA" ) )
2403 d = 22.0 ;
2404 else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) )
2405 d = 64.0 ;
2406 else if ( name.contains( "TIDBINBILLA" ) )
2407 d = 70.0 ;
2408 else if ( name.contains( "CEDUNA" ) )
2409 d = 30.0 ;
2410 else if ( name.contains( "HOBART" ) )
2411 d = 26.0 ;
2412 else if ( name.contains( "APEX" ) )
2413 d = 12.0 ;
2414 else if ( name.contains( "ASTE" ) )
2415 d = 10.0 ;
2416 else if ( name.contains( "NRO" ) )
2417 d = 45.0 ;
2418 else
2419 d = 1.0 ;
2420}
2421
2422void MSWriter::fillAntenna()
2423{
2424 //double startSec = mathutil::gettimeofday_sec() ;
2425 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
2426
2427 // only 1 row
2428 Table anttab = mstable_->antenna() ;
2429 anttab.addRow( 1, True ) ;
2430
2431 Table &table = table_->table() ;
2432 const TableRecord &keys = table.keywordSet() ;
2433 String hAntName = keys.asString( "AntennaName" ) ;
2434 String::size_type pos = hAntName.find( "//" ) ;
2435 String antennaName ;
2436 String stationName ;
2437 if ( pos != String::npos ) {
2438 stationName = hAntName.substr( 0, pos ) ;
2439 hAntName = hAntName.substr( pos+2 ) ;
2440 }
2441 pos = hAntName.find( "@" ) ;
2442 if ( pos != String::npos ) {
2443 antennaName = hAntName.substr( 0, pos ) ;
2444 stationName = hAntName.substr( pos+1 ) ;
2445 }
2446 else {
2447 antennaName = hAntName ;
2448 }
2449 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
2450
2451 String mount, atype ;
2452 Double diameter ;
2453 antennaProperty( antennaName, mount, atype, diameter ) ;
2454
2455 TableRow tr( anttab ) ;
2456 TableRecord &r = tr.record() ;
2457 RecordFieldPtr<String> nameRF( r, "NAME" ) ;
2458 RecordFieldPtr<String> stationRF( r, "STATION" ) ;
2459 RecordFieldPtr<String> mountRF( r, "MOUNT" ) ;
2460 RecordFieldPtr<String> typeRF( r, "TYPE" ) ;
2461 RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ;
2462 RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ;
2463 *nameRF = antennaName ;
2464 *mountRF = mount ;
2465 *typeRF = atype ;
2466 *dishDiameterRF = diameter ;
2467 *positionRF = antpos ;
2468 *stationRF = stationName ;
2469
2470 tr.put( 0 ) ;
2471
2472 //double endSec = mathutil::gettimeofday_sec() ;
2473 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2474}
2475
2476void MSWriter::fillProcessor()
2477{
2478// double startSec = mathutil::gettimeofday_sec() ;
2479// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
2480
2481 // only add empty 1 row
2482 MSProcessor msProc = mstable_->processor() ;
2483 msProc.addRow( 1, True ) ;
2484
2485// double endSec = mathutil::gettimeofday_sec() ;
2486// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2487}
2488
2489void MSWriter::fillSource()
2490{
2491// double startSec = mathutil::gettimeofday_sec() ;
2492// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
2493
2494 // access to MS SOURCE subtable
2495 MSSource msSrc = mstable_->source() ;
2496
2497 // access to MOLECULE subtable
2498 STMolecules stm = table_->molecules() ;
2499
2500 Int srcId = 0 ;
2501 Vector<Double> restFreq ;
2502 Vector<String> molName ;
2503 Vector<String> fMolName ;
2504
2505 // row based
2506 TableRow row( msSrc ) ;
2507 TableRecord &rec = row.record() ;
2508 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
2509 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
2510 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
2511 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
2512 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
2513 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
2514 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
2515 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
2516 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
2517 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
2518 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
2519
2520 //
2521 // ITERATION: SRCNAME
2522 //
2523 TableIterator iter0( table_->table(), "SRCNAME" ) ;
2524 while( !iter0.pastEnd() ) {
2525 //Table t0( iter0.table() ) ;
2526 Table t0 = iter0.table() ;
2527
2528 // get necessary information
2529 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
2530 String srcName = srcNameCol( 0 ) ;
2531 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
2532 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
2533 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
2534 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
2535 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
2536 Double srcVel = srcVelCol( 0 ) ;
2537 srcRec_.define( srcName, srcId ) ;
2538
2539 // NAME
2540 *nameRF = srcName ;
2541
2542 // SOURCE_ID
2543 *srcidRF = srcId ;
2544
2545 // PROPER_MOTION
2546 *srcpmRF = srcPM ;
2547
2548 // DIRECTION
2549 *srcdirRF = srcDir ;
2550
2551 //
2552 // ITERATION: MOLECULE_ID
2553 //
2554 TableIterator iter1( t0, "MOLECULE_ID" ) ;
2555 while( !iter1.pastEnd() ) {
2556 //Table t1( iter1.table() ) ;
2557 Table t1 = iter1.table() ;
2558
2559 // get necessary information
2560 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
2561 uInt molId = molIdCol( 0 ) ;
2562 stm.getEntry( restFreq, molName, fMolName, molId ) ;
2563
2564 uInt numFreq = restFreq.size() ;
2565
2566 // NUM_LINES
2567 *numlineRF = numFreq ;
2568
2569 // REST_FREQUENCY
2570 *restfreqRF = restFreq ;
2571
2572 // TRANSITION
2573 //*transitionRF = fMolName ;
2574 Vector<String> transition ;
2575 if ( fMolName.size() != 0 ) {
2576 transition = fMolName ;
2577 }
2578 else if ( molName.size() != 0 ) {
2579 transition = molName ;
2580 }
2581 else {
2582 transition.resize( numFreq ) ;
2583 transition = "" ;
2584 }
2585 *transitionRF = transition ;
2586
2587 // SYSVEL
2588 Vector<Double> sysvelArr( numFreq, srcVel ) ;
2589 *sysvelRF = sysvelArr ;
2590
2591 //
2592 // ITERATION: IFNO
2593 //
2594 TableIterator iter2( t1, "IFNO" ) ;
2595 while( !iter2.pastEnd() ) {
2596 //Table t2( iter2.table() ) ;
2597 Table t2 = iter2.table() ;
2598 uInt nrow = msSrc.nrow() ;
2599
2600 // get necessary information
2601 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
2602 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
2603 Double midTime ;
2604 Double interval ;
2605 getValidTimeRange( midTime, interval, t2 ) ;
2606
2607 // fill SPECTRAL_WINDOW_ID
2608 *spwidRF = ifno ;
2609
2610 // fill TIME, INTERVAL
2611 *timeRF = midTime ;
2612 *intervalRF = interval ;
2613
2614 // add row
2615 msSrc.addRow( 1, True ) ;
2616 row.put( nrow ) ;
2617
2618 iter2.next() ;
2619 }
2620
2621 iter1.next() ;
2622 }
2623
2624 // increment srcId if SRCNAME changed
2625 srcId++ ;
2626
2627 iter0.next() ;
2628 }
2629
2630// double endSec = mathutil::gettimeofday_sec() ;
2631// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2632}
2633
2634void MSWriter::fillWeather()
2635{
2636// double startSec = mathutil::gettimeofday_sec() ;
2637// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
2638
2639 // access to MS WEATHER subtable
2640 MSWeather msw = mstable_->weather() ;
2641
2642 // access to WEATHER subtable
2643 Table stw = table_->weather().table() ;
2644 uInt nrow = stw.nrow() ;
2645
2646 if ( nrow == 0 )
2647 return ;
2648
2649 msw.addRow( nrow, True ) ;
2650 MSWeatherColumns mswCols( msw ) ;
2651
2652 // ANTENNA_ID is always 0
2653 Vector<Int> antIdArr( nrow, 0 ) ;
2654 mswCols.antennaId().putColumn( antIdArr ) ;
2655
2656 // fill weather status
2657 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
2658 mswCols.temperature().putColumn( sharedFloatCol ) ;
2659 sharedFloatCol.attach( stw, "PRESSURE" ) ;
2660 mswCols.pressure().putColumn( sharedFloatCol ) ;
2661 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
2662 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
2663 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
2664 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
2665 sharedFloatCol.attach( stw, "WINDAZ" ) ;
2666 mswCols.windDirection().putColumn( sharedFloatCol ) ;
2667
2668 // fill TIME and INTERVAL
2669 Double midTime ;
2670 Double interval ;
2671 Vector<Double> intervalArr( nrow, 0.0 ) ;
2672 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2673 while( !iter.pastEnd() ) {
2674 //Table tab( iter.table() ) ;
2675 Table tab = iter.table() ;
2676
2677 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2678 uInt wid = widCol( 0 ) ;
2679
2680 getValidTimeRange( midTime, interval, tab ) ;
2681 mswCols.time().put( wid, midTime ) ;
2682 intervalArr[wid] = interval ;
2683
2684 iter.next() ;
2685 }
2686 mswCols.interval().putColumn( intervalArr ) ;
2687
2688// double endSec = mathutil::gettimeofday_sec() ;
2689// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2690}
2691
2692void MSWriter::fillSysCal()
2693{
2694 Table mssc = mstable_->sysCal() ;
2695
2696 {
2697 static const char *cols[] = {
2698 "BEAMNO", "IFNO", "TIME", "POLNO",
2699 NULL
2700 };
2701 static const TypeManagerImpl<uInt> tmUInt;
2702 static const TypeManagerImpl<Double> tmDouble;
2703 static const TypeManager *const tms[] = {
2704 &tmUInt, &tmUInt, &tmDouble, &tmUInt, NULL
2705 };
2706 //double t0 = mathutil::gettimeofday_sec() ;
2707 MSSysCalVisitor myVisitor(table_->table(),mssc);
2708 //double t1 = mathutil::gettimeofday_sec() ;
2709 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
2710 traverseTable(table_->table(), cols, tms, &myVisitor);
2711 //double t3 = mathutil::gettimeofday_sec() ;
2712 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
2713 }
2714
2715}
2716
2717void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
2718{
2719// double startSec = mathutil::gettimeofday_sec() ;
2720// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2721
2722 // sort table
2723 //Table stab = tab.sort( "TIME" ) ;
2724
2725 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2726 Vector<Double> timeArr = timeCol.getColumn() ;
2727 Double minTime ;
2728 Double maxTime ;
2729 minMax( minTime, maxTime, timeArr ) ;
2730 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2731 // unit for TIME
2732 // Scantable: "d"
2733 // MS: "s"
2734 me = midTime ;
2735 interval = ( maxTime - minTime ) * 86400.0 ;
2736
2737// double endSec = mathutil::gettimeofday_sec() ;
2738// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2739}
2740
2741void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2742{
2743// double startSec = mathutil::gettimeofday_sec() ;
2744// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2745
2746 // sort table
2747 //Table stab = tab.sort( "TIME" ) ;
2748
2749 Double minTime ;
2750 Double maxTime ;
2751 minMax( minTime, maxTime, atime ) ;
2752 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2753 // unit for TIME
2754 // Scantable: "d"
2755 // MS: "s"
2756 me = midTime ;
2757 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
2758
2759// double endSec = mathutil::gettimeofday_sec() ;
2760// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2761}
2762
2763}
Note: See TracBrowser for help on using the repository browser.