source: trunk/src/MSWriter.cpp@ 2887

Last change on this file since 2887 was 2869, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5841

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_sdsave, test_importasdm_sd, all sd regressions

Put in Release Notes: Yes

Module(s): sd

Description: Describe your changes here...

Filler/writer are changed so that SCANNO is consistent with
SCAN_NUMBER in MS and scanNumber in ASDM.

It affects several regressions and unit tests so that tests
are necessary to be updated. This change should be verified
by the updated unit tests/regression tests.


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