source: trunk/src/MSWriter.cpp@ 2706

Last change on this file since 2706 was 2643, checked in by ShinnosukeKawakami, 12 years ago

hpc34 to trunk 24th Aug. 2012

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