source: branches/hpc34/src/MSWriter.cpp@ 3137

Last change on this file since 3137 was 2597, 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: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Use Block instead of Matrix to store processed pair of
(feedId,spwId) and (polId,spwId) for addFeed and addDataDescription.


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