source: trunk/src/MSWriter.cpp@ 2760

Last change on this file since 2760 was 2746, 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...

Fixed a known bug that SYSCAL table is not properly handled when
import/export.


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