source: trunk/src/MSWriter.cpp @ 2291

Last change on this file since 2291 was 2291, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: MSFillerUtils and MSWriterUtils added

Test Programs: sd regressions, test_sdsave

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote MSFiller and MSWriter based on TableTraverse?.
TableTraverse? is changed a bit since it doesn't work on plain table.


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