source: trunk/src/MSWriter.cpp@ 2608

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave[test101]

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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