source: trunk/src/MSWriter.cpp@ 3004

Last change on this file since 3004 was 2975, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-5545, CAS-6597

Ready for Test: Yes

Interface Changes: Yes/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...

Fix for a bug in row flag handling.


File size: 79.3 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( idx < 4 && !isFilled[idx] ) ++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 < 4; ++ipol) {
353 if (isFilled[ipol] and 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 < 4; ++ipol) {
399 if (isFilled[ipol] and 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 addScanrate_(False)
679 {
680 rowidx = 0 ;
681 fieldName = "" ;
682 defaultFieldId = 0 ;
683 spwId = -1 ;
684 subscan = 1 ;
685 ptName = "" ;
686 srcId = 0 ;
687
688 row = TableRow( ms ) ;
689
690 initPolarization() ;
691 initFrequencies() ;
692
693 //
694 // add rows to MS
695 //
696 uInt addrow = table.nrow() ;
697 ms.addRow( addrow ) ;
698
699 // attach to Scantable columns
700 spectraCol.attach( table, "SPECTRA" ) ;
701 flagtraCol.attach( table, "FLAGTRA" ) ;
702 flagRowCol.attach( table, "FLAGROW" ) ;
703 tcalIdCol.attach( table, "TCAL_ID" ) ;
704 intervalCol.attach( table, "INTERVAL" ) ;
705 directionCol.attach( table, "DIRECTION" ) ;
706 scanRateCol.attach( table, "SCANRATE" ) ;
707 timeCol.attach( table, "TIME" ) ;
708 freqIdCol.attach( table, "FREQ_ID" ) ;
709 sourceNameCol.attach( table, "SRCNAME" ) ;
710 sourceDirectionCol.attach( table, "SRCDIRECTION" ) ;
711 fieldNameCol.attach( table, "FIELDNAME" ) ;
712
713 // MS subtables
714 attachSubtables() ;
715
716 // attach to MS columns
717 attachMain() ;
718 attachPointing() ;
719
720 // decide whether put SCANRATE to MS POINTING column
721 if (anyNE(scanRateCol.getColumn(), (Double)0.0)) {
722 addScanrate_ = True;
723 }
724 }
725
726 virtual void enterFieldName(const uInt recordNo, const String &/*columnValue*/) {
727 //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str());
728 fieldName = fieldNameCol.asString( recordNo ) ;
729 String::size_type pos = fieldName.find( "__" ) ;
730 if ( pos != String::npos ) {
731 fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
732 fieldName = fieldName.substr( 0, pos ) ;
733 }
734 else {
735 fieldId = defaultFieldId ;
736 defaultFieldId++ ;
737 }
738 Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ;
739 Vector<Double> srcDir = sourceDirectionCol( recordNo ) ;
740 Vector<Double> srate = scanRateCol( recordNo ) ;
741 String srcName = sourceNameCol.asString( recordNo ) ;
742
743 addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ;
744
745 // put value
746 *fieldIdRF = fieldId ;
747 }
748 virtual void leaveFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
749 }
750 virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue) {
751 //printf("%u: BeamNo: %u\n", recordNo, columnValue);
752
753 feedId = (Int)columnValue ;
754
755 // put value
756 *feed1RF = feedId ;
757 *feed2RF = feedId ;
758 }
759 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
760 }
761 virtual void enterScanNo(const uInt /*recordNo*/, uInt columnValue) {
762 //printf("%u: ScanNo: %u\n", recordNo, columnValue);
763
764 // put value
765 // CAS-5841: SCANNO should be consistent with MS SCAN_NUMBER
766 *scanNumberRF = (Int)columnValue ;
767 }
768 virtual void leaveScanNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
769 subscan = 1 ;
770 }
771 virtual void enterIfNo(const uInt recordNo, uInt columnValue) {
772 //printf("%u: IfNo: %u\n", recordNo, columnValue);
773
774 spwId = (Int)columnValue ;
775 uInt freqId = freqIdCol.asuInt( recordNo ) ;
776
777 Vector<Float> sp = spectraCol( recordNo ) ;
778 uInt nchan = sp.nelements() ;
779 holder->setNchan( nchan ) ;
780
781 addSpectralWindow( spwId, freqId ) ;
782
783 addFeed( feedId, spwId ) ;
784 }
785 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
786 }
787 virtual void enterSrcType(const uInt /*recordNo*/, Int columnValue) {
788 //printf("%u: SrcType: %d\n", recordNo, columnValue);
789
790 Int stateId = addState( columnValue ) ;
791
792 // put value
793 *stateIdRF = stateId ;
794 }
795 virtual void leaveSrcType(const uInt /*recordNo*/, Int /*columnValue*/) {
796 }
797 virtual void enterCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
798 //printf("%u: CycleNo: %u\n", recordNo, columnValue);
799 }
800 virtual void leaveCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
801 }
802 virtual void enterTime(const uInt recordNo, Double columnValue) {
803 //printf("%u: Time: %f\n", recordNo, columnValue);
804
805 Double timeSec = columnValue * 86400.0 ;
806 Double interval = intervalCol.asdouble( recordNo ) ;
807
808 if ( ptName.empty() ) {
809 Vector<Double> dir = directionCol( recordNo ) ;
810 Vector<Double> rate = scanRateCol( recordNo ) ;
811 //if ( anyNE( rate, 0.0 ) ) {
812 if (addScanrate_) {
813 Matrix<Double> msdir( 2, 2 ) ;
814 msdir.column( 0 ) = dir ;
815 msdir.column( 1 ) = rate ;
816 addPointing( timeSec, interval, msdir ) ;
817 }
818 else {
819 Matrix<Double> msdir( 2, 1 ) ;
820 msdir.column( 0 ) = dir ;
821 addPointing( timeSec, interval, msdir ) ;
822 }
823 }
824
825 // put value
826 *timeRF = timeSec ;
827 *timeCentroidRF = timeSec ;
828 *intervalRF = interval ;
829 *exposureRF = interval ;
830 }
831 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) {
832 if ( holder->nPol() > 0 ) {
833 Int polId = addPolarization() ;
834 Int ddId = addDataDescription( polId, spwId ) ;
835
836 // put field
837 *dataDescIdRF = ddId ;
838 holder->post() ;
839
840 // commit row
841 row.put( rowidx ) ;
842 rowidx++ ;
843
844 // reset holder
845 holder->reset() ;
846 }
847 }
848 virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
849 //printf("%u: PolNo: %d\n", recordNo, columnValue);
850 }
851 virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
852 }
853
854 virtual Bool visitRecord(const uInt recordNo,
855 const String &/*fieldName*/,
856 const uInt /*beamNo*/,
857 const uInt /*scanNo*/,
858 const uInt /*ifNo*/,
859 const Int /*srcType*/,
860 const uInt /*cycleNo*/,
861 const Double /*time*/,
862 const uInt polNo) {
863 //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo,
864 // fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo);
865
866 Vector<Float> sp = spectraCol( recordNo ) ;
867 Vector<uChar> tmp = flagtraCol( recordNo ) ;
868 Vector<Bool> fl( tmp.shape() ) ;
869 convertArray( fl, tmp ) ;
870 Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ;
871 holder->accumulate( polNo, sp, fl, flr ) ;
872
873 return True ;
874 }
875
876 virtual void finish() {
877 BaseMSWriterVisitor::finish();
878 //printf("Total: %u\n", count);
879
880 // remove rows
881 if ( ms.nrow() > rowidx ) {
882 uInt numRemove = ms.nrow() - rowidx ;
883 //cout << "numRemove = " << numRemove << endl ;
884 Vector<uInt> rows( numRemove ) ;
885 indgen( rows, rowidx ) ;
886 ms.removeRow( rows ) ;
887 }
888
889 // fill empty SPECTRAL_WINDOW rows
890 infillSpectralWindow() ;
891
892 // fill empty FIELD rows
893 infillField() ;
894 }
895
896 void dataColumnName( String name )
897 {
898 if ( name == "DATA" )
899 holder = new ComplexDataHolder( row, poltype ) ;
900 else if ( name == "FLOAT_DATA" )
901 holder = new FloatDataHolder( row, poltype ) ;
902 }
903 void pointingTableName( String name ) {
904 ptName = name ;
905 }
906 void setSourceRecord( Record &r ) {
907 srcRec = r ;
908 }
909private:
910 void addField( Int &fid, String &fname, String &srcName,
911 Vector<Double> &sdir, Vector<Double> &srate,
912 Double &tSec )
913 {
914 uInt nrow = fieldtab.nrow() ;
915 while( (Int)nrow <= fid ) {
916 fieldtab.addRow( 1, True ) ;
917 nrow++ ;
918 }
919
920 Matrix<Double> dir ;
921 Int numPoly = 0 ;
922 if ( anyNE( srate, 0.0 ) ) {
923 dir.resize( 2, 2 ) ;
924 dir.column( 0 ) = sdir ;
925 dir.column( 1 ) = srate ;
926 numPoly = 1 ;
927 }
928 else {
929 dir.resize( 2, 1 ) ;
930 dir.column( 0 ) = sdir ;
931 }
932 srcId = srcRec.asInt( srcName ) ;
933
934 TableRow tr( fieldtab ) ;
935 TableRecord &r = tr.record() ;
936 putField( "NAME", r, fname ) ;
937 putField( "NUM_POLY", r, numPoly ) ;
938 putField( "TIME", r, tSec ) ;
939 putField( "SOURCE_ID", r, srcId ) ;
940 defineField( "DELAY_DIR", r, dir ) ;
941 defineField( "REFERENCE_DIR", r, dir ) ;
942 defineField( "PHASE_DIR", r, dir ) ;
943 tr.put( fid ) ;
944
945 // for POINTING table
946 *poNameRF = fname ;
947 }
948 Int addState( Int &id )
949 {
950 String obsMode ;
951 Bool isSignal ;
952 Double tnoise ;
953 Double tload ;
954 queryType( id, obsMode, isSignal, tnoise, tload ) ;
955
956 String key = obsMode+"_"+String::toString( subscan ) ;
957 Int idx = -1 ;
958 uInt nEntry = stateEntry.nelements() ;
959 for ( uInt i = 0 ; i < nEntry ; i++ ) {
960 if ( stateEntry[i] == key ) {
961 idx = i ;
962 break ;
963 }
964 }
965 if ( idx == -1 ) {
966 uInt nrow = statetab.nrow() ;
967 statetab.addRow( 1, True ) ;
968 TableRow tr( statetab ) ;
969 TableRecord &r = tr.record() ;
970 putField( "OBS_MODE", r, obsMode ) ;
971 putField( "SIG", r, isSignal ) ;
972 isSignal = !isSignal ;
973 putField( "REF", r, isSignal ) ;
974 putField( "CAL", r, tnoise ) ;
975 putField( "LOAD", r, tload ) ;
976 tr.put( nrow ) ;
977 idx = nrow ;
978
979 stateEntry.resize( nEntry+1, True ) ;
980 stateEntry[nEntry] = key ;
981 }
982 subscan++ ;
983
984 return idx ;
985 }
986 void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir )
987 {
988 uInt nrow = potab.nrow() ;
989 potab.addRow() ;
990
991 *poNumPolyRF = dir.ncolumn() - 1 ;
992 *poTimeRF = tSec ;
993 *poTimeOriginRF = tSec ;
994 *poIntervalRF = interval ;
995 poDirectionRF.define( dir ) ;
996 poTargetRF.define( dir ) ;
997 porow.put( nrow ) ;
998 }
999 Int addPolarization()
1000 {
1001 Int idx = -1 ;
1002 Vector<Int> corrType = holder->corrTypeInt() ;
1003 uInt nEntry = polEntry.size() ;
1004 for ( uInt i = 0 ; i < nEntry ; i++ ) {
1005 if ( polEntry[i].conform( corrType ) && allEQ( polEntry[i], corrType ) ) {
1006 idx = i ;
1007 break ;
1008 }
1009 }
1010
1011 Int numCorr = holder->nPol() ;
1012 Matrix<Int> corrProduct = corrProductTemplate[numCorr].copy() ;
1013 if ( numCorr == 1 && (corrType[0] == Stokes::YY || corrType[0] == Stokes::LL ) ) {
1014 corrProduct = 1;
1015 }
1016
1017 if ( idx == -1 ) {
1018 uInt nrow = poltab.nrow() ;
1019 poltab.addRow( 1, True ) ;
1020 TableRow tr( poltab ) ;
1021 TableRecord &r = tr.record() ;
1022 putField( "NUM_CORR", r, numCorr ) ;
1023 defineField( "CORR_TYPE", r, corrType ) ;
1024 defineField( "CORR_PRODUCT", r, corrProduct ) ;
1025 tr.put( nrow ) ;
1026 idx = nrow ;
1027
1028 polEntry.resize( nEntry+1 ) ;
1029 polEntry[nEntry] = corrType ;
1030 }
1031
1032 return idx ;
1033 }
1034 Int addDataDescription( Int pid, Int sid )
1035 {
1036 Int idx = -1 ;
1037 uInt nItem = 2 ;
1038 uInt len = ddEntry.nelements() ;
1039 uInt nEntry = len / nItem ;
1040 const Int *dd_p = ddEntry.storage() ;
1041 for ( uInt i = 0 ; i < nEntry ; i++ ) {
1042 Int pol = *dd_p ;
1043 dd_p++ ;
1044 Int spw = *dd_p ;
1045 dd_p++ ;
1046 if ( pid == pol && sid == spw ) {
1047 idx = i ;
1048 break ;
1049 }
1050 }
1051
1052 if ( idx == -1 ) {
1053 uInt nrow = ddtab.nrow() ;
1054 ddtab.addRow( 1, True ) ;
1055 TableRow tr( ddtab ) ;
1056 TableRecord &r = tr.record() ;
1057 putField( "POLARIZATION_ID", r, pid ) ;
1058 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
1059 tr.put( nrow ) ;
1060 idx = nrow ;
1061
1062 ddEntry.resize( len+nItem ) ;
1063 ddEntry[len] = pid ;
1064 ddEntry[len+1] = sid ;
1065 }
1066
1067 return idx ;
1068 }
1069 void infillSpectralWindow()
1070 {
1071 ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ;
1072 Vector<Int> nchan = nchanCol.getColumn() ;
1073 TableRow tr( spwtab ) ;
1074 TableRecord &r = tr.record() ;
1075 Int mfr = freqframe ;
1076 Int oneChan = 1 ;
1077 Vector<Double> dummy( 1, 0.0 ) ;
1078 putField( "MEAS_FREQ_REF", r, mfr ) ;
1079 defineField( "CHAN_FREQ", r, dummy ) ;
1080 defineField( "CHAN_WIDTH", r, dummy ) ;
1081 defineField( "EFFECTIVE_BW", r, dummy ) ;
1082 defineField( "RESOLUTION", r, dummy ) ;
1083 putField( "NUM_CHAN", r, oneChan ) ;
1084
1085 for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
1086 if ( nchan[i] == 0 )
1087 tr.put( i ) ;
1088 }
1089 }
1090 void infillField()
1091 {
1092 ScalarColumn<Int> sourceIdCol(fieldtab, "SOURCE_ID");
1093 ArrayColumn<Double> delayDirCol(fieldtab, "DELAY_DIR");
1094 ArrayColumn<Double> phaseDirCol(fieldtab, "PHASE_DIR");
1095 ArrayColumn<Double> referenceDirCol(fieldtab, "REFERENCE_DIR");
1096 uInt nrow = fieldtab.nrow();
1097 Matrix<Double> dummy(IPosition(2, 2, 1), 0.0);
1098 for (uInt irow = 0; irow < nrow; ++irow) {
1099 if (!phaseDirCol.isDefined(irow)) {
1100 delayDirCol.put(irow, dummy);
1101 phaseDirCol.put(irow, dummy);
1102 referenceDirCol.put(irow, dummy);
1103 sourceIdCol.put(irow, -1);
1104 }
1105 }
1106 }
1107 void addSpectralWindow( Int sid, uInt fid )
1108 {
1109 if (processedFreqId.find((uInt)fid) == processedFreqId.end()
1110 || processedIFNO.find((uInt)sid) == processedIFNO.end() ) {
1111 uInt nrow = spwtab.nrow() ;
1112 while( (Int)nrow <= sid ) {
1113 spwtab.addRow( 1, True ) ;
1114 nrow++ ;
1115 }
1116 processedFreqId.insert((uInt)fid);
1117 processedIFNO.insert((uInt)sid);
1118 }
1119 else {
1120 return ;
1121 }
1122
1123
1124 Double rp = refpix[fid] ;
1125 Double rv = refval[fid] ;
1126 Double ic = increment[fid] ;
1127
1128 Int mfrInt = (Int)freqframe ;
1129 Int nchan = holder->nChan() ;
1130 Double bw = nchan * abs( ic ) ;
1131 Double reffreq = rv - rp * ic ;
1132 Int netsb = 0 ; // USB->0, LSB->1
1133 if ( ic < 0 )
1134 netsb = 1 ;
1135 Vector<Double> res( nchan, abs(ic) ) ;
1136 Vector<Double> cw( nchan, ic ) ;
1137 Vector<Double> chanf( nchan ) ;
1138 indgen( chanf, reffreq, ic ) ;
1139
1140 TableRow tr( spwtab ) ;
1141 TableRecord &r = tr.record() ;
1142 putField( "MEAS_FREQ_REF", r, mfrInt ) ;
1143 putField( "NUM_CHAN", r, nchan ) ;
1144 putField( "TOTAL_BANDWIDTH", r, bw ) ;
1145 putField( "REF_FREQUENCY", r, reffreq ) ;
1146 putField( "NET_SIDEBAND", r, netsb ) ;
1147 defineField( "RESOLUTION", r, res ) ;
1148// defineField( "CHAN_WIDTH", r, res ) ;
1149 defineField( "CHAN_WIDTH", r, cw ) ;
1150 defineField( "EFFECTIVE_BW", r, res ) ;
1151 defineField( "CHAN_FREQ", r, chanf ) ;
1152 tr.put( sid ) ;
1153 }
1154 void addFeed( Int fid, Int sid )
1155 {
1156 Int idx = -1 ;
1157 uInt nItem = 2 ;
1158 uInt len = feedEntry.nelements() ;
1159 uInt nEntry = len / nItem ;
1160 const Int *fe_p = feedEntry.storage() ;
1161 for ( uInt i = 0 ; i < nEntry ; i++ ) {
1162 Int feed = *fe_p ;
1163 fe_p++ ;
1164 Int spw = *fe_p ;
1165 fe_p++ ;
1166 if ( fid == feed && sid == spw ) {
1167 idx = i ;
1168 break ;
1169 }
1170 }
1171
1172
1173 if ( idx == -1 ) {
1174 uInt nrow = feedtab.nrow() ;
1175 feedtab.addRow( 1, True ) ;
1176 Int numReceptors = 2 ;
1177 Vector<String> polType( numReceptors ) ;
1178 Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ;
1179 Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
1180 if ( poltype == "linear" ) {
1181 polType[0] = "X" ;
1182 polType[1] = "Y" ;
1183 }
1184 else if ( poltype == "circular" ) {
1185 polType[0] = "R" ;
1186 polType[1] = "L" ;
1187 }
1188 else {
1189 polType[0] = "X" ;
1190 polType[1] = "Y" ;
1191 }
1192 Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
1193
1194 TableRow tr( feedtab ) ;
1195 TableRecord &r = tr.record() ;
1196 putField( "FEED_ID", r, fid ) ;
1197 putField( "BEAM_ID", r, fid ) ;
1198 Int tmp = 0 ;
1199 putField( "ANTENNA_ID", r, tmp ) ;
1200 putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
1201 putField( "NUM_RECEPTORS", r, numReceptors ) ;
1202 defineField( "POLARIZATION_TYPE", r, polType ) ;
1203 defineField( "BEAM_OFFSET", r, beamOffset ) ;
1204 defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ;
1205 defineField( "POL_RESPONSE", r, polResponse ) ;
1206 tr.put( nrow ) ;
1207
1208 feedEntry.resize( len+nItem ) ;
1209 feedEntry[len] = fid ;
1210 feedEntry[len+1] = sid ;
1211 }
1212 }
1213 void initPolarization()
1214 {
1215 const TableRecord &keys = table.keywordSet() ;
1216 poltype = keys.asString( "POLTYPE" ) ;
1217
1218 initCorrProductTemplate() ;
1219 }
1220 void initFrequencies()
1221 {
1222 const TableRecord &keys = table.keywordSet() ;
1223 Table tab = keys.asTable( "FREQUENCIES" ) ;
1224 ROScalarColumn<uInt> idcol( tab, "ID" ) ;
1225 ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ;
1226 ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ;
1227 ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ;
1228 Vector<uInt> id = idcol.getColumn() ;
1229 Vector<Double> rp = rpcol.getColumn() ;
1230 Vector<Double> rv = rvcol.getColumn() ;
1231 Vector<Double> ic = iccol.getColumn() ;
1232 for ( uInt i = 0 ; i < id.nelements() ; i++ ) {
1233 refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ;
1234 refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ;
1235 increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ;
1236 }
1237 String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ;
1238 MFrequency::getType( freqframe, frameStr ) ;
1239 }
1240 void attachSubtables()
1241 {
1242 //const TableRecord &keys = table.keywordSet() ;
1243 TableRecord &mskeys = ms.rwKeywordSet() ;
1244
1245 // FIELD table
1246 fieldtab = mskeys.asTable( "FIELD" ) ;
1247
1248 // SPECTRAL_WINDOW table
1249 spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ;
1250
1251 // POINTING table
1252 potab = mskeys.asTable( "POINTING" ) ;
1253
1254 // POLARIZATION table
1255 poltab = mskeys.asTable( "POLARIZATION" ) ;
1256
1257 // DATA_DESCRIPTION table
1258 ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ;
1259
1260 // STATE table
1261 statetab = mskeys.asTable( "STATE" ) ;
1262
1263 // FEED table
1264 feedtab = mskeys.asTable( "FEED" ) ;
1265 }
1266 void attachMain()
1267 {
1268 TableRecord &r = row.record() ;
1269 dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ;
1270 timeRF.attachToRecord( r, "TIME" ) ;
1271 timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ;
1272 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1273 exposureRF.attachToRecord( r, "EXPOSURE" ) ;
1274 fieldIdRF.attachToRecord( r, "FIELD_ID" ) ;
1275 feed1RF.attachToRecord( r, "FEED1" ) ;
1276 feed2RF.attachToRecord( r, "FEED2" ) ;
1277 scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ;
1278 stateIdRF.attachToRecord( r, "STATE_ID" ) ;
1279
1280 // constant values
1281 //Int id = 0 ;
1282 RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ;
1283 *intRF = 0 ;
1284 intRF.attachToRecord( r, "ANTENNA1" ) ;
1285 *intRF = 0 ;
1286 intRF.attachToRecord( r, "ANTENNA2" ) ;
1287 *intRF = 0 ;
1288 intRF.attachToRecord( r, "ARRAY_ID" ) ;
1289 *intRF = 0 ;
1290 intRF.attachToRecord( r, "PROCESSOR_ID" ) ;
1291 *intRF = 0 ;
1292 RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ;
1293 arrayRF.define( Vector<Double>( 3, 0.0 ) ) ;
1294 }
1295 void attachPointing()
1296 {
1297 porow = TableRow( potab ) ;
1298 TableRecord &r = porow.record() ;
1299 poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ;
1300 poTimeRF.attachToRecord( r, "TIME" ) ;
1301 poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ;
1302 poIntervalRF.attachToRecord( r, "INTERVAL" ) ;
1303 poNameRF.attachToRecord( r, "NAME" ) ;
1304 poDirectionRF.attachToRecord( r, "DIRECTION" ) ;
1305 poTargetRF.attachToRecord( r, "TARGET" ) ;
1306
1307 // constant values
1308 RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ;
1309 *antIdRF = 0 ;
1310 RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ;
1311 *trackingRF = True ;
1312 }
1313 void queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1314 {
1315 t = 0.0 ;
1316 l = 0.0 ;
1317
1318 String sep1="#" ;
1319 String sep2="," ;
1320 String target="OBSERVE_TARGET" ;
1321 String atmcal="CALIBRATE_TEMPERATURE" ;
1322 String onstr="ON_SOURCE" ;
1323 String offstr="OFF_SOURCE" ;
1324 String pswitch="POSITION_SWITCH" ;
1325 String nod="NOD" ;
1326 String fswitch="FREQUENCY_SWITCH" ;
1327 String sigstr="SIG" ;
1328 String refstr="REF" ;
1329 String unspecified="UNSPECIFIED" ;
1330 String ftlow="LOWER" ;
1331 String fthigh="HIGHER" ;
1332 switch ( type ) {
1333 case SrcType::PSON:
1334 stype = target+sep1+onstr+sep2+pswitch ;
1335 b = True ;
1336 break ;
1337 case SrcType::PSOFF:
1338 stype = target+sep1+offstr+sep2+pswitch ;
1339 b = False ;
1340 break ;
1341 case SrcType::NOD:
1342 stype = target+sep1+onstr+sep2+nod ;
1343 b = True ;
1344 break ;
1345 case SrcType::FSON:
1346 stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1347 b = True ;
1348 break ;
1349 case SrcType::FSOFF:
1350 stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
1351 b = False ;
1352 break ;
1353 case SrcType::SKY:
1354 stype = atmcal+sep1+offstr+sep2+unspecified ;
1355 b = False ;
1356 break ;
1357 case SrcType::HOT:
1358 stype = atmcal+sep1+offstr+sep2+unspecified ;
1359 b = False ;
1360 break ;
1361 case SrcType::WARM:
1362 stype = atmcal+sep1+offstr+sep2+unspecified ;
1363 b = False ;
1364 break ;
1365 case SrcType::COLD:
1366 stype = atmcal+sep1+offstr+sep2+unspecified ;
1367 b = False ;
1368 break ;
1369 case SrcType::PONCAL:
1370 stype = atmcal+sep1+onstr+sep2+pswitch ;
1371 b = True ;
1372 break ;
1373 case SrcType::POFFCAL:
1374 stype = atmcal+sep1+offstr+sep2+pswitch ;
1375 b = False ;
1376 break ;
1377 case SrcType::NODCAL:
1378 stype = atmcal+sep1+onstr+sep2+nod ;
1379 b = True ;
1380 break ;
1381 case SrcType::FONCAL:
1382 stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
1383 b = True ;
1384 break ;
1385 case SrcType::FOFFCAL:
1386 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
1387 b = False ;
1388 break ;
1389 case SrcType::FSLO:
1390 stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
1391 b = True ;
1392 break ;
1393 case SrcType::FLOOFF:
1394 stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1395 b = False ;
1396 break ;
1397 case SrcType::FLOSKY:
1398 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1399 b = False ;
1400 break ;
1401 case SrcType::FLOHOT:
1402 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1403 b = False ;
1404 break ;
1405 case SrcType::FLOWARM:
1406 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1407 b = False ;
1408 break ;
1409 case SrcType::FLOCOLD:
1410 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
1411 b = False ;
1412 break ;
1413 case SrcType::FSHI:
1414 stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
1415 b = True ;
1416 break ;
1417 case SrcType::FHIOFF:
1418 stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1419 b = False ;
1420 break ;
1421 case SrcType::FHISKY:
1422 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1423 b = False ;
1424 break ;
1425 case SrcType::FHIHOT:
1426 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1427 b = False ;
1428 break ;
1429 case SrcType::FHIWARM:
1430 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1431 b = False ;
1432 break ;
1433 case SrcType::FHICOLD:
1434 stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
1435 b = False ;
1436 break ;
1437 case SrcType::SIG:
1438 stype = target+sep1+onstr+sep2+unspecified ;
1439 b = True ;
1440 break ;
1441 case SrcType::REF:
1442 stype = target+sep1+offstr+sep2+unspecified ;
1443 b = False ;
1444 break ;
1445 default:
1446 stype = unspecified ;
1447 b = True ;
1448 break ;
1449 }
1450 }
1451 void initCorrProductTemplate()
1452 {
1453 Int n = 1 ;
1454 {
1455 Matrix<Int> c( 2, n, 0 ) ;
1456 corrProductTemplate[n] = c ;
1457 }
1458 n = 2 ;
1459 {
1460 Matrix<Int> c( 2, n, 0 ) ;
1461 c.column( 1 ) = 1 ;
1462 corrProductTemplate[n] = c ;
1463 }
1464 n = 4 ;
1465 {
1466 Matrix<Int> c( 2, n, 0 ) ;
1467 c( 0, 2 ) = 1 ;
1468 c( 0, 3 ) = 1 ;
1469 c( 1, 1 ) = 1 ;
1470 c( 1, 3 ) = 1 ;
1471 corrProductTemplate[n] = c ;
1472 }
1473 }
1474
1475 Table &ms;
1476 TableRow row;
1477 uInt rowidx;
1478 String fieldName;
1479 Int fieldId;
1480 Int srcId;
1481 Int defaultFieldId;
1482 Int spwId;
1483 Int feedId;
1484 Int subscan;
1485 CountedPtr<DataHolder> holder;
1486 String ptName;
1487 Bool useFloat;
1488 String poltype;
1489 Bool addScanrate_;
1490
1491 // MS subtables
1492 Table spwtab;
1493 Table statetab;
1494 Table ddtab;
1495 Table poltab;
1496 Table fieldtab;
1497 Table feedtab;
1498 Table potab;
1499
1500 // Scantable MAIN columns
1501 ROArrayColumn<Float> spectraCol;
1502 ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol;
1503 ROArrayColumn<uChar> flagtraCol;
1504 ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol,
1505 sourceNameCol,fieldNameCol;
1506
1507 // MS MAIN columns
1508 RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF,
1509 scanNumberRF,stateIdRF;
1510 RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF;
1511
1512 // MS POINTING columns
1513 TableRow porow;
1514 RecordFieldPtr<Int> poNumPolyRF ;
1515 RecordFieldPtr<Double> poTimeRF,
1516 poTimeOriginRF,
1517 poIntervalRF ;
1518 RecordFieldPtr<String> poNameRF ;
1519 RecordFieldPtr< Matrix<Double> > poDirectionRF,
1520 poTargetRF ;
1521
1522 Vector<String> stateEntry;
1523 Block<Int> ddEntry;
1524 Block<Int> feedEntry;
1525 vector< Vector<Int> > polEntry;
1526 set<uInt> processedFreqId;
1527 set<uInt> processedIFNO;
1528 map<uInt,Double> refpix;
1529 map<uInt,Double> refval;
1530 map<uInt,Double> increment;
1531 MFrequency::Types freqframe;
1532 Record srcRec;
1533 map< Int, Matrix<Int> > corrProductTemplate;
1534};
1535
1536class BaseMSSysCalVisitor: public TableVisitor {
1537 uInt lastRecordNo;
1538 uInt lastBeamNo, lastIfNo, lastPolNo;
1539 Double lastTime;
1540protected:
1541 const Table &table;
1542 uInt count;
1543public:
1544 BaseMSSysCalVisitor(const Table &table)
1545 : table(table)
1546 {
1547 count = 0;
1548 }
1549
1550 virtual void enterBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1551 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1552 virtual void enterIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1553 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1554 virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1555 virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) { }
1556 virtual void enterTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1557 virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) { }
1558
1559 virtual Bool visitRecord(const uInt /*recordNo*/,
1560 const uInt /*beamNo*/,
1561 const uInt /*ifNo*/,
1562 const uInt /*polNo*/,
1563 const Double /*time*/) { return True ;}
1564
1565 virtual Bool visit(Bool isFirst, const uInt recordNo,
1566 const uInt nCols, void const *const colValues[]) {
1567 uInt beamNo, ifNo, polNo;
1568 Double time;
1569 { // prologue
1570 uInt i = 0;
1571 {
1572 const uInt *col = (const uInt *)colValues[i++];
1573 beamNo = col[recordNo];
1574 }
1575 {
1576 const uInt *col = (const uInt *)colValues[i++];
1577 ifNo = col[recordNo];
1578 }
1579 {
1580 const Double *col = (const Double *)colValues[i++];
1581 time = col[recordNo];
1582 }
1583 {
1584 const uInt *col = (const uInt *)colValues[i++];
1585 polNo = col[recordNo];
1586 }
1587 assert(nCols == i);
1588 }
1589
1590 if (isFirst) {
1591 enterBeamNo(recordNo, beamNo);
1592 enterIfNo(recordNo, ifNo);
1593 enterTime(recordNo, time);
1594 enterPolNo(recordNo, polNo);
1595 } else {
1596 if (lastBeamNo != beamNo) {
1597 leavePolNo(lastRecordNo, lastPolNo);
1598 leaveTime(lastRecordNo, lastTime);
1599 leaveIfNo(lastRecordNo, lastIfNo);
1600 leaveBeamNo(lastRecordNo, lastBeamNo);
1601
1602 enterBeamNo(recordNo, beamNo);
1603 enterIfNo(recordNo, ifNo);
1604 enterTime(recordNo, time);
1605 enterPolNo(recordNo, polNo);
1606 } else if (lastIfNo != ifNo) {
1607 leavePolNo(lastRecordNo, lastPolNo);
1608 leaveTime(lastRecordNo, lastTime);
1609 leaveIfNo(lastRecordNo, lastIfNo);
1610
1611 enterIfNo(recordNo, ifNo);
1612 enterTime(recordNo, time);
1613 enterPolNo(recordNo, polNo);
1614 } else if (lastTime != time) {
1615 leavePolNo(lastRecordNo, lastPolNo);
1616 leaveTime(lastRecordNo, lastTime);
1617
1618 enterTime(recordNo, time);
1619 enterPolNo(recordNo, polNo);
1620 } else if (lastPolNo != polNo) {
1621 leavePolNo(lastRecordNo, lastPolNo);
1622 enterPolNo(recordNo, polNo);
1623 }
1624 }
1625 count++;
1626 Bool result = visitRecord(recordNo, beamNo, ifNo, polNo, time);
1627
1628 { // epilogue
1629 lastRecordNo = recordNo;
1630
1631 lastBeamNo = beamNo;
1632 lastIfNo = ifNo;
1633 lastPolNo = polNo;
1634 lastTime = time;
1635 }
1636 return result ;
1637 }
1638
1639 virtual void finish() {
1640 if (count > 0) {
1641 leavePolNo(lastRecordNo, lastPolNo);
1642 leaveTime(lastRecordNo, lastTime);
1643 leaveIfNo(lastRecordNo, lastIfNo);
1644 leaveBeamNo(lastRecordNo, lastBeamNo);
1645 }
1646 }
1647};
1648
1649class BaseTsysHolder
1650{
1651public:
1652 BaseTsysHolder( ROArrayColumn<Float> &tsysCol )
1653 : col( tsysCol ),
1654 nchan(0)
1655 {
1656 reset() ;
1657 }
1658 virtual ~BaseTsysHolder() {}
1659 virtual Array<Float> getTsys() = 0 ;
1660 void setNchan( uInt n ) { nchan = n ; }
1661 void appendTsys( uInt row )
1662 {
1663 Vector<Float> v = col( row ) ;
1664 uInt len = tsys.nrow() ;
1665 tsys.resize( len+1, nchan, True ) ;
1666 if ( v.nelements() == nchan )
1667 tsys.row( len ) = v ;
1668 else
1669 tsys.row( len ) = v[0] ;
1670 }
1671 void setTsys( uInt row, uInt idx )
1672 {
1673 if ( idx >= nrow() )
1674 appendTsys( row ) ;
1675 else {
1676 Vector<Float> v = col( row ) ;
1677 if ( v.nelements() == nchan )
1678 tsys.row( idx ) = v ;
1679 else
1680 tsys.row( idx ) = v[0] ;
1681 }
1682 }
1683 void reset()
1684 {
1685 tsys.resize() ;
1686 }
1687 uInt nrow() { return tsys.nrow() ; }
1688 Bool isEffective()
1689 {
1690 return ( !(tsys.empty()) && anyNE( tsys, (Float)1.0 ) ) ;
1691 }
1692 BaseTsysHolder &operator= ( const BaseTsysHolder &v )
1693 {
1694 if ( this != &v )
1695 tsys.assign( v.tsys ) ;
1696 return *this ;
1697 }
1698protected:
1699 ROArrayColumn<Float> col ;
1700 Matrix<Float> tsys ;
1701 uInt nchan ;
1702};
1703
1704class TsysHolder : public BaseTsysHolder
1705{
1706public:
1707 TsysHolder( ROArrayColumn<Float> &tsysCol )
1708 : BaseTsysHolder( tsysCol )
1709 {}
1710 virtual ~TsysHolder() {}
1711 virtual Array<Float> getTsys()
1712 {
1713 return tsys.column( 0 ) ;
1714 }
1715};
1716
1717class TsysSpectrumHolder : public BaseTsysHolder
1718{
1719public:
1720 TsysSpectrumHolder( ROArrayColumn<Float> &tsysCol )
1721 : BaseTsysHolder( tsysCol )
1722 {}
1723 virtual ~TsysSpectrumHolder() {}
1724 virtual Array<Float> getTsys()
1725 {
1726 return tsys ;
1727 }
1728};
1729
1730class BaseTcalProcessor
1731{
1732public:
1733 BaseTcalProcessor( ROArrayColumn<Float> &tcalCol )
1734 : col_( tcalCol )
1735 {}
1736 virtual ~BaseTcalProcessor() {}
1737 void setTcalId( Vector<uInt> &tcalId ) { id_.assign( tcalId ) ; }
1738 virtual Array<Float> getTcal() = 0 ;
1739protected:
1740 ROArrayColumn<Float> col_ ;
1741 Vector<uInt> id_ ;
1742};
1743
1744class TcalProcessor : public BaseTcalProcessor
1745{
1746public:
1747 TcalProcessor( ROArrayColumn<Float> &tcalCol )
1748 : BaseTcalProcessor( tcalCol )
1749 {}
1750 virtual ~TcalProcessor() {}
1751 virtual Array<Float> getTcal()
1752 {
1753 uInt npol = id_.nelements() ;
1754 Vector<Float> tcal( npol ) ;
1755 for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
1756 tcal[ipol] = col_( id_[ipol] ).data()[0] ;
1757 //cout << "TcalProcessor: tcal = " << tcal << endl ;
1758 return tcal ;
1759 }
1760};
1761
1762class TcalSpectrumProcessor : public BaseTcalProcessor
1763{
1764public:
1765 TcalSpectrumProcessor( ROArrayColumn<Float> &tcalCol )
1766 : BaseTcalProcessor( tcalCol )
1767 {}
1768 virtual ~TcalSpectrumProcessor() {}
1769 virtual Array<Float> getTcal()
1770 {
1771 uInt npol = id_.nelements() ;
1772 //Vector<Float> tcal0 = col_( 0 ) ;
1773 Vector<Float> tcal0 = col_( id_[0] ) ;
1774 uInt nchan = tcal0.nelements() ;
1775 Matrix<Float> tcal( npol, nchan ) ;
1776 tcal.row( 0 ) = tcal0 ;
1777 for ( uInt ipol = 1 ; ipol < npol ; ipol++ )
1778 tcal.row( ipol ) = col_( id_[ipol] ) ;
1779 return tcal ;
1780 }
1781};
1782
1783class MSSysCalVisitor : public BaseMSSysCalVisitor
1784{
1785public:
1786 MSSysCalVisitor( const Table &from, Table &to )
1787 : BaseMSSysCalVisitor( from ),
1788 sctab( to ),
1789 rowidx( 0 ),
1790 polno_()
1791 {
1792 scrow = TableRow( sctab ) ;
1793
1794 lastTcalId.resize() ;
1795 theTcalId.resize() ;
1796 startTime = 0.0 ;
1797 endTime = 0.0 ;
1798
1799 const TableRecord &keys = table.keywordSet() ;
1800 Table tcalTable = keys.asTable( "TCAL" ) ;
1801 tcalCol.attach( tcalTable, "TCAL" ) ;
1802 tsysCol.attach( table, "TSYS" ) ;
1803 tcalIdCol.attach( table, "TCAL_ID" ) ;
1804 intervalCol.attach( table, "INTERVAL" ) ;
1805 effectiveTcal.resize( tcalTable.nrow() ) ;
1806 for ( uInt irow = 0 ; irow < tcalTable.nrow() ; irow++ ) {
1807 if ( allEQ( tcalCol( irow ), (Float)1.0 ) )
1808 effectiveTcal[irow] = False ;
1809 else
1810 effectiveTcal[irow] = True ;
1811 }
1812
1813 TableRecord &r = scrow.record() ;
1814 RecordFieldPtr<Int> antennaIdRF( r, "ANTENNA_ID" ) ;
1815 *antennaIdRF = 0 ;
1816 feedIdRF.attachToRecord( r, "FEED_ID" ) ;
1817 specWinIdRF.attachToRecord( r, "SPECTRAL_WINDOW_ID" ) ;
1818 timeRF.attachToRecord( r, "TIME" ) ;
1819 intervalRF.attachToRecord( r, "INTERVAL" ) ;
1820 if ( r.isDefined( "TCAL" ) ) {
1821 tcalRF.attachToRecord( r, "TCAL" ) ;
1822 tcalProcessor = new TcalProcessor( tcalCol ) ;
1823 }
1824 else if ( r.isDefined( "TCAL_SPECTRUM" ) ) {
1825 tcalRF.attachToRecord( r, "TCAL_SPECTRUM" ) ;
1826 tcalProcessor = new TcalSpectrumProcessor( tcalCol ) ;
1827 }
1828 if ( r.isDefined( "TSYS" ) ) {
1829 tsysRF.attachToRecord( r, "TSYS" ) ;
1830 theTsys = new TsysHolder( tsysCol ) ;
1831 lastTsys = new TsysHolder( tsysCol ) ;
1832 }
1833 else {
1834 tsysRF.attachToRecord( r, "TSYS_SPECTRUM" ) ;
1835 theTsys = new TsysSpectrumHolder( tsysCol ) ;
1836 lastTsys = new TsysSpectrumHolder( tsysCol ) ;
1837 }
1838
1839 }
1840
1841 virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue)
1842 {
1843 *feedIdRF = (Int)columnValue ;
1844 }
1845 virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/)
1846 {
1847 }
1848 virtual void enterIfNo(const uInt recordNo, uInt columnValue)
1849 {
1850 //cout << "enterIfNo" << endl ;
1851 ROArrayColumn<Float> sp( table, "SPECTRA" ) ;
1852 uInt nchan = sp( recordNo ).nelements() ;
1853 theTsys->setNchan( nchan ) ;
1854 lastTsys->setNchan( nchan ) ;
1855
1856 *specWinIdRF = (Int)columnValue ;
1857 }
1858 virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/)
1859 {
1860 //cout << "leaveIfNo" << endl ;
1861 post() ;
1862 reset(true);
1863 startTime = 0.0 ;
1864 endTime = 0.0 ;
1865 }
1866 virtual void enterTime(const uInt recordNo, Double columnValue)
1867 {
1868 //cout << "enterTime" << endl ;
1869 interval = intervalCol.asdouble( recordNo ) ;
1870 // start time and end time
1871 if ( startTime == 0.0 ) {
1872 startTime = columnValue * 86400.0 - 0.5 * interval ;
1873 endTime = columnValue * 86400.0 + 0.5 * interval ;
1874 }
1875 }
1876 virtual void leaveTime(const uInt /*recordNo*/, Double columnValue)
1877 {
1878 //cout << "leaveTime" << endl ;
1879 if ( isUpdated() ) {
1880 post() ;
1881 *lastTsys = *theTsys ;
1882 lastTcalId = theTcalId ;
1883 reset(false);
1884 startTime = columnValue * 86400.0 - 0.5 * interval ;
1885 endTime = columnValue * 86400.0 + 0.5 * interval ;
1886 }
1887 else {
1888 endTime = columnValue * 86400.0 + 0.5 * interval ;
1889 }
1890 }
1891 virtual void enterPolNo(const uInt recordNo, uInt columnValue)
1892 {
1893 //cout << "enterPolNo" << endl ;
1894 Vector<Float> tsys = tsysCol( recordNo ) ;
1895 uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
1896 polno_.insert( columnValue ) ;
1897 uInt numPol = polno_.size() ;
1898 if ( lastTsys->nrow() < numPol )
1899 lastTsys->appendTsys( recordNo ) ;
1900 if ( lastTcalId.nelements() <= numPol )
1901 appendTcalId( lastTcalId, tcalId, numPol-1 ) ;
1902 if ( theTsys->nrow() < numPol )
1903 theTsys->appendTsys( recordNo ) ;
1904 else {
1905 theTsys->setTsys( recordNo, numPol-1 ) ;
1906 }
1907 if ( theTcalId.nelements() < numPol )
1908 appendTcalId( theTcalId, tcalId, numPol-1 ) ;
1909 else
1910 setTcalId( theTcalId, tcalId, numPol-1 ) ;
1911 }
1912 virtual void leavePolNo( const uInt /*recordNo*/, uInt /*columnValue*/ )
1913 {
1914 }
1915
1916private:
1917 void reset(bool completely)
1918 {
1919 if (completely) {
1920 lastTsys->reset() ;
1921 lastTcalId.resize() ;
1922 }
1923 theTsys->reset() ;
1924 theTcalId.resize() ;
1925 polno_.clear();
1926 }
1927 void appendTcalId( Vector<uInt> &v, uInt &elem, uInt polId )
1928 {
1929 v.resize( polId+1, True ) ;
1930 v[polId] = elem ;
1931 }
1932 void setTcalId( Vector<uInt> &v, uInt &elem, uInt polId )
1933 {
1934 v[polId] = elem ;
1935 }
1936 void post()
1937 {
1938 // check if given Tcal and Tsys is effective
1939 Bool isEffective = False ;
1940 for ( uInt ipol = 0 ; ipol < lastTcalId.nelements() ; ipol++ ) {
1941 if ( effectiveTcal[lastTcalId[ipol]] ) {
1942 isEffective = True ;
1943 break ;
1944 }
1945 }
1946 if ( !isEffective ) {
1947 if ( !(lastTsys->isEffective()) )
1948 return ;
1949 }
1950
1951 //cout << " interval: " << (endTime-startTime) << " lastTcalId = " << lastTcalId << endl ;
1952 Double midTime = 0.5 * ( startTime + endTime ) ;
1953 Double interval = endTime - startTime ;
1954 *timeRF = midTime ;
1955 *intervalRF = interval ;
1956 tcalProcessor->setTcalId( lastTcalId ) ;
1957 Array<Float> tcal = tcalProcessor->getTcal() ;
1958 Array<Float> tsys = lastTsys->getTsys() ;
1959 tcalRF.define( tcal ) ;
1960 tsysRF.define( tsys ) ;
1961 sctab.addRow( 1, True ) ;
1962 scrow.put( rowidx ) ;
1963 rowidx++ ;
1964 }
1965
1966 Bool isUpdated()
1967 {
1968 Bool ret = (anyNE( theTcalId, lastTcalId ) || anyNE( theTsys->getTsys(), lastTsys->getTsys() )) ;
1969 return ret ;
1970 }
1971
1972 Table &sctab;
1973 TableRow scrow;
1974 uInt rowidx;
1975
1976 Double startTime,endTime,interval;
1977
1978 CountedPtr<BaseTsysHolder> lastTsys,theTsys;
1979 Vector<uInt> lastTcalId,theTcalId;
1980 set<uInt> polno_;
1981 CountedPtr<BaseTcalProcessor> tcalProcessor ;
1982 Vector<Bool> effectiveTcal;
1983
1984 RecordFieldPtr<Int> feedIdRF,specWinIdRF;
1985 RecordFieldPtr<Double> timeRF,intervalRF;
1986 RecordFieldPtr< Array<Float> > tcalRF,tsysRF;
1987
1988 ROArrayColumn<Float> tsysCol,tcalCol;
1989 ROTableColumn tcalIdCol,intervalCol;
1990};
1991
1992MSWriter::MSWriter(CountedPtr<Scantable> stable)
1993 : table_(stable),
1994 mstable_(NULL),
1995 isWeather_(False),
1996 tcalSpec_(False),
1997 tsysSpec_(False),
1998 ptTabName_("")
1999{
2000 os_ = LogIO() ;
2001 os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
2002// os_ << "MSWriter::MSWriter()" << LogIO::POST ;
2003
2004 // initialize writer
2005 init() ;
2006}
2007
2008MSWriter::~MSWriter()
2009{
2010 os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
2011// os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
2012
2013 if ( mstable_ != 0 )
2014 delete mstable_ ;
2015}
2016
2017bool MSWriter::write(const string& filename, const Record& rec)
2018{
2019 os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
2020 //double startSec = mathutil::gettimeofday_sec() ;
2021 //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
2022
2023 filename_ = filename ;
2024
2025 // parsing MS options
2026 Bool overwrite = False ;
2027 if ( rec.isDefined( "ms" ) ) {
2028 Record msrec = rec.asRecord( "ms" ) ;
2029 if ( msrec.isDefined( "overwrite" ) ) {
2030 overwrite = msrec.asBool( "overwrite" ) ;
2031 }
2032 }
2033
2034 os_ << "Parsing MS options" << endl ;
2035 os_ << " overwrite = " << overwrite << LogIO::POST ;
2036
2037 File file( filename_ ) ;
2038 if ( file.exists() ) {
2039 if ( overwrite ) {
2040 os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
2041 if ( file.isRegular() ) RegularFile(file).remove() ;
2042 else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
2043 else SymLink(file).remove() ;
2044 }
2045 else {
2046 os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
2047 return False ;
2048 }
2049 }
2050
2051 // set up MS
2052 setupMS() ;
2053
2054 // subtables
2055 // OBSERVATION
2056 fillObservation() ;
2057
2058 // ANTENNA
2059 fillAntenna() ;
2060
2061 // PROCESSOR
2062 fillProcessor() ;
2063
2064 // SOURCE
2065 fillSource() ;
2066
2067 // WEATHER
2068 if ( isWeather_ )
2069 fillWeather() ;
2070
2071 // SYSCAL
2072 fillSysCal() ;
2073
2074 /***
2075 * Start iteration using TableVisitor
2076 ***/
2077 {
2078 static const char *cols[] = {
2079 "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME",
2080 "POLNO",
2081 NULL
2082 };
2083 static const TypeManagerImpl<uInt> tmUInt;
2084 static const TypeManagerImpl<Int> tmInt;
2085 static const TypeManagerImpl<Double> tmDouble;
2086 static const TypeManagerImpl<String> tmString;
2087 static const TypeManager *const tms[] = {
2088 &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL
2089 };
2090 //double t0 = mathutil::gettimeofday_sec() ;
2091 MSWriterVisitor myVisitor(table_->table(),*mstable_);
2092 //double t1 = mathutil::gettimeofday_sec() ;
2093 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
2094 String dataColName = "FLOAT_DATA" ;
2095 if ( useData_ )
2096 dataColName = "DATA" ;
2097 myVisitor.dataColumnName( dataColName ) ;
2098 myVisitor.pointingTableName( ptTabName_ ) ;
2099 myVisitor.setSourceRecord( srcRec_ ) ;
2100 //double t2 = mathutil::gettimeofday_sec() ;
2101 traverseTable(table_->table(), cols, tms, &myVisitor);
2102 //double t3 = mathutil::gettimeofday_sec() ;
2103 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
2104 }
2105 /***
2106 * End iteration using TableVisitor
2107 ***/
2108
2109 // ASDM tables
2110 const TableRecord &stKeys = table_->table().keywordSet() ;
2111 TableRecord &msKeys = mstable_->rwKeywordSet() ;
2112 uInt nfields = stKeys.nfields() ;
2113 for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
2114 String kname = stKeys.name( ifield ) ;
2115 if ( kname.find( "ASDM" ) != String::npos ) {
2116 String asdmpath = stKeys.asString( ifield ) ;
2117 os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
2118 if ( Table::isReadable( asdmpath ) ) {
2119 Table newAsdmTab( asdmpath, Table::Old ) ;
2120 newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
2121 os_ << "add subtable: " << kname << LogIO::POST ;
2122 msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
2123 }
2124 }
2125 }
2126
2127 // replace POINTING table with original one if exists
2128 if ( ptTabName_ != "" ) {
2129 delete mstable_ ;
2130 mstable_ = 0 ;
2131 Table newPtTab( ptTabName_, Table::Old ) ;
2132 newPtTab.copy( filename_+"/POINTING", Table::New ) ;
2133 }
2134
2135 //double endSec = mathutil::gettimeofday_sec() ;
2136 //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2137
2138 os_ << "Exported data as MS" << LogIO::POST ;
2139
2140 return True ;
2141}
2142
2143void MSWriter::init()
2144{
2145// os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
2146// double startSec = mathutil::gettimeofday_sec() ;
2147// os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
2148
2149 // access to scantable
2150 header_ = table_->getHeader() ;
2151
2152 // FLOAT_DATA? or DATA?
2153 if ( header_.npol > 2 ) {
2154 useFloatData_ = False ;
2155 useData_ = True ;
2156 }
2157 else {
2158 useFloatData_ = True ;
2159 useData_ = False ;
2160 }
2161
2162 // polarization type
2163 polType_ = header_.poltype ;
2164 if ( polType_ == "" )
2165 polType_ = "stokes" ;
2166 else if ( polType_.find( "linear" ) != String::npos )
2167 polType_ = "linear" ;
2168 else if ( polType_.find( "circular" ) != String::npos )
2169 polType_ = "circular" ;
2170 else if ( polType_.find( "stokes" ) != String::npos )
2171 polType_ = "stokes" ;
2172 else if ( polType_.find( "linpol" ) != String::npos )
2173 polType_ = "linpol" ;
2174 else
2175 polType_ = "notype" ;
2176
2177 // Check if some subtables are exists
2178 Bool isTcal = False ;
2179 if ( table_->tcal().table().nrow() != 0 ) {
2180 ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
2181 if ( col.isDefined( 0 ) ) {
2182 os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
2183 isTcal = True ;
2184 }
2185 else {
2186 os_ << "No TCAL rows" << LogIO::POST ;
2187 }
2188 }
2189 else {
2190 os_ << "No TCAL rows" << LogIO::POST ;
2191 }
2192 if ( table_->weather().table().nrow() != 0 ) {
2193 ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
2194 if ( col.isDefined( 0 ) ) {
2195 os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
2196 isWeather_ =True ;
2197 }
2198 else {
2199 os_ << "No WEATHER rows" << LogIO::POST ;
2200 }
2201 }
2202 else {
2203 os_ << "No WEATHER rows" << LogIO::POST ;
2204 }
2205
2206 // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
2207 if ( header_.nchan != 1 ) {
2208 if ( isTcal ) {
2209 // examine TCAL subtable
2210 Table tcaltab = table_->tcal().table() ;
2211 ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
2212 for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
2213 if ( tcalCol( irow ).size() != 1 )
2214 tcalSpec_ = True ;
2215 }
2216 }
2217 // examine spectral data
2218 TableIterator iter0( table_->table(), "IFNO" ) ;
2219 while( !iter0.pastEnd() ) {
2220 Table t0( iter0.table() ) ;
2221 ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
2222 uInt len = sharedFloatArrCol( 0 ).size() ;
2223 if ( len != 1 ) {
2224 sharedFloatArrCol.attach( t0, "TSYS" ) ;
2225 if ( sharedFloatArrCol( 0 ).size() != 1 )
2226 tsysSpec_ = True ;
2227 }
2228 iter0.next() ;
2229 }
2230 }
2231
2232 // check if reference for POINTING table exists
2233 const TableRecord &rec = table_->table().keywordSet() ;
2234 if ( rec.isDefined( "POINTING" ) ) {
2235 ptTabName_ = rec.asString( "POINTING" ) ;
2236 if ( !Table::isReadable( ptTabName_ ) ) {
2237 ptTabName_ = "" ;
2238 }
2239 }
2240
2241// double endSec = mathutil::gettimeofday_sec() ;
2242// os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2243}
2244
2245void MSWriter::setupMS()
2246{
2247// os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
2248// double startSec = mathutil::gettimeofday_sec() ;
2249// os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
2250
2251 String dunit = table_->getHeader().fluxunit ;
2252
2253 TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
2254 if ( useFloatData_ )
2255 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
2256 else if ( useData_ )
2257 MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
2258
2259 SetupNewTable newtab( filename_, msDesc, Table::New ) ;
2260
2261 mstable_ = new MeasurementSet( newtab ) ;
2262
2263 TableColumn col ;
2264 if ( useFloatData_ )
2265 col.attach( *mstable_, "FLOAT_DATA" ) ;
2266 else if ( useData_ )
2267 col.attach( *mstable_, "DATA" ) ;
2268 col.rwKeywordSet().define( "UNIT", dunit ) ;
2269
2270 // create subtables
2271 TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
2272 SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
2273 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
2274
2275 TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
2276 SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
2277 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
2278
2279 TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
2280 SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
2281 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
2282
2283 TableDesc feedDesc = MSFeed::requiredTableDesc() ;
2284 SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
2285 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
2286
2287 TableDesc fieldDesc = MSField::requiredTableDesc() ;
2288 SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
2289 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
2290
2291 TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
2292 SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
2293 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
2294
2295 TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
2296 SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
2297 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
2298
2299 TableDesc historyDesc = MSHistory::requiredTableDesc() ;
2300 SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
2301 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
2302
2303 TableDesc observationDesc = MSObservation::requiredTableDesc() ;
2304 SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
2305 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
2306
2307 TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
2308 SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
2309 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
2310
2311 TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
2312 SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
2313 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
2314
2315 TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
2316 SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
2317 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
2318
2319 TableDesc sourceDesc = MSSource::requiredTableDesc() ;
2320 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
2321 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
2322 MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
2323 SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
2324 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
2325
2326 TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
2327 SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
2328 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
2329
2330 TableDesc stateDesc = MSState::requiredTableDesc() ;
2331 SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
2332 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
2333
2334 TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
2335 if ( tcalSpec_ )
2336 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
2337 else
2338 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ;
2339 if ( tsysSpec_ )
2340 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
2341 else
2342 MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ;
2343 SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
2344 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
2345
2346 TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
2347 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
2348 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
2349 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
2350 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
2351 MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
2352 SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
2353 mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
2354
2355 mstable_->initRefs() ;
2356
2357// double endSec = mathutil::gettimeofday_sec() ;
2358// os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2359}
2360
2361void MSWriter::fillObservation()
2362{
2363 //double startSec = mathutil::gettimeofday_sec() ;
2364 //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
2365
2366 // only 1 row
2367 mstable_->observation().addRow( 1, True ) ;
2368 MSObservationColumns msObsCols( mstable_->observation() ) ;
2369 msObsCols.observer().put( 0, header_.observer ) ;
2370 // tentatively put antennaname (from ANTENNA subtable)
2371 String hAntennaName = header_.antennaname ;
2372 String::size_type pos = hAntennaName.find( "//" ) ;
2373 String telescopeName ;
2374 if ( pos != String::npos ) {
2375 telescopeName = hAntennaName.substr( 0, pos ) ;
2376 }
2377 else {
2378 pos = hAntennaName.find( "@" ) ;
2379 telescopeName = hAntennaName.substr( 0, pos ) ;
2380 }
2381// os_ << "telescopeName = " << telescopeName << LogIO::POST ;
2382 msObsCols.telescopeName().put( 0, telescopeName ) ;
2383 msObsCols.project().put( 0, header_.project ) ;
2384 //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
2385 Table sortedtable = table_->table().sort("TIME") ;
2386 ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
2387 Vector<MEpoch> trange( 2 ) ;
2388 trange[0] = timeCol( 0 ) ;
2389 trange[1] = timeCol( table_->nrow()-1 ) ;
2390 msObsCols.timeRangeMeas().put( 0, trange ) ;
2391
2392 //double endSec = mathutil::gettimeofday_sec() ;
2393 //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2394}
2395
2396void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d )
2397{
2398 name.upcase() ;
2399
2400 m = "ALT-AZ" ;
2401 t = "GROUND-BASED" ;
2402 if ( name.matches( Regex( "DV[0-9]+$" ) )
2403 || name.matches( Regex( "DA[0-9]+$" ) )
2404 || name.matches( Regex( "PM[0-9]+$" ) ) )
2405 d = 12.0 ;
2406 else if ( name.matches( Regex( "CM[0-9]+$" ) ) )
2407 d = 7.0 ;
2408 else if ( name.contains( "GBT" ) )
2409 d = 104.9 ;
2410 else if ( name.contains( "MOPRA" ) )
2411 d = 22.0 ;
2412 else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) )
2413 d = 64.0 ;
2414 else if ( name.contains( "TIDBINBILLA" ) )
2415 d = 70.0 ;
2416 else if ( name.contains( "CEDUNA" ) )
2417 d = 30.0 ;
2418 else if ( name.contains( "HOBART" ) )
2419 d = 26.0 ;
2420 else if ( name.contains( "APEX" ) )
2421 d = 12.0 ;
2422 else if ( name.contains( "ASTE" ) )
2423 d = 10.0 ;
2424 else if ( name.contains( "NRO" ) )
2425 d = 45.0 ;
2426 else
2427 d = 1.0 ;
2428}
2429
2430void MSWriter::fillAntenna()
2431{
2432 //double startSec = mathutil::gettimeofday_sec() ;
2433 //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
2434
2435 // only 1 row
2436 Table anttab = mstable_->antenna() ;
2437 anttab.addRow( 1, True ) ;
2438
2439 Table &table = table_->table() ;
2440 const TableRecord &keys = table.keywordSet() ;
2441 String hAntName = keys.asString( "AntennaName" ) ;
2442 String::size_type pos = hAntName.find( "//" ) ;
2443 String antennaName ;
2444 String stationName ;
2445 if ( pos != String::npos ) {
2446 stationName = hAntName.substr( 0, pos ) ;
2447 hAntName = hAntName.substr( pos+2 ) ;
2448 }
2449 pos = hAntName.find( "@" ) ;
2450 if ( pos != String::npos ) {
2451 antennaName = hAntName.substr( 0, pos ) ;
2452 stationName = hAntName.substr( pos+1 ) ;
2453 }
2454 else {
2455 antennaName = hAntName ;
2456 }
2457 Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
2458
2459 String mount, atype ;
2460 Double diameter ;
2461 antennaProperty( antennaName, mount, atype, diameter ) ;
2462
2463 TableRow tr( anttab ) ;
2464 TableRecord &r = tr.record() ;
2465 RecordFieldPtr<String> nameRF( r, "NAME" ) ;
2466 RecordFieldPtr<String> stationRF( r, "STATION" ) ;
2467 RecordFieldPtr<String> mountRF( r, "MOUNT" ) ;
2468 RecordFieldPtr<String> typeRF( r, "TYPE" ) ;
2469 RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ;
2470 RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ;
2471 *nameRF = antennaName ;
2472 *mountRF = mount ;
2473 *typeRF = atype ;
2474 *dishDiameterRF = diameter ;
2475 *positionRF = antpos ;
2476 *stationRF = stationName ;
2477
2478 tr.put( 0 ) ;
2479
2480 //double endSec = mathutil::gettimeofday_sec() ;
2481 //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2482}
2483
2484void MSWriter::fillProcessor()
2485{
2486// double startSec = mathutil::gettimeofday_sec() ;
2487// os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
2488
2489 // only add empty 1 row
2490 MSProcessor msProc = mstable_->processor() ;
2491 msProc.addRow( 1, True ) ;
2492
2493// double endSec = mathutil::gettimeofday_sec() ;
2494// os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2495}
2496
2497void MSWriter::fillSource()
2498{
2499// double startSec = mathutil::gettimeofday_sec() ;
2500// os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
2501
2502 // access to MS SOURCE subtable
2503 MSSource msSrc = mstable_->source() ;
2504
2505 // access to MOLECULE subtable
2506 STMolecules stm = table_->molecules() ;
2507
2508 Int srcId = 0 ;
2509
2510 // row based
2511 TableRow row( msSrc ) ;
2512 TableRecord &rec = row.record() ;
2513 RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
2514 RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
2515 RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
2516 RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
2517 RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
2518 RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
2519 RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
2520 RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
2521 RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
2522 RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
2523 RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
2524
2525 //
2526 // ITERATION: SRCNAME
2527 //
2528 TableIterator iter0( table_->table(), "SRCNAME" ) ;
2529 while( !iter0.pastEnd() ) {
2530 //Table t0( iter0.table() ) ;
2531 Table t0 = iter0.table() ;
2532
2533 // get necessary information
2534 ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
2535 String srcName = srcNameCol( 0 ) ;
2536 ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
2537 Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
2538 sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
2539 Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
2540 ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
2541 Double srcVel = srcVelCol( 0 ) ;
2542 srcRec_.define( srcName, srcId ) ;
2543
2544 // NAME
2545 *nameRF = srcName ;
2546
2547 // SOURCE_ID
2548 *srcidRF = srcId ;
2549
2550 // PROPER_MOTION
2551 *srcpmRF = srcPM ;
2552
2553 // DIRECTION
2554 *srcdirRF = srcDir ;
2555
2556 //
2557 // ITERATION: MOLECULE_ID
2558 //
2559 TableIterator iter1( t0, "MOLECULE_ID" ) ;
2560 while( !iter1.pastEnd() ) {
2561 //Table t1( iter1.table() ) ;
2562 Table t1 = iter1.table() ;
2563
2564 // get necessary information
2565 Vector<Double> restFreq ;
2566 Vector<String> molName ;
2567 Vector<String> fMolName ;
2568 ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
2569 uInt molId = molIdCol( 0 ) ;
2570 stm.getEntry( restFreq, molName, fMolName, molId ) ;
2571
2572 uInt numFreq = restFreq.size() ;
2573
2574 // NUM_LINES
2575 *numlineRF = numFreq ;
2576
2577 // REST_FREQUENCY
2578 restfreqRF.define(restFreq);
2579
2580 // TRANSITION
2581 Vector<String> transition ;
2582 if ( fMolName.size() != 0 ) {
2583 transition = fMolName ;
2584 }
2585 else if ( molName.size() != 0 ) {
2586 transition = molName ;
2587 }
2588 else {
2589 transition.resize( numFreq ) ;
2590 transition = "" ;
2591 }
2592 transitionRF.define(transition);
2593
2594 // SYSVEL
2595 Vector<Double> sysvelArr( numFreq, srcVel ) ;
2596 sysvelRF.define(sysvelArr);
2597
2598 //
2599 // ITERATION: IFNO
2600 //
2601 TableIterator iter2( t1, "IFNO" ) ;
2602 while( !iter2.pastEnd() ) {
2603 //Table t2( iter2.table() ) ;
2604 Table t2 = iter2.table() ;
2605 uInt nrow = msSrc.nrow() ;
2606
2607 // get necessary information
2608 ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
2609 uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
2610 Double midTime ;
2611 Double interval ;
2612 getValidTimeRange( midTime, interval, t2 ) ;
2613
2614 // fill SPECTRAL_WINDOW_ID
2615 *spwidRF = ifno ;
2616
2617 // fill TIME, INTERVAL
2618 *timeRF = midTime ;
2619 *intervalRF = interval ;
2620
2621 // add row
2622 msSrc.addRow( 1, True ) ;
2623 row.put( nrow ) ;
2624
2625 iter2.next() ;
2626 }
2627
2628 iter1.next() ;
2629 }
2630
2631 // increment srcId if SRCNAME changed
2632 srcId++ ;
2633
2634 iter0.next() ;
2635 }
2636
2637// double endSec = mathutil::gettimeofday_sec() ;
2638// os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2639}
2640
2641void MSWriter::fillWeather()
2642{
2643// double startSec = mathutil::gettimeofday_sec() ;
2644// os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
2645
2646 // access to MS WEATHER subtable
2647 MSWeather msw = mstable_->weather() ;
2648
2649 // access to WEATHER subtable
2650 Table stw = table_->weather().table() ;
2651 uInt nrow = stw.nrow() ;
2652
2653 if ( nrow == 0 )
2654 return ;
2655
2656 msw.addRow( nrow, True ) ;
2657 MSWeatherColumns mswCols( msw ) ;
2658
2659 // ANTENNA_ID is always 0
2660 Vector<Int> antIdArr( nrow, 0 ) ;
2661 mswCols.antennaId().putColumn( antIdArr ) ;
2662
2663 // fill weather status
2664 ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
2665 mswCols.temperature().putColumn( sharedFloatCol ) ;
2666 sharedFloatCol.attach( stw, "PRESSURE" ) ;
2667 mswCols.pressure().putColumn( sharedFloatCol ) ;
2668 sharedFloatCol.attach( stw, "HUMIDITY" ) ;
2669 mswCols.relHumidity().putColumn( sharedFloatCol ) ;
2670 sharedFloatCol.attach( stw, "WINDSPEED" ) ;
2671 mswCols.windSpeed().putColumn( sharedFloatCol ) ;
2672 sharedFloatCol.attach( stw, "WINDAZ" ) ;
2673 mswCols.windDirection().putColumn( sharedFloatCol ) ;
2674
2675 // fill TIME and INTERVAL
2676 Double midTime ;
2677 Double interval ;
2678 Vector<Double> intervalArr( nrow, 0.0 ) ;
2679 TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2680 while( !iter.pastEnd() ) {
2681 //Table tab( iter.table() ) ;
2682 Table tab = iter.table() ;
2683
2684 ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2685 uInt wid = widCol( 0 ) ;
2686
2687 getValidTimeRange( midTime, interval, tab ) ;
2688 mswCols.time().put( wid, midTime ) ;
2689 intervalArr[wid] = interval ;
2690
2691 iter.next() ;
2692 }
2693 mswCols.interval().putColumn( intervalArr ) ;
2694
2695// double endSec = mathutil::gettimeofday_sec() ;
2696// os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2697}
2698
2699void MSWriter::fillSysCal()
2700{
2701 Table mssc = mstable_->sysCal() ;
2702
2703 {
2704 static const char *cols[] = {
2705 "BEAMNO", "IFNO", "TIME", "POLNO",
2706 NULL
2707 };
2708 static const TypeManagerImpl<uInt> tmUInt;
2709 static const TypeManagerImpl<Double> tmDouble;
2710 static const TypeManager *const tms[] = {
2711 &tmUInt, &tmUInt, &tmDouble, &tmUInt, NULL
2712 };
2713 //double t0 = mathutil::gettimeofday_sec() ;
2714 MSSysCalVisitor myVisitor(table_->table(),mssc);
2715 //double t1 = mathutil::gettimeofday_sec() ;
2716 //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
2717 traverseTable(table_->table(), cols, tms, &myVisitor);
2718 //double t3 = mathutil::gettimeofday_sec() ;
2719 //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
2720 }
2721
2722}
2723
2724void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
2725{
2726// double startSec = mathutil::gettimeofday_sec() ;
2727// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2728
2729 // sort table
2730 //Table stab = tab.sort( "TIME" ) ;
2731
2732 ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2733 Vector<Double> timeArr = timeCol.getColumn() ;
2734 Double minTime ;
2735 Double maxTime ;
2736 minMax( minTime, maxTime, timeArr ) ;
2737 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2738 // unit for TIME
2739 // Scantable: "d"
2740 // MS: "s"
2741 me = midTime ;
2742 interval = ( maxTime - minTime ) * 86400.0 ;
2743
2744// double endSec = mathutil::gettimeofday_sec() ;
2745// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2746}
2747
2748void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2749{
2750// double startSec = mathutil::gettimeofday_sec() ;
2751// os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2752
2753 // sort table
2754 //Table stab = tab.sort( "TIME" ) ;
2755
2756 Double minTime ;
2757 Double maxTime ;
2758 minMax( minTime, maxTime, atime ) ;
2759 Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
2760 // unit for TIME
2761 // Scantable: "d"
2762 // MS: "s"
2763 me = midTime ;
2764 interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
2765
2766// double endSec = mathutil::gettimeofday_sec() ;
2767// os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2768}
2769
2770}
Note: See TracBrowser for help on using the repository browser.