source: trunk/src/MSWriter.cpp @ 2793

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

New Development: No

JIRA Issue: Yes CSV-2333

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Fixed an issue that FIELD table sometimes has empty (invalid) rows.


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