source: branches/plotter2/src/MSWriter.cpp@ 2894

Last change on this file since 2894 was 2805, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CSV-2594

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on filling SPECTRAL_WINDOW table.


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