source: trunk/src/MSWriter.cpp @ 2902

Last change on this file since 2902 was 2869, checked in by Takeshi Nakazato, 10 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
Line 
1//
2// C++ Interface: MSWriter
3//
4// Description:
5//
6// This class is specific writer for MS format
7//
8// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
9//
10// Copyright: See COPYING file that comes with this distribution
11//
12//
13#include <assert.h>
14
15#include <set>
16
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>
22#include <casa/Arrays/Cube.h>
23#include <casa/Containers/RecordField.h>
24
25#include <tables/Tables/ExprNode.h>
26#include <tables/Tables/TableDesc.h>
27#include <tables/Tables/SetupNewTab.h>
28#include <tables/Tables/TableIter.h>
29#include <tables/Tables/RefRows.h>
30#include <tables/Tables/TableRow.h>
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>
36#include <ms/MeasurementSets/MSSourceIndex.h>
37
38#include "MSWriter.h"
39#include "STHeader.h"
40#include "STFrequencies.h"
41#include "STMolecules.h"
42#include "STTcal.h"
43#include "MathUtils.h"
44#include "TableTraverse.h"
45
46using namespace casa ;
47using namespace std ;
48
49namespace asap {
50
51class CorrTypeHandler {
52public:
53  CorrTypeHandler()
54  {}
55  virtual ~CorrTypeHandler() {}
56  virtual Vector<Stokes::StokesTypes> corrType() = 0 ;
57  virtual void reset()
58  {
59    npol = 0 ;
60  }
61  void append( uInt polno )
62  {
63    polnos[npol] = polno ;
64    npol++ ;
65  }
66  uInt nPol() { return npol ; }
67protected:
68  Vector<Stokes::StokesTypes> polmap ;
69  uInt polnos[4] ;
70  uInt npol ;
71};
72
73class LinearHandler : public CorrTypeHandler {
74public:
75  LinearHandler()
76    : CorrTypeHandler()
77  {
78    initMap() ;
79  }
80  virtual ~LinearHandler() {}
81  virtual Vector<Stokes::StokesTypes> corrType()
82  {
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]] ;
87    }
88    else if ( npol == 4 ) {
89      ret[0] = polmap[0] ;
90      ret[1] = polmap[2] ;
91      ret[2] = polmap[3] ;
92      ret[3] = polmap[1] ;
93    }
94    else {
95      throw( AipsError("npol > 4") ) ;
96    }
97    return ret ;
98  }
99protected:
100  void initMap()
101  {
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]] ;
123    }
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] ;
129    }
130    else {
131      throw( AipsError("npol > 4") ) ;
132    }
133    return ret ;
134  }
135private:
136  void initMap()
137  {
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]] ;
159    }
160    else {
161      throw( AipsError("npol > 4") ) ;
162    }
163    return ret ;
164  }
165private:
166  void initMap()
167  {
168    polmap.resize( 4 ) ;
169    polmap[0] = Stokes::I ;
170    polmap[1] = Stokes::Q ;
171    polmap[2] = Stokes::U ;
172    polmap[3] = Stokes::V ;
173  }
174};
175class LinPolHandler : public CorrTypeHandler {
176public:
177  LinPolHandler()
178    : CorrTypeHandler()
179  {
180    initMap() ;
181  }
182  virtual ~LinPolHandler() {}
183  virtual Vector<Stokes::StokesTypes> corrType()
184  {
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 ;
194  }
195private:
196  void initMap()
197  {
198    polmap.resize( 2 ) ;
199    polmap[0] = Stokes::Plinear ;
200    polmap[1] = Stokes::Pangle ;
201  }
202};
203
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() {}
219  void post()
220  {
221    postData() ;
222    postFlag() ;
223    postFlagRow() ;
224    postAuxiliary() ;
225  }
226  virtual void reset()
227  {
228    corr->reset() ;
229    flagRow = False ;
230    npol = 0 ;
231    for ( uInt i = 0 ; i < 4 ; i++ )
232      isFilled[i] = False ;   
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 ) ;
240    isFilled[id] = True;
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    // CAS-5545 FLAG_ROW must always be set False
277    //*flagRowRF = anyEQ( flagRow, True ) ;
278    *flagRowRF = False;
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  }
289  uInt filledIndex()
290  {
291    uInt idx = 0;
292    while( !isFilled[idx] && idx < 4 ) ++idx;
293    return idx;
294  }
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 ;
309  Bool isFilled[4] ;
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};
339
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  {
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   
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    }
366  }
367  virtual void postData()
368  {
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    }
377  }
378private:
379  void attachData()
380  {
381    dataRF.attachToRecord( row.record(), "FLOAT_DATA" ) ;
382  }
383  RecordFieldPtr< Matrix<Float> > dataRF;
384};
385
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  {
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   
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    }
411    else if ( npol == 2 ) {
412      flagRF.define( flag( IPosition( 2, 0, 0 ), IPosition( 2, npol-1, nchan-1 ) ) ) ;
413    }
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    }
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
448class BaseMSWriterVisitor: public TableVisitor {
449  const String *lastFieldName;
450  uInt lastRecordNo;
451  uInt lastBeamNo, lastScanNo, lastIfNo, lastPolNo;
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 
467  virtual void enterFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
468  }
469  virtual void leaveFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
470  }
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*/) { }
485
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 ;}
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;
503    uInt polNo;
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
689    initPolarization() ;
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 
720  virtual void enterFieldName(const uInt recordNo, const String &/*columnValue*/) {
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  }
742  virtual void leaveFieldName(const uInt /*recordNo*/, const String &/*columnValue*/) {
743  }
744  virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue) {
745    //printf("%u: BeamNo: %u\n", recordNo, columnValue);
746   
747    feedId = (Int)columnValue ;
748
749    // put value
750    *feed1RF = feedId ;
751    *feed2RF = feedId ;
752  }
753  virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
754  }
755  virtual void enterScanNo(const uInt /*recordNo*/, uInt columnValue) {
756    //printf("%u: ScanNo: %u\n", recordNo, columnValue);
757
758    // put value
759    // CAS-5841: SCANNO should be consistent with MS SCAN_NUMBER
760    *scanNumberRF = (Int)columnValue ;
761  }
762  virtual void leaveScanNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
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() ;
773    holder->setNchan( nchan ) ;
774
775    addSpectralWindow( spwId, freqId ) ;
776
777    addFeed( feedId, spwId ) ;
778  }
779  virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
780  }
781  virtual void enterSrcType(const uInt /*recordNo*/, Int columnValue) {
782    //printf("%u: SrcType: %d\n", recordNo, columnValue);
783
784    Int stateId = addState( columnValue ) ;
785
786    // put value
787    *stateIdRF = stateId ;
788  }
789  virtual void leaveSrcType(const uInt /*recordNo*/, Int /*columnValue*/) {
790  }
791  virtual void enterCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
792    //printf("%u: CycleNo: %u\n", recordNo, columnValue);
793  }
794  virtual void leaveCycleNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
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  }
824  virtual void leaveTime(const uInt /*recordNo*/, Double /*columnValue*/) {
825    if ( holder->nPol() > 0 ) {
826      Int polId = addPolarization() ;
827      Int ddId = addDataDescription( polId, spwId ) ;
828       
829      // put field
830      *dataDescIdRF = ddId ;
831      holder->post() ;
832     
833      // commit row
834      row.put( rowidx ) ;
835      rowidx++ ;
836
837      // reset holder
838      holder->reset() ;
839    }
840  }
841  virtual void enterPolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
842    //printf("%u: PolNo: %d\n", recordNo, columnValue);
843  }
844  virtual void leavePolNo(const uInt /*recordNo*/, uInt /*columnValue*/) {
845  }
846
847  virtual Bool visitRecord(const uInt recordNo,
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*/,
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 ) ;
864    holder->accumulate( polNo, sp, fl, flr ) ;
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() ;
884
885    // fill empty FIELD rows
886    infillField() ;
887  }
888
889  void dataColumnName( String name )
890  {
891    if ( name == "DATA" )
892      holder = new ComplexDataHolder( row, poltype ) ;
893    else if ( name == "FLOAT_DATA" )
894      holder = new FloatDataHolder( row, poltype ) ;
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() ;
982    potab.addRow() ;
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  }
992  Int addPolarization()
993  {
994    Int idx = -1 ;
995    Vector<Int> corrType = holder->corrTypeInt() ;
996    uInt nEntry = polEntry.size() ;
997    for ( uInt i = 0 ; i < nEntry ; i++ ) {
998      if ( polEntry[i].conform( corrType ) && allEQ( polEntry[i], corrType ) ) {
999        idx = i ;
1000        break ;
1001      }
1002    }
1003   
1004    Int numCorr = holder->nPol() ;
1005    Matrix<Int> corrProduct = corrProductTemplate[numCorr].copy() ;
1006    if ( numCorr == 1 && (corrType[0] == Stokes::YY || corrType[0] == Stokes::LL ) ) {
1007      corrProduct = 1;
1008    }
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 ) ;
1022      polEntry[nEntry] = corrType ;
1023    }
1024
1025    return idx ;
1026  }
1027  Int addDataDescription( Int pid, Int sid )
1028  {
1029    Int idx = -1 ;
1030    uInt nItem = 2 ;
1031    uInt len = ddEntry.nelements() ;
1032    uInt nEntry = len / nItem ;
1033    const Int *dd_p = ddEntry.storage() ;
1034    for ( uInt i = 0 ; i < nEntry ; i++ ) {
1035      Int pol = *dd_p ;
1036      dd_p++ ;
1037      Int spw = *dd_p ;
1038      dd_p++ ;
1039      if ( pid == pol && sid == spw ) {
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
1055      ddEntry.resize( len+nItem ) ;
1056      ddEntry[len] = pid ;
1057      ddEntry[len+1] = sid ;
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() ;
1068    Int mfr = freqframe ;
1069    Int oneChan = 1 ;
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 ) ;
1076    putField( "NUM_CHAN", r, oneChan ) ;
1077
1078    for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
1079      if ( nchan[i] == 0 )
1080        tr.put( i ) ;
1081    }
1082  }
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  }
1100  void addSpectralWindow( Int sid, uInt fid )
1101  {
1102    if (processedFreqId.find((uInt)fid) == processedFreqId.end()
1103        || processedIFNO.find((uInt)sid) == processedIFNO.end() ) {
1104      uInt nrow = spwtab.nrow() ;
1105      while( (Int)nrow <= sid ) {
1106        spwtab.addRow( 1, True ) ;
1107        nrow++ ;
1108      }
1109      processedFreqId.insert((uInt)fid);
1110      processedIFNO.insert((uInt)sid);
1111    }
1112    else {
1113      return ;
1114    }
1115     
1116
1117    Double rp = refpix[fid] ;
1118    Double rv = refval[fid] ;
1119    Double ic = increment[fid] ;
1120
1121    Int mfrInt = (Int)freqframe ;
1122    Int nchan = holder->nChan() ;
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) ) ;
1129    Vector<Double> cw( nchan, ic ) ;
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 ) ;
1141//     defineField( "CHAN_WIDTH", r, res ) ;
1142    defineField( "CHAN_WIDTH", r, cw ) ;
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 ;
1150    uInt nItem = 2 ;
1151    uInt len = feedEntry.nelements() ;
1152    uInt nEntry = len / nItem ;
1153    const Int *fe_p = feedEntry.storage() ;
1154    for ( uInt i = 0 ; i < nEntry ; i++ ) {
1155      Int feed = *fe_p ;
1156      fe_p++ ;
1157      Int spw = *fe_p ;
1158      fe_p++ ;
1159      if ( fid == feed && sid == spw ) {
1160        idx = i ;
1161        break ;
1162      }
1163    }
1164
1165
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
1201      feedEntry.resize( len+nItem ) ;
1202      feedEntry[len] = fid ;
1203      feedEntry[len+1] = sid ;
1204    }
1205  }
1206  void initPolarization()
1207  {
1208    const TableRecord &keys = table.keywordSet() ;
1209    poltype = keys.asString( "POLTYPE" ) ;
1210
1211    initCorrProductTemplate() ;
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  {
1235    //const TableRecord &keys = table.keywordSet() ;
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
1274    //Int id = 0 ;
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  }
1444  void initCorrProductTemplate()
1445  {
1446    Int n = 1 ;
1447    {
1448      Matrix<Int> c( 2, n, 0 ) ;
1449      corrProductTemplate[n] = c ;
1450    }
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    }
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;
1478  CountedPtr<DataHolder> holder;
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;
1515  Block<Int> ddEntry;
1516  Block<Int> feedEntry;
1517  vector< Vector<Int> > polEntry;
1518  set<uInt> processedFreqId;
1519  set<uInt> processedIFNO;
1520  map<uInt,Double> refpix;
1521  map<uInt,Double> refval;
1522  map<uInt,Double> increment;
1523  MFrequency::Types freqframe;
1524  Record srcRec;
1525  map< Int, Matrix<Int> > corrProductTemplate;
1526};
1527
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 
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*/) { }
1550
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 ;}
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  }
1650  virtual ~BaseTsysHolder() {}
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  {}
1702  virtual ~TsysHolder() {}
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  {}
1715  virtual ~TsysSpectrumHolder() {}
1716  virtual Array<Float> getTsys()
1717  {
1718    return tsys ;
1719  }
1720};
1721
1722class BaseTcalProcessor
1723{
1724public:
1725  BaseTcalProcessor( ROArrayColumn<Float> &tcalCol )
1726    : col_( tcalCol )
1727  {}
1728  virtual ~BaseTcalProcessor() {}
1729  void setTcalId( Vector<uInt> &tcalId ) { id_.assign( tcalId ) ; }
1730  virtual Array<Float> getTcal() = 0 ;
1731protected:
1732  ROArrayColumn<Float> col_ ;
1733  Vector<uInt> id_ ;
1734};
1735
1736class TcalProcessor : public BaseTcalProcessor
1737{
1738public:
1739  TcalProcessor( ROArrayColumn<Float> &tcalCol )
1740    : BaseTcalProcessor( tcalCol )
1741  {}
1742  virtual ~TcalProcessor() {}
1743  virtual Array<Float> getTcal()
1744  {
1745    uInt npol = id_.nelements() ;
1746    Vector<Float> tcal( npol ) ;
1747    for ( uInt ipol = 0 ; ipol < npol ; ipol++ )
1748      tcal[ipol] = col_( id_[ipol] ).data()[0] ;
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  {}
1760  virtual ~TcalSpectrumProcessor() {}
1761  virtual Array<Float> getTcal()
1762  {
1763    uInt npol = id_.nelements() ;
1764    //Vector<Float> tcal0 = col_( 0 ) ;
1765    Vector<Float> tcal0 = col_( id_[0] ) ;
1766    uInt nchan = tcal0.nelements() ;
1767    Matrix<Float> tcal( npol, nchan ) ;
1768    tcal.row( 0 ) = tcal0 ;
1769    for ( uInt ipol = 1 ; ipol < npol ; ipol++ )
1770      tcal.row( ipol ) = col_( id_[ipol] ) ;
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 ),
1781      rowidx( 0 ),
1782      polno_()
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
1833  virtual void enterBeamNo(const uInt /*recordNo*/, uInt columnValue)
1834  {
1835    *feedIdRF = (Int)columnValue ;
1836  }
1837  virtual void leaveBeamNo(const uInt /*recordNo*/, uInt /*columnValue*/)
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  }
1850  virtual void leaveIfNo(const uInt /*recordNo*/, uInt /*columnValue*/)
1851  {
1852    //cout << "leaveIfNo" << endl ;
1853    post() ;
1854    reset(true);
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  }
1868  virtual void leaveTime(const uInt /*recordNo*/, Double columnValue)
1869  {
1870    //cout << "leaveTime" << endl ;
1871    if ( isUpdated() ) {
1872      post() ;
1873      *lastTsys = *theTsys ;
1874      lastTcalId = theTcalId ;
1875      reset(false);
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 ) ;
1888    polno_.insert( columnValue ) ;
1889    uInt numPol = polno_.size() ;
1890    if ( lastTsys->nrow() < numPol )
1891      lastTsys->appendTsys( recordNo ) ;
1892    if ( lastTcalId.nelements() <= numPol )
1893      appendTcalId( lastTcalId, tcalId, numPol-1 ) ;
1894    if ( theTsys->nrow() < numPol )
1895      theTsys->appendTsys( recordNo ) ;
1896    else {
1897      theTsys->setTsys( recordNo, numPol-1 ) ;
1898    }
1899    if ( theTcalId.nelements() < numPol )
1900      appendTcalId( theTcalId, tcalId, numPol-1 ) ;
1901    else
1902      setTcalId( theTcalId, tcalId, numPol-1 ) ;
1903  }
1904  virtual void leavePolNo( const uInt /*recordNo*/, uInt /*columnValue*/ )
1905  {
1906  }
1907   
1908private:
1909  void reset(bool completely)
1910  {
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  {
1921    v.resize( polId+1, True ) ;
1922    v[polId] = elem ;
1923  }
1924  void setTcalId( Vector<uInt> &v, uInt &elem, uInt polId )
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() ;
1950    Array<Float> tsys = lastTsys->getTsys() ;
1951    tcalRF.define( tcal ) ;
1952    tsysRF.define( tsys ) ;
1953    sctab.addRow( 1, True ) ;
1954    scrow.put( rowidx ) ;
1955    rowidx++ ;
1956  }
1957 
1958  Bool isUpdated()
1959  {
1960    Bool ret = (anyNE( theTcalId, lastTcalId ) || anyNE( theTsys->getTsys(), lastTsys->getTsys() )) ;
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;
1972  set<uInt> polno_;
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
1984MSWriter::MSWriter(CountedPtr<Scantable> stable)
1985  : table_(stable),
1986    mstable_(NULL),
1987    isWeather_(False),
1988    tcalSpec_(False),
1989    tsysSpec_(False),
1990    ptTabName_("")
1991{
1992  os_ = LogIO() ;
1993  os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
1994//   os_ << "MSWriter::MSWriter()" << LogIO::POST ;
1995
1996  // initialize writer
1997  init() ;
1998}
1999
2000MSWriter::~MSWriter()
2001{
2002  os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
2003//   os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
2004
2005  if ( mstable_ != 0 )
2006    delete mstable_ ;
2007}
2008
2009bool MSWriter::write(const string& filename, const Record& rec)
2010{
2011  os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
2012  //double startSec = mathutil::gettimeofday_sec() ;
2013  //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
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
2053  // PROCESSOR
2054  fillProcessor() ;
2055
2056  // SOURCE
2057  fillSource() ;
2058
2059  // WEATHER
2060  if ( isWeather_ )
2061    fillWeather() ;
2062
2063  // SYSCAL
2064  fillSysCal() ;
2065
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[] = {
2080      &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmUInt, NULL
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 ;
2096  }
2097  /***
2098   * End iteration using TableVisitor
2099   ***/
2100
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
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  }
2126
2127  //double endSec = mathutil::gettimeofday_sec() ;
2128  //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2129
2130  os_ << "Exported data as MS" << LogIO::POST ;
2131
2132  return True ;
2133}
2134
2135void MSWriter::init()
2136{
2137//   os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
2138//   double startSec = mathutil::gettimeofday_sec() ;
2139//   os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
2140 
2141  // access to scantable
2142  header_ = table_->getHeader() ;
2143
2144  // FLOAT_DATA? or DATA?
2145  if ( header_.npol > 2 ) {
2146    useFloatData_ = False ;
2147    useData_ = True ;
2148  }
2149  else {
2150    useFloatData_ = True ;
2151    useData_ = False ;
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
2169  // Check if some subtables are exists
2170  Bool isTcal = False ;
2171  if ( table_->tcal().table().nrow() != 0 ) {
2172    ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
2173    if ( col.isDefined( 0 ) ) {
2174      os_ << "TCAL table exists: nrow=" << table_->tcal().table().nrow() << LogIO::POST ;
2175      isTcal = True ;
2176    }
2177    else {
2178      os_ << "No TCAL rows" << LogIO::POST ;
2179    }
2180  }
2181  else {
2182    os_ << "No TCAL rows" << LogIO::POST ;
2183  }
2184  if ( table_->weather().table().nrow() != 0 ) {
2185    ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
2186    if ( col.isDefined( 0 ) ) {
2187      os_ << "WEATHER table exists: nrow=" << table_->weather().table().nrow() << LogIO::POST ;
2188      isWeather_ =True ;
2189    }
2190    else {
2191      os_ << "No WEATHER rows" << LogIO::POST ;
2192    }
2193  }
2194  else {
2195    os_ << "No WEATHER rows" << LogIO::POST ;
2196  }
2197
2198  // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
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      }
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  }
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
2233//   double endSec = mathutil::gettimeofday_sec() ;
2234//   os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2235}
2236
2237void MSWriter::setupMS()
2238{
2239//   os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
2240//   double startSec = mathutil::gettimeofday_sec() ;
2241//   os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
2242 
2243  String dunit = table_->getHeader().fluxunit ;
2244
2245  TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
2246  if ( useFloatData_ )
2247    MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
2248  else if ( useData_ )
2249    MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
2250
2251  SetupNewTable newtab( filename_, msDesc, Table::New ) ;
2252
2253  mstable_ = new MeasurementSet( newtab ) ;
2254
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
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() ;
2312  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
2313  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
2314  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
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() ;
2327  if ( tcalSpec_ )
2328    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
2329  else
2330    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 1 ) ;
2331  if ( tsysSpec_ )
2332    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
2333  else
2334    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 1 ) ;
2335  SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
2336  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
2337
2338  TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
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 ) ;
2344  SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
2345  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
2346
2347  mstable_->initRefs() ;
2348
2349//   double endSec = mathutil::gettimeofday_sec() ;
2350//   os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2351}
2352
2353void MSWriter::fillObservation()
2354{
2355  //double startSec = mathutil::gettimeofday_sec() ;
2356  //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
2357
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  }
2373//   os_ << "telescopeName = " << telescopeName << LogIO::POST ;
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 ) ;
2383
2384  //double endSec = mathutil::gettimeofday_sec() ;
2385  //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2386}
2387
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
2422void MSWriter::fillAntenna()
2423{
2424  //double startSec = mathutil::gettimeofday_sec() ;
2425  //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
2426
2427  // only 1 row
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( "//" ) ;
2435  String antennaName ;
2436  String stationName ;
2437  if ( pos != String::npos ) {
2438    stationName = hAntName.substr( 0, pos ) ;
2439    hAntName = hAntName.substr( pos+2 ) ;
2440  }
2441  pos = hAntName.find( "@" ) ;
2442  if ( pos != String::npos ) {
2443    antennaName = hAntName.substr( 0, pos ) ;
2444    stationName = hAntName.substr( pos+1 ) ;
2445  }
2446  else {
2447    antennaName = hAntName ;
2448  }
2449  Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
2450 
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" ) ;
2459  RecordFieldPtr<String> mountRF( r, "MOUNT" ) ;
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 ;
2468  *stationRF = stationName ;
2469 
2470  tr.put( 0 ) ;
2471
2472  //double endSec = mathutil::gettimeofday_sec() ;
2473  //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2474}
2475 
2476void MSWriter::fillProcessor()
2477{
2478//   double startSec = mathutil::gettimeofday_sec() ;
2479//   os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
2480 
2481  // only add empty 1 row
2482  MSProcessor msProc = mstable_->processor() ;
2483  msProc.addRow( 1, True ) ;
2484
2485//   double endSec = mathutil::gettimeofday_sec() ;
2486//   os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2487}
2488
2489void MSWriter::fillSource()
2490{
2491//   double startSec = mathutil::gettimeofday_sec() ;
2492//   os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
2493 
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 ;
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
2517  //
2518  // ITERATION: SRCNAME
2519  //
2520  TableIterator iter0( table_->table(), "SRCNAME" ) ;
2521  while( !iter0.pastEnd() ) {
2522    //Table t0( iter0.table() ) ;
2523    Table t0 =  iter0.table()  ;
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 ) ;
2534    srcRec_.define( srcName, srcId ) ;
2535
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
2548    //
2549    // ITERATION: MOLECULE_ID
2550    //
2551    TableIterator iter1( t0, "MOLECULE_ID" ) ;
2552    while( !iter1.pastEnd() ) {
2553      //Table t1( iter1.table() ) ;
2554      Table t1 = iter1.table() ;
2555
2556      // get necessary information
2557      Vector<Double> restFreq ;
2558      Vector<String> molName ;
2559      Vector<String> fMolName ;
2560      ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
2561      uInt molId = molIdCol( 0 ) ;
2562      stm.getEntry( restFreq, molName, fMolName, molId ) ;
2563
2564      uInt numFreq = restFreq.size() ;
2565     
2566      // NUM_LINES
2567      *numlineRF = numFreq ;
2568
2569      // REST_FREQUENCY
2570      restfreqRF.define(restFreq);
2571
2572      // TRANSITION
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      }
2584      transitionRF.define(transition);
2585
2586      // SYSVEL
2587      Vector<Double> sysvelArr( numFreq, srcVel ) ;
2588      sysvelRF.define(sysvelArr);
2589
2590      //
2591      // ITERATION: IFNO
2592      //
2593      TableIterator iter2( t1, "IFNO" ) ;
2594      while( !iter2.pastEnd() ) {
2595        //Table t2( iter2.table() ) ;
2596        Table t2 = iter2.table() ;
2597        uInt nrow = msSrc.nrow() ;
2598
2599        // get necessary information
2600        ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
2601        uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
2602        Double midTime ;
2603        Double interval ;
2604        getValidTimeRange( midTime, interval, t2 ) ;
2605
2606        // fill SPECTRAL_WINDOW_ID
2607        *spwidRF = ifno ;
2608
2609        // fill TIME, INTERVAL
2610        *timeRF = midTime ;
2611        *intervalRF = interval ;
2612
2613        // add row
2614        msSrc.addRow( 1, True ) ;
2615        row.put( nrow ) ;
2616
2617        iter2.next() ;
2618      }
2619
2620      iter1.next() ;
2621    }
2622
2623    // increment srcId if SRCNAME changed
2624    srcId++ ;
2625
2626    iter0.next() ;
2627  }
2628
2629//   double endSec = mathutil::gettimeofday_sec() ;
2630//   os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2631}
2632
2633void MSWriter::fillWeather()
2634{
2635//   double startSec = mathutil::gettimeofday_sec() ;
2636//   os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
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
2668  Double midTime ;
2669  Double interval ;
2670  Vector<Double> intervalArr( nrow, 0.0 ) ;
2671  TableIterator iter( table_->table(), "WEATHER_ID" ) ;
2672  while( !iter.pastEnd() ) {
2673    //Table tab( iter.table() ) ;
2674    Table tab = iter.table() ;
2675
2676    ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
2677    uInt wid = widCol( 0 ) ;
2678
2679    getValidTimeRange( midTime, interval, tab ) ;
2680    mswCols.time().put( wid, midTime ) ;
2681    intervalArr[wid] = interval ;
2682
2683    iter.next() ;
2684  }
2685  mswCols.interval().putColumn( intervalArr ) ;
2686
2687//   double endSec = mathutil::gettimeofday_sec() ;
2688//   os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2689}
2690
2691void MSWriter::fillSysCal()
2692{
2693  Table mssc = mstable_->sysCal() ;
2694
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 ;
2712  }
2713 
2714}
2715
2716void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
2717{
2718//   double startSec = mathutil::gettimeofday_sec() ;
2719//   os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
2720
2721  // sort table
2722  //Table stab = tab.sort( "TIME" ) ;
2723
2724  ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
2725  Vector<Double> timeArr = timeCol.getColumn() ;
2726  Double minTime ;
2727  Double maxTime ;
2728  minMax( minTime, maxTime, timeArr ) ;
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
2736//   double endSec = mathutil::gettimeofday_sec() ;
2737//   os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2738}
2739
2740void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
2741{
2742//   double startSec = mathutil::gettimeofday_sec() ;
2743//   os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
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 ;
2756  interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
2757
2758//   double endSec = mathutil::gettimeofday_sec() ;
2759//   os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
2760}
2761
2762}
Note: See TracBrowser for help on using the repository browser.