source: trunk/src/MSWriter.cpp @ 2022

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

New Development: No

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Copy ASDM subtable that are optionally created by importasdm task.
For filler, their root names are stored as String in table keyword,
while for writer, existing tables are copied and are stored as Table
in table keyword.
File size: 61.9 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
14#include <casa/OS/File.h>
15#include <casa/OS/RegularFile.h>
16#include <casa/OS/Directory.h>
17#include <casa/OS/SymLink.h>
18#include <casa/BasicSL/String.h>
19#include <casa/Containers/RecordField.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
41#include <ctime>
42#include <sys/time.h>
43
44// Boost
45#include <boost/pool/object_pool.hpp>
46
47using namespace casa ;
48using namespace std ;
49
50namespace asap {
51double MSWriter::gettimeofday_sec()
52{
53  struct timeval tv ;
54  gettimeofday( &tv, NULL ) ;
55  return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
56}
57
58MSWriter::MSWriter(CountedPtr<Scantable> stable)
59  : table_(stable),
60    isTcal_(False),
61    isWeather_(False),
62    tcalSpec_(False),
63    tsysSpec_(False),
64    ptTabName_("")
65{
66  os_ = LogIO() ;
67  os_.origin( LogOrigin( "MSWriter", "MSWriter()", WHERE ) ) ;
68  os_ << "MSWriter::MSWriter()" << LogIO::POST ;
69
70  // initialize writer
71  init() ;
72}
73
74MSWriter::~MSWriter()
75{
76  os_.origin( LogOrigin( "MSWriter", "~MSWriter()", WHERE ) ) ;
77  os_ << "MSWriter::~MSWriter()" << LogIO::POST ;
78
79  if ( mstable_ != 0 )
80    delete mstable_ ;
81}
82 
83bool MSWriter::write(const string& filename, const Record& rec)
84{
85  os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
86  double startSec = gettimeofday_sec() ;
87  os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
88
89  filename_ = filename ;
90
91  // parsing MS options
92  Bool overwrite = False ;
93  if ( rec.isDefined( "ms" ) ) {
94    Record msrec = rec.asRecord( "ms" ) ;
95    if ( msrec.isDefined( "overwrite" ) ) {
96      overwrite = msrec.asBool( "overwrite" ) ;
97    }
98  }
99
100  os_ << "Parsing MS options" << endl ;
101  os_ << "   overwrite = " << overwrite << LogIO::POST ;
102
103  File file( filename_ ) ;
104  if ( file.exists() ) {
105    if ( overwrite ) {
106      os_ << filename_ << " exists. Overwrite existing data... " << LogIO::POST ;
107      if ( file.isRegular() ) RegularFile(file).remove() ;
108      else if ( file.isDirectory() ) Directory(file).removeRecursive() ;
109      else SymLink(file).remove() ;
110    }
111    else {
112      os_ << LogIO::SEVERE << "ERROR: " << filename_ << " exists..." << LogIO::POST ;
113      return False ;
114    }
115  }
116
117  // TEST
118  // memory allocation by boost::object_pool
119//   boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ;
120//   ROTableColumn *tcol = 0 ;
121  //
122
123  // set up MS
124  setupMS() ;
125 
126  // subtables
127  // OBSERVATION
128  fillObservation() ;
129
130  // ANTENNA
131  fillAntenna() ;
132
133  // PROCESSOR
134  fillProcessor() ;
135
136  // SOURCE
137  fillSource() ;
138
139  // WEATHER
140  if ( isWeather_ )
141    fillWeather() ;
142
143  // MAIN
144  // Iterate over several ids
145  Vector<uInt> processedFreqId( 0 ) ;
146  Int defaultFieldId = 0 ;
147
148  // row based
149  TableRow row( *mstable_ ) ;
150  TableRecord &trec = row.record() ;
151  NoticeTarget *dataRF = 0 ;
152  if ( useFloatData_ )
153    dataRF = new RecordFieldPtr< Array<Float> >( trec, "FLOAT_DATA" ) ;
154  else if ( useData_ )
155    dataRF = new RecordFieldPtr< Array<Complex> >( trec, "DATA" ) ;
156  RecordFieldPtr< Array<Bool> > flagRF( trec, "FLAG" ) ;
157  RecordFieldPtr<Bool> flagrowRF( trec, "FLAG_ROW" ) ;
158  RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
159  RecordFieldPtr<Double> timecRF( trec, "TIME_CENTROID" ) ;
160  RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
161  RecordFieldPtr<Double> exposureRF( trec, "EXPOSURE" ) ;
162  RecordFieldPtr< Array<Float> > weightRF( trec, "WEIGHT" ) ;
163  RecordFieldPtr< Array<Float> > sigmaRF( trec, "SIGMA" ) ;
164  RecordFieldPtr<Int> ddidRF( trec, "DATA_DESC_ID" ) ;
165  RecordFieldPtr<Int> stateidRF( trec, "STATE_ID" ) ;
166  RecordFieldPtr< Array<Bool> > flagcatRF( trec, "FLAG_CATEGORY" ) ;
167
168  // OBSERVATION_ID is always 0
169  RecordFieldPtr<Int> intRF( trec, "OBSERVATION_ID" ) ;
170  *intRF = 0 ;
171 
172  // ANTENNA1 and ANTENNA2 are always 0
173  intRF.attachToRecord( trec, "ANTENNA1" ) ;
174  *intRF = 0 ;
175  intRF.attachToRecord( trec, "ANTENNA2" ) ;
176  *intRF = 0 ;
177 
178  // ARRAY_ID is tentatively set to 0
179  intRF.attachToRecord( trec, "ARRAY_ID" ) ;
180  *intRF = 0 ;
181
182  // PROCESSOR_ID is tentatively set to 0
183  intRF.attachToRecord( trec, "PROCESSOR_ID" ) ;
184  *intRF = 0 ;
185
186  // UVW is always [0,0,0]
187  RecordFieldPtr< Array<Double> > uvwRF( trec, "UVW" ) ;
188  *uvwRF = Vector<Double>( 3, 0.0 ) ;
189
190  //
191  // ITERATION: FIELDNAME
192  //
193  TableIterator iter0( table_->table(), "FIELDNAME" ) ;
194  while( !iter0.pastEnd() ) {
195    //Table t0( iter0.table() ) ;
196    Table t0 = iter0.table() ;
197    ROTableColumn sharedCol( t0, "FIELDNAME" ) ;
198    String fieldName = sharedCol.asString( 0 ) ;
199//     tcol = tpoolr->construct( t0, "FIELDNAME" ) ;
200//     String fieldName = tcol->asString( 0 ) ;
201//     tpoolr->destroy( tcol ) ;
202    sharedCol.attach( t0, "SRCNAME" ) ;
203    String srcName = sharedCol.asString( 0 ) ;
204//     tcol = tpoolr->construct( t0, "SRCNAME" ) ;
205//     String srcName = tcol->asString( 0 ) ;
206//     tpoolr->destroy( tcol ) ;
207    sharedCol.attach( t0, "TIME" ) ;
208    Double minTime = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day->sec
209//     tcol = tpoolr->construct( t0, "TIME" ) ;
210//     Double minTime = (Double)(tcol->asdouble( 0 )) * 86400.0 ; // day -> sec
211//     tpoolr->destroy( tcol ) ;
212    ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ;
213    //Vector<Double> scanRate = scanRateCol.getColumn()[0] ;
214    Vector<Double> scanRate = scanRateCol( 0 ) ;
215    String::size_type pos = fieldName.find( "__" ) ;
216    Int fieldId = -1 ;
217    if ( pos != String::npos ) {
218      os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
219      fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
220      fieldName = fieldName.substr( 0, pos ) ;
221    }
222    else {
223      os_ << "use default field id" << LogIO::POST ;
224      fieldId = defaultFieldId ;
225      defaultFieldId++ ;
226    }
227    os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
228
229    // FIELD_ID
230    intRF.attachToRecord( trec, "FIELD_ID" ) ;
231    *intRF = fieldId ;
232
233    //
234    // ITERATION: BEAMNO
235    //
236    TableIterator iter1( t0, "BEAMNO" ) ;
237    while( !iter1.pastEnd() ) {
238      Table t1 = iter1.table() ;
239//       tcol = tpoolr->construct( t1, "BEAMNO" ) ;
240//       uInt beamNo = tcol->asuInt( 0 ) ;
241//       tpoolr->destroy( tcol ) ;
242      sharedCol.attach( t1, "BEAMNO" ) ;
243      uInt beamNo = sharedCol.asuInt( 0 ) ;
244      os_ << "beamNo = " << beamNo << LogIO::POST ;
245
246      // FEED1 and FEED2
247      intRF.attachToRecord( trec, "FEED1" ) ;
248      *intRF = beamNo ;
249      intRF.attachToRecord( trec, "FEED2" ) ;
250      *intRF = beamNo ;
251
252      //
253      // ITERATION: SCANNO
254      //
255      TableIterator iter2( t1, "SCANNO" ) ;
256      while( !iter2.pastEnd() ) {
257        Table t2 = iter2.table() ;
258//         tcol = tpoolr->construct( t2, "SCANNO" ) ;
259//         uInt scanNo = tcol->asuInt( 0 ) ;
260//         tpoolr->destroy( tcol ) ;
261        sharedCol.attach( t2, "SCANNO" ) ;
262        uInt scanNo = sharedCol.asuInt( 0 ) ;
263        os_ << "scanNo = " << scanNo << LogIO::POST ;
264
265        // SCAN_NUMBER
266        // MS: 1-based
267        // Scantable: 0-based
268        intRF.attachToRecord( trec, "SCAN_NUMBER" ) ;
269        *intRF = scanNo + 1 ;
270
271        //
272        // ITERATION: IFNO
273        //
274        TableIterator iter3( t2, "IFNO" ) ;
275        while( !iter3.pastEnd() ) {
276          Table t3 = iter3.table() ;
277//           tcol = tpoolr->construct( t3, "IFNO" ) ;
278//           uInt ifNo = tcol->asuInt( 0 ) ;
279//           tpoolr->destroy( tcol ) ;
280          sharedCol.attach( t3, "IFNO" ) ;
281          uInt ifNo = sharedCol.asuInt( 0 ) ;
282          os_ << "ifNo = " << ifNo << LogIO::POST ;
283//           tcol = tpoolr->construct( t3, "FREQ_ID" ) ;
284//           uInt freqId = tcol->asuInt( 0 ) ;
285//           tpoolr->destroy( tcol ) ;
286          sharedCol.attach( t3, "FREQ_ID" ) ;
287          uInt freqId = sharedCol.asuInt( 0 ) ;
288          os_ << "freqId = " << freqId << LogIO::POST ;
289          Int subscan = 1 ; // 1-base
290          //
291          // ITERATION: SRCTYPE
292          //
293          TableIterator iter4( t3, "SRCTYPE" ) ;
294          while( !iter4.pastEnd() ) {
295            Table t4 = iter4.table() ;
296//             tcol = tpoolr->construct( t4, "SRCTYPE" ) ;
297//             uInt srcType = tcol->asInt( 0 ) ;
298//             tpoolr->destroy( tcol ) ;
299            sharedCol.attach( t4, "SRCTYPE" ) ;
300            Int srcType = sharedCol.asInt( 0 ) ;
301            Int stateId = addState( srcType, subscan ) ;
302            *stateidRF = stateId ;
303            //
304            // ITERATION: CYCLENO and TIME
305            //
306            Block<String> cols( 2 ) ;
307            cols[0] = "CYCLENO" ;
308            cols[1] = "TIME" ;
309            TableIterator iter5( t4, cols ) ;
310            while( !iter5.pastEnd() ) {
311              Table t5 =  iter5.table().sort("POLNO") ;
312              //sharedCol.attach( t5, "CYCLENO" ) ;
313              //uInt cycleNo = sharedCol.asuInt( 0 ) ;
314              Int nrow = t5.nrow() ;
315              os_ << "nrow = " << nrow << LogIO::POST ;
316             
317              Vector<Int> polnos( nrow ) ;
318              indgen( polnos, 0 ) ;
319              Int polid = addPolarization( polnos ) ;
320              os_ << "polid = " << polid << LogIO::POST ;
321             
322              // DATA/FLOAT_DATA
323              ROArrayColumn<Float> specCol( t5, "SPECTRA" ) ;
324              ROArrayColumn<uChar> flagCol( t5, "FLAGTRA" ) ;
325              uInt nchan = specCol( 0 ).size() ;
326              IPosition cellshape( 2, nrow, nchan ) ;
327              if ( useFloatData_ ) {
328                // FLOAT_DATA
329                Matrix<Float> dataArr( cellshape ) ;
330                Matrix<Bool> flagArr( cellshape ) ;
331                Vector<Bool> tmpB ;
332                for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
333                  dataArr.row( ipol ) = specCol( ipol ) ;
334                  tmpB.reference( flagArr.row( ipol ) ) ;
335                  convertArray( tmpB, flagCol( ipol ) ) ;
336                }
337                //*(*((RecordFieldPtr< Array<Float> > *)dataRF)) = dataArr ;
338                ((RecordFieldPtr< Array<Float> > *)dataRF)->define( dataArr ) ;
339               
340                // FLAG
341                //*flagRF = flagArr ;
342                flagRF.define( flagArr ) ;
343              }
344              else if ( useData_ ) {
345                // DATA
346                // assume nrow = 4
347                Matrix<Complex> dataArr( cellshape ) ;
348                Vector<Float> zeroIm( nchan, 0 ) ;
349                Matrix<Float> dummy( IPosition( 2, 2, nchan ) ) ;
350                dummy.row( 0 ) = specCol( 0 ) ;
351                dummy.row( 1 ) = zeroIm ;
352                dataArr.row( 0 ) = RealToComplex( dummy ) ;
353                dummy.row( 0 ) = specCol( 1 ) ;
354                dataArr.row( 3 ) = RealToComplex( dummy ) ;
355                dummy.row( 0 ) = specCol( 2 ) ;
356                dummy.row( 1 ) = specCol( 3 ) ;
357                dataArr.row( 1 ) = RealToComplex( dummy ) ;
358                dataArr.row( 2 ) = conj( dataArr.row( 1 ) ) ;
359                //*(*((RecordFieldPtr< Array<Complex> > *)dataRF)) = dataArr ;
360                ((RecordFieldPtr< Array<Complex> > *)dataRF)->define( dataArr ) ;
361               
362               
363                // FLAG
364                Matrix<Bool> flagArr( cellshape ) ;
365                Vector<Bool> tmpB ;
366                tmpB.reference( flagArr.row( 0 ) ) ;
367                convertArray( tmpB, flagCol( 0 ) ) ;
368                tmpB.reference( flagArr.row( 3 ) ) ;
369                convertArray( tmpB, flagCol( 3 ) ) ;
370                tmpB.reference( flagArr.row( 1 ) ) ;
371                convertArray( tmpB, ( flagCol( 2 ) | flagCol( 3 ) ) ) ;
372                flagArr.row( 2 ) = flagArr.row( 1 ) ;
373                //*flagRF = flagArr ;
374                flagRF.define( flagArr ) ;
375              }
376
377              // FLAG_ROW
378//               tcol = tpoolr->construct( t5, "FLAGROW" ) ;
379              sharedCol.attach( t5, "FLAGROW" ) ;
380              Vector<uInt> flagRowArr( nrow ) ;
381              for ( Int irow = 0 ; irow < nrow ; irow++ )
382                flagRowArr[irow] = sharedCol.asuInt( irow ) ;
383//                 flagRowArr[irow] = tcol->asuInt( irow ) ;
384//               tpoolr->destroy( tcol ) ;
385              *flagrowRF = anyNE( flagRowArr, (uInt)0 ) ;
386
387              // TIME and TIME_CENTROID
388//               tcol = tpoolr->construct( t5, "TIME" ) ;
389//               Double mTimeV = (Double)(tcol->asdouble(0)) * 86400.0 ; // day->sec
390//               tpoolr->destroy( tcol ) ;
391              sharedCol.attach( t5, "TIME" ) ;
392              Double mTimeV = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day -> sec
393              *timeRF = mTimeV ;
394              *timecRF = mTimeV ;
395
396              // INTERVAL and EXPOSURE
397//               tcol = tpoolr->construct( t5, "INTERVAL" ) ;
398//               Double interval = (Double)(tcol->asdouble( 0 )) ;
399//               tpoolr->destroy( tcol ) ;
400              sharedCol.attach( t5, "INTERVAL" ) ;
401              Double interval = (Double)sharedCol.asdouble( 0 ) ;
402              *intervalRF = interval ;
403              *exposureRF = interval ;
404             
405              // WEIGHT and SIGMA
406              // always 1 at the moment
407              Vector<Float> wArr( nrow, 1.0 ) ;
408//               *weightRF = wArr ;
409//               *sigmaRF = wArr ;
410              weightRF.define( wArr ) ;
411              sigmaRF.define( wArr ) ;
412             
413              // add DATA_DESCRIPTION row
414              Int ddid = addDataDescription( polid, ifNo ) ;
415              os_ << "ddid = " << ddid << LogIO::POST ;
416              *ddidRF = ddid ;
417             
418              // for SYSCAL table
419//               tcol = tpoolr->construct( t5, "TCAL_ID" ) ;
420              sharedCol.attach( t5, "TCAL_ID" ) ;
421              Vector<uInt> tcalIdArr( nrow ) ;
422              for ( Int irow = 0 ; irow < nrow ; irow++ )
423                tcalIdArr[irow] = sharedCol.asuInt( irow ) ;
424//                 tcalIdArr[irow] = tcol->asuInt( irow ) ;
425//               tpoolr->destroy( tcol ) ;
426              os_ << "tcalIdArr = " << tcalIdArr << LogIO::POST ;
427              String key = String::toString( tcalIdArr[0] ) ;
428              if ( !tcalIdRec_.isDefined( key ) ) {
429                tcalIdRec_.define( key, tcalIdArr ) ;
430                tcalRowRec_.define( key, t5.rowNumbers() ) ;
431              }
432              else {
433                Vector<uInt> pastrows = tcalRowRec_.asArrayuInt( key ) ;
434                tcalRowRec_.define( key, concatenateArray( pastrows, t5.rowNumbers() ) ) ;
435              }
436             
437              // fill STATE_ID
438              //ROScalarColumn<Int> srcTypeCol( t5, "SRCTYPE" ) ;
439              //Int srcType = srcTypeCol( 0 ) ;
440              //Int stateId = addState( srcType, subscan ) ;
441              //*stateidRF = stateId ;
442             
443              // for POINTING table
444              if ( ptTabName_ == "" ) {
445                ROArrayColumn<Double> dirCol( t5, "DIRECTION" ) ;
446                Vector<Double> dir = dirCol( 0 ) ;
447                dirCol.attach( t5, "SCANRATE" ) ;
448                Vector<Double> rate = dirCol( 0 ) ;
449                Matrix<Double> msDir( 2, 1 ) ;
450                msDir.column( 0 ) = dir ;
451                if ( anyNE( rate, 0.0 ) ) {
452                  msDir.resize( 2, 2 ) ;
453                  msDir.column( 1 ) = rate ;
454                }
455                addPointing( fieldName, mTimeV, interval, msDir ) ;
456              }
457             
458              // FLAG_CATEGORY is tentatively set
459              //*flagcatRF = Cube<Bool>( nrow, nchan, 1, False ) ;
460              flagcatRF.define( Cube<Bool>( nrow, nchan, 1, False ) ) ;
461             
462              // add row
463              mstable_->addRow( 1, True ) ;
464              row.put( mstable_->nrow()-1 ) ;
465             
466              iter5.next() ;
467            }
468           
469            iter4.next() ;
470          }
471
472          // add SPECTRAL_WINDOW row
473          if ( allNE( processedFreqId, freqId ) ) {
474            uInt vsize = processedFreqId.size() ;
475            processedFreqId.resize( vsize+1, True ) ;
476            processedFreqId[vsize] = freqId ;
477            addSpectralWindow( ifNo, freqId ) ;
478          }
479         
480          iter3.next() ;
481        }
482       
483        iter2.next() ;
484      }
485     
486      // add FEED row
487      addFeed( beamNo ) ;
488     
489      iter1.next() ;
490    }
491   
492    // add FIELD row
493    addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
494
495    iter0.next() ;
496  }
497
498//   delete tpoolr ;
499  delete dataRF ;
500
501  // SYSCAL
502  if ( isTcal_ )
503    fillSysCal() ;
504
505  // ASDM tables
506  const TableRecord &stKeys = table_->table().keywordSet() ;
507  TableRecord &msKeys = mstable_->rwKeywordSet() ;
508  uInt nfields = stKeys.nfields() ;
509  for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
510    String kname = stKeys.name( ifield ) ;
511    if ( kname.find( "ASDM" ) != String::npos ) {
512      String asdmpath = stKeys.asString( ifield ) ;
513      os_ << "found ASDM table: " << asdmpath << LogIO::POST ;
514      if ( Table::isReadable( asdmpath ) ) {
515        Table newAsdmTab( asdmpath, Table::Old ) ;
516        newAsdmTab.copy( filename_+"/"+kname, Table::New ) ;
517        os_ << "add subtable: " << kname << LogIO::POST ;
518        msKeys.defineTable( kname, Table( filename_+"/"+kname, Table::Old ) ) ;
519      }
520    }
521  }
522
523  // replace POINTING table with original one if exists
524  if ( ptTabName_ != "" ) {
525    delete mstable_ ;
526    mstable_ = 0 ;
527    Table newPtTab( ptTabName_, Table::Old ) ;
528    newPtTab.copy( filename_+"/POINTING", Table::New ) ;
529  }
530
531  double endSec = gettimeofday_sec() ;
532  os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
533
534  return True ;
535}
536 
537void MSWriter::init()
538{
539//   os_.origin( LogOrigin( "MSWriter", "init()", WHERE ) ) ;
540  double startSec = gettimeofday_sec() ;
541  os_ << "start MSWriter::init() startSec=" << startSec << LogIO::POST ;
542 
543  // access to scantable
544  header_ = table_->getHeader() ;
545
546  // FLOAT_DATA? or DATA?
547  if ( header_.npol > 2 ) {
548    useFloatData_ = False ;
549    useData_ = True ;
550  }
551  else {
552    useFloatData_ = True ;
553    useData_ = False ;
554  }
555
556  // polarization type
557  polType_ = header_.poltype ;
558  if ( polType_ == "" )
559    polType_ = "stokes" ;
560  else if ( polType_.find( "linear" ) != String::npos )
561    polType_ = "linear" ;
562  else if ( polType_.find( "circular" ) != String::npos )
563    polType_ = "circular" ;
564  else if ( polType_.find( "stokes" ) != String::npos )
565    polType_ = "stokes" ;
566  else if ( polType_.find( "linpol" ) != String::npos )
567    polType_ = "linpol" ;
568  else
569    polType_ = "notype" ;
570
571  // Check if some subtables are exists
572  if ( table_->tcal().table().nrow() != 0 ) {
573    ROTableColumn col( table_->tcal().table(), "TCAL" ) ;
574    if ( col.isDefined( 0 ) ) {
575      os_ << "nrow=" << table_->tcal().table().nrow() << ": TCAL table exists" << LogIO::POST ;
576      isTcal_ = True ;
577    }
578    else {
579      os_ << "no TCAL rows" << LogIO::POST ;
580    }
581  }
582  else {
583    os_ << "no TCAL rows" << LogIO::POST ;
584  }
585  if ( table_->weather().table().nrow() != 0 ) {
586    ROTableColumn col( table_->weather().table(), "TEMPERATURE" ) ;
587    if ( col.isDefined( 0 ) ) {
588      os_ << "nrow=" << table_->weather().table().nrow() << ": WEATHER table exists" << LogIO::POST ;
589      isWeather_ =True ;
590    }
591    else {
592      os_ << "no WEATHER rows" << LogIO::POST ;
593    }
594  }
595  else {
596    os_ << "no WEATHER rows" << LogIO::POST ;
597  }
598
599  // Are TCAL_SPECTRUM and TSYS_SPECTRUM necessary?
600  if ( isTcal_ && header_.nchan != 1 ) {
601    // examine TCAL subtable
602    Table tcaltab = table_->tcal().table() ;
603    ROArrayColumn<Float> tcalCol( tcaltab, "TCAL" ) ;
604    for ( uInt irow = 0 ; irow < tcaltab.nrow() ; irow++ ) {
605      if ( tcalCol( irow ).size() != 1 )
606        tcalSpec_ = True ;
607    }
608    // examine spectral data
609    TableIterator iter0( table_->table(), "IFNO" ) ;
610    while( !iter0.pastEnd() ) {
611      Table t0( iter0.table() ) ;
612      ROArrayColumn<Float> sharedFloatArrCol( t0, "SPECTRA" ) ;
613      uInt len = sharedFloatArrCol( 0 ).size() ;
614      if ( len != 1 ) {
615        sharedFloatArrCol.attach( t0, "TSYS" ) ;
616        if ( sharedFloatArrCol( 0 ).size() != 1 )
617          tsysSpec_ = True ;
618      }
619      iter0.next() ;
620    }
621  }
622
623  // check if reference for POINTING table exists
624  const TableRecord &rec = table_->table().keywordSet() ;
625  if ( rec.isDefined( "POINTING" ) ) {
626    ptTabName_ = rec.asString( "POINTING" ) ;
627    if ( !Table::isReadable( ptTabName_ ) ) {
628      ptTabName_ = "" ;
629    }
630  }
631
632  double endSec = gettimeofday_sec() ;
633  os_ << "end MSWriter::init() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
634}
635
636void MSWriter::setupMS()
637{
638//   os_.origin( LogOrigin( "MSWriter", "setupMS()", WHERE ) ) ;
639  double startSec = gettimeofday_sec() ;
640  os_ << "start MSWriter::setupMS() startSec=" << startSec << LogIO::POST ;
641
642  TableDesc msDesc = MeasurementSet::requiredTableDesc() ;
643  if ( useFloatData_ )
644    MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::FLOAT_DATA, 2 ) ;
645  else if ( useData_ )
646    MeasurementSet::addColumnToDesc( msDesc, MSMainEnums::DATA, 2 ) ;
647
648  SetupNewTable newtab( filename_, msDesc, Table::New ) ;
649
650  mstable_ = new MeasurementSet( newtab ) ;
651
652  // create subtables
653  TableDesc antennaDesc = MSAntenna::requiredTableDesc() ;
654  SetupNewTable antennaTab( mstable_->antennaTableName(), antennaDesc, Table::New ) ;
655  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::ANTENNA ), Table( antennaTab ) ) ;
656
657  TableDesc dataDescDesc = MSDataDescription::requiredTableDesc() ;
658  SetupNewTable dataDescTab( mstable_->dataDescriptionTableName(), dataDescDesc, Table::New ) ;
659  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DATA_DESCRIPTION ), Table( dataDescTab ) ) ;
660
661  TableDesc dopplerDesc = MSDoppler::requiredTableDesc() ;
662  SetupNewTable dopplerTab( mstable_->dopplerTableName(), dopplerDesc, Table::New ) ;
663  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::DOPPLER ), Table( dopplerTab ) ) ;
664
665  TableDesc feedDesc = MSFeed::requiredTableDesc() ;
666  SetupNewTable feedTab( mstable_->feedTableName(), feedDesc, Table::New ) ;
667  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FEED ), Table( feedTab ) ) ;
668
669  TableDesc fieldDesc = MSField::requiredTableDesc() ;
670  SetupNewTable fieldTab( mstable_->fieldTableName(), fieldDesc, Table::New ) ;
671  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FIELD ), Table( fieldTab ) ) ;
672
673  TableDesc flagCmdDesc = MSFlagCmd::requiredTableDesc() ;
674  SetupNewTable flagCmdTab( mstable_->flagCmdTableName(), flagCmdDesc, Table::New ) ;
675  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FLAG_CMD ), Table( flagCmdTab ) ) ;
676
677  TableDesc freqOffsetDesc = MSFreqOffset::requiredTableDesc() ;
678  SetupNewTable freqOffsetTab( mstable_->freqOffsetTableName(), freqOffsetDesc, Table::New ) ;
679  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::FREQ_OFFSET ), Table( freqOffsetTab ) ) ;
680
681  TableDesc historyDesc = MSHistory::requiredTableDesc() ;
682  SetupNewTable historyTab( mstable_->historyTableName(), historyDesc, Table::New ) ;
683  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::HISTORY ), Table( historyTab ) ) ;
684
685  TableDesc observationDesc = MSObservation::requiredTableDesc() ;
686  SetupNewTable observationTab( mstable_->observationTableName(), observationDesc, Table::New ) ;
687  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::OBSERVATION ), Table( observationTab ) ) ;
688
689  TableDesc pointingDesc = MSPointing::requiredTableDesc() ;
690  SetupNewTable pointingTab( mstable_->pointingTableName(), pointingDesc, Table::New ) ;
691  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POINTING ), Table( pointingTab ) ) ;
692
693  TableDesc polarizationDesc = MSPolarization::requiredTableDesc() ;
694  SetupNewTable polarizationTab( mstable_->polarizationTableName(), polarizationDesc, Table::New ) ;
695  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::POLARIZATION ), Table( polarizationTab ) ) ;
696
697  TableDesc processorDesc = MSProcessor::requiredTableDesc() ;
698  SetupNewTable processorTab( mstable_->processorTableName(), processorDesc, Table::New ) ;
699  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::PROCESSOR ), Table( processorTab ) ) ;
700
701  TableDesc sourceDesc = MSSource::requiredTableDesc() ;
702  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::TRANSITION, 1 ) ;
703  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::REST_FREQUENCY, 1 ) ;
704  MSSource::addColumnToDesc( sourceDesc, MSSourceEnums::SYSVEL, 1 ) ;
705  SetupNewTable sourceTab( mstable_->sourceTableName(), sourceDesc, Table::New ) ;
706  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SOURCE ), Table( sourceTab ) ) ;
707
708  TableDesc spwDesc = MSSpectralWindow::requiredTableDesc() ;
709  SetupNewTable spwTab( mstable_->spectralWindowTableName(), spwDesc, Table::New ) ;
710  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SPECTRAL_WINDOW ), Table( spwTab ) ) ;
711
712  TableDesc stateDesc = MSState::requiredTableDesc() ;
713  SetupNewTable stateTab( mstable_->stateTableName(), stateDesc, Table::New ) ;
714  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::STATE ), Table( stateTab ) ) ;
715
716  // TODO: add TCAL_SPECTRUM and TSYS_SPECTRUM if necessary
717  TableDesc sysCalDesc = MSSysCal::requiredTableDesc() ;
718  MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL, 2 ) ;
719  MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS, 2 ) ;
720  if ( tcalSpec_ )
721    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TCAL_SPECTRUM, 2 ) ;
722  if ( tsysSpec_ )
723    MSSysCal::addColumnToDesc( sysCalDesc, MSSysCalEnums::TSYS_SPECTRUM, 2 ) ;
724  SetupNewTable sysCalTab( mstable_->sysCalTableName(), sysCalDesc, Table::New ) ;
725  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::SYSCAL ), Table( sysCalTab ) ) ;
726
727  TableDesc weatherDesc = MSWeather::requiredTableDesc() ;
728  MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::TEMPERATURE ) ;
729  MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::PRESSURE ) ;
730  MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::REL_HUMIDITY ) ;
731  MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_SPEED ) ;
732  MSWeather::addColumnToDesc( weatherDesc, MSWeatherEnums::WIND_DIRECTION ) ;
733  SetupNewTable weatherTab( mstable_->weatherTableName(), weatherDesc, Table::New ) ;
734  mstable_->rwKeywordSet().defineTable( MeasurementSet::keywordName( MeasurementSet::WEATHER ), Table( weatherTab ) ) ;
735
736  mstable_->initRefs() ;
737
738  double endSec = gettimeofday_sec() ;
739  os_ << "end MSWriter::setupMS() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
740}
741
742void MSWriter::fillObservation()
743{
744  double startSec = gettimeofday_sec() ;
745  os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
746
747  // only 1 row
748  mstable_->observation().addRow( 1, True ) ;
749  MSObservationColumns msObsCols( mstable_->observation() ) ;
750  msObsCols.observer().put( 0, header_.observer ) ;
751  // tentatively put antennaname (from ANTENNA subtable)
752  String hAntennaName = header_.antennaname ;
753  String::size_type pos = hAntennaName.find( "//" ) ;
754  String telescopeName ;
755  if ( pos != String::npos ) {
756    telescopeName = hAntennaName.substr( 0, pos ) ;
757  }
758  else {
759    pos = hAntennaName.find( "@" ) ;
760    telescopeName = hAntennaName.substr( 0, pos ) ;
761  }
762  os_ << "telescopeName = " << telescopeName << LogIO::POST ;
763  msObsCols.telescopeName().put( 0, telescopeName ) ;
764  msObsCols.project().put( 0, header_.project ) ;
765  //ScalarMeasColumn<MEpoch> timeCol( table_->table().sort("TIME"), "TIME" ) ;
766  Table sortedtable = table_->table().sort("TIME") ;
767  ScalarMeasColumn<MEpoch> timeCol( sortedtable, "TIME" ) ;
768  Vector<MEpoch> trange( 2 ) ;
769  trange[0] = timeCol( 0 ) ;
770  trange[1] = timeCol( table_->nrow()-1 ) ;
771  msObsCols.timeRangeMeas().put( 0, trange ) ;
772
773  double endSec = gettimeofday_sec() ;
774  os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
775}
776
777void MSWriter::fillAntenna()
778{
779  double startSec = gettimeofday_sec() ;
780  os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
781
782  // only 1 row
783  mstable_->antenna().addRow( 1, True ) ;
784  MSAntennaColumns msAntCols( mstable_->antenna() ) ;
785
786  String hAntennaName = header_.antennaname ;
787  String::size_type pos = hAntennaName.find( "//" ) ;
788  String antennaName ;
789  String stationName ;
790  if ( pos != String::npos ) {
791    hAntennaName = hAntennaName.substr( pos+2 ) ;
792  }
793  pos = hAntennaName.find( "@" ) ;
794  if ( pos != String::npos ) {
795    antennaName = hAntennaName.substr( 0, pos ) ;
796    stationName = hAntennaName.substr( pos+1 ) ;
797  }
798  else {
799    antennaName = hAntennaName ;
800    stationName = hAntennaName ;
801  }
802  os_ << "antennaName = " << antennaName << LogIO::POST ;
803  os_ << "stationName = " << stationName << LogIO::POST ;
804 
805  msAntCols.name().put( 0, antennaName ) ;
806  msAntCols.station().put( 0, stationName ) ;
807
808  os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
809 
810  msAntCols.position().put( 0, header_.antennaposition ) ;
811
812  double endSec = gettimeofday_sec() ;
813  os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
814}
815
816void MSWriter::fillProcessor()
817{
818  double startSec = gettimeofday_sec() ;
819  os_ << "start MSWriter::fillProcessor() startSec=" << startSec << LogIO::POST ;
820 
821  // only add empty 1 row
822  MSProcessor msProc = mstable_->processor() ;
823  msProc.addRow( 1, True ) ;
824
825  double endSec = gettimeofday_sec() ;
826  os_ << "end MSWriter::fillProcessor() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
827}
828
829void MSWriter::fillSource()
830{
831  double startSec = gettimeofday_sec() ;
832  os_ << "start MSWriter::fillSource() startSec=" << startSec << LogIO::POST ;
833 
834  // access to MS SOURCE subtable
835  MSSource msSrc = mstable_->source() ;
836
837  // access to MOLECULE subtable
838  STMolecules stm = table_->molecules() ;
839
840  Int srcId = 0 ;
841  Vector<Double> restFreq ;
842  Vector<String> molName ;
843  Vector<String> fMolName ;
844
845  // row based
846  TableRow row( msSrc ) ;
847  TableRecord &rec = row.record() ;
848  RecordFieldPtr<Int> srcidRF( rec, "SOURCE_ID" ) ;
849  RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
850  RecordFieldPtr< Array<Double> > srcpmRF( rec, "PROPER_MOTION" ) ;
851  RecordFieldPtr< Array<Double> > srcdirRF( rec, "DIRECTION" ) ;
852  RecordFieldPtr<Int> numlineRF( rec, "NUM_LINES" ) ;
853  RecordFieldPtr< Array<Double> > restfreqRF( rec, "REST_FREQUENCY" ) ;
854  RecordFieldPtr< Array<Double> > sysvelRF( rec, "SYSVEL" ) ;
855  RecordFieldPtr< Array<String> > transitionRF( rec, "TRANSITION" ) ;
856  RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
857  RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
858  RecordFieldPtr<Int> spwidRF( rec, "SPECTRAL_WINDOW_ID" ) ;
859
860  //
861  // ITERATION: SRCNAME
862  //
863  TableIterator iter0( table_->table(), "SRCNAME" ) ;
864  while( !iter0.pastEnd() ) {
865    //Table t0( iter0.table() ) ;
866    Table t0 =  iter0.table()  ;
867
868    // get necessary information
869    ROScalarColumn<String> srcNameCol( t0, "SRCNAME" ) ;
870    String srcName = srcNameCol( 0 ) ;
871    ROArrayColumn<Double> sharedDArrRCol( t0, "SRCPROPERMOTION" ) ;
872    Vector<Double> srcPM = sharedDArrRCol( 0 ) ;
873    sharedDArrRCol.attach( t0, "SRCDIRECTION" ) ;
874    Vector<Double> srcDir = sharedDArrRCol( 0 ) ;
875    ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
876    Double srcVel = srcVelCol( 0 ) ;
877
878    // NAME
879    *nameRF = srcName ;
880
881    // SOURCE_ID
882    *srcidRF = srcId ;
883
884    // PROPER_MOTION
885    *srcpmRF = srcPM ;
886   
887    // DIRECTION
888    *srcdirRF = srcDir ;
889
890    //
891    // ITERATION: MOLECULE_ID
892    //
893    TableIterator iter1( t0, "MOLECULE_ID" ) ;
894    while( !iter1.pastEnd() ) {
895      //Table t1( iter1.table() ) ;
896      Table t1 = iter1.table() ;
897
898      // get necessary information
899      ROScalarColumn<uInt> molIdCol( t1, "MOLECULE_ID" ) ;
900      uInt molId = molIdCol( 0 ) ;
901      stm.getEntry( restFreq, molName, fMolName, molId ) ;
902
903      uInt numFreq = restFreq.size() ;
904     
905      // NUM_LINES
906      *numlineRF = numFreq ;
907
908      // REST_FREQUENCY
909      *restfreqRF = restFreq ;
910
911      // TRANSITION
912      //*transitionRF = fMolName ;
913      Vector<String> transition ;
914      if ( fMolName.size() != 0 ) {
915        transition = fMolName ;
916      }
917      else if ( molName.size() != 0 ) {
918        transition = molName ;
919      }
920      else {
921        transition.resize( numFreq ) ;
922        transition = "" ;
923      }
924      *transitionRF = transition ;
925
926      // SYSVEL
927      Vector<Double> sysvelArr( numFreq, srcVel ) ;
928      *sysvelRF = sysvelArr ;
929
930      //
931      // ITERATION: IFNO
932      //
933      TableIterator iter2( t1, "IFNO" ) ;
934      while( !iter2.pastEnd() ) {
935        //Table t2( iter2.table() ) ;
936        Table t2 = iter2.table() ;
937        uInt nrow = msSrc.nrow() ;
938
939        // get necessary information
940        ROScalarColumn<uInt> ifNoCol( t2, "IFNO" ) ;
941        uInt ifno = ifNoCol( 0 ) ; // IFNO = SPECTRAL_WINDOW_ID
942        Double midTime ;
943        Double interval ;
944        getValidTimeRange( midTime, interval, t2 ) ;
945
946        // fill SPECTRAL_WINDOW_ID
947        *spwidRF = ifno ;
948
949        // fill TIME, INTERVAL
950        *timeRF = midTime ;
951        *intervalRF = interval ;
952
953        // add row
954        msSrc.addRow( 1, True ) ;
955        row.put( nrow ) ;
956
957        iter2.next() ;
958      }
959
960      iter1.next() ;
961    }
962
963    // increment srcId if SRCNAME changed
964    srcId++ ;
965
966    iter0.next() ;
967  }
968
969  double endSec = gettimeofday_sec() ;
970  os_ << "end MSWriter::fillSource() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
971}
972
973void MSWriter::fillWeather()
974{
975  double startSec = gettimeofday_sec() ;
976  os_ << "start MSWriter::fillWeather() startSec=" << startSec << LogIO::POST ;
977
978  // access to MS WEATHER subtable
979  MSWeather msw = mstable_->weather() ;
980
981  // access to WEATHER subtable
982  Table stw = table_->weather().table() ;
983  uInt nrow = stw.nrow() ;
984
985  if ( nrow == 0 )
986    return ;
987
988  msw.addRow( nrow, True ) ;
989  MSWeatherColumns mswCols( msw ) ;
990
991  // ANTENNA_ID is always 0
992  Vector<Int> antIdArr( nrow, 0 ) ;
993  mswCols.antennaId().putColumn( antIdArr ) ;
994
995  // fill weather status
996  ROScalarColumn<Float> sharedFloatCol( stw, "TEMPERATURE" ) ;
997  mswCols.temperature().putColumn( sharedFloatCol ) ;
998  sharedFloatCol.attach( stw, "PRESSURE" ) ;
999  mswCols.pressure().putColumn( sharedFloatCol ) ;
1000  sharedFloatCol.attach( stw, "HUMIDITY" ) ;
1001  mswCols.relHumidity().putColumn( sharedFloatCol ) ;
1002  sharedFloatCol.attach( stw, "WINDSPEED" ) ;
1003  mswCols.windSpeed().putColumn( sharedFloatCol ) ;
1004  sharedFloatCol.attach( stw, "WINDAZ" ) ;
1005  mswCols.windDirection().putColumn( sharedFloatCol ) ;
1006
1007  // fill TIME and INTERVAL
1008  Double midTime ;
1009  Double interval ;
1010  Vector<Double> intervalArr( nrow, 0.0 ) ;
1011  TableIterator iter( table_->table(), "WEATHER_ID" ) ;
1012  while( !iter.pastEnd() ) {
1013    //Table tab( iter.table() ) ;
1014    Table tab = iter.table() ;
1015
1016    ROScalarColumn<uInt> widCol( tab, "WEATHER_ID" ) ;
1017    uInt wid = widCol( 0 ) ;
1018
1019    getValidTimeRange( midTime, interval, tab ) ;
1020    mswCols.time().put( wid, midTime ) ;
1021    intervalArr[wid] = interval ;
1022
1023    iter.next() ;
1024  }
1025  mswCols.interval().putColumn( intervalArr ) ;
1026
1027  double endSec = gettimeofday_sec() ;
1028  os_ << "end MSWriter::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1029}
1030
1031void MSWriter::fillSysCal()
1032{
1033  double startSec = gettimeofday_sec() ;
1034  os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
1035
1036  //tcalIdRec_.print( cout ) ;
1037
1038  // access to MS SYSCAL subtable
1039  MSSysCal mssc = mstable_->sysCal() ;
1040
1041  // access to TCAL subtable
1042  Table stt = table_->tcal().table() ;
1043  uInt nrow = stt.nrow() ;
1044
1045  // access to MAIN table
1046  Block<String> cols( 6 ) ;
1047  cols[0] = "TIME" ;
1048  cols[1] = "TCAL_ID" ;
1049  cols[2] = "TSYS" ;
1050  cols[3] = "BEAMNO" ;
1051  cols[4] = "IFNO" ;
1052  cols[5] = "INTERVAL" ;
1053  Table tab = table_->table().project( cols ) ;
1054
1055  if ( nrow == 0 )
1056    return ;
1057
1058  nrow = tcalIdRec_.nfields() ;
1059
1060  Double midTime ;
1061  Double interval ;
1062  String timeStr ;
1063
1064  // row base
1065  TableRow row( mssc ) ;
1066  TableRecord &trec = row.record() ;
1067  RecordFieldPtr<Int> antennaRF( trec, "ANTENNA_ID" ) ;
1068  RecordFieldPtr<Int> feedRF( trec, "FEED_ID" ) ;
1069  RecordFieldPtr<Int> spwRF( trec, "SPECTRAL_WINDOW_ID" ) ;
1070  RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
1071  RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
1072  RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
1073  RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
1074  RecordFieldPtr< Array<Float> > tsysspRF ;
1075  RecordFieldPtr< Array<Float> > tcalspRF ;
1076  if ( tsysSpec_ )
1077    tsysspRF.attachToRecord( trec, "TSYS_SPECTRUM" ) ;
1078  if ( tcalSpec_ )
1079    tcalspRF.attachToRecord( trec, "TCAL_SPECTRUM" ) ;
1080
1081  // ANTENNA_ID is always 0
1082  *antennaRF = 0 ;
1083
1084  Table sortedstt = stt.sort( "ID" ) ;
1085  ROArrayColumn<Float> tcalCol( sortedstt, "TCAL" ) ;
1086  ROTableColumn idCol( sortedstt, "ID" ) ;
1087  ROArrayColumn<Float> tsysCol( tab, "TSYS" ) ;
1088  ROTableColumn tcalidCol( tab, "TCAL_ID" ) ;
1089  ROTableColumn timeCol( tab, "TIME" ) ;
1090  ROTableColumn intervalCol( tab, "INTERVAL" ) ;
1091  ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
1092  ROTableColumn ifnoCol( tab, "IFNO" ) ;
1093  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1094    double t1 = gettimeofday_sec() ;
1095    Vector<uInt> ids = tcalIdRec_.asArrayuInt( irow ) ;
1096    os_ << "ids = " << ids << LogIO::POST ;
1097    uInt npol = ids.size() ;
1098    Vector<uInt> rows = tcalRowRec_.asArrayuInt( irow ) ;
1099    os_ << "rows = " << rows << LogIO::POST ;
1100    Vector<Double> atime( rows.nelements() ) ;
1101    Vector<Double> ainterval( rows.nelements() ) ;
1102    Vector<uInt> atcalid( rows.nelements() ) ;
1103    for( uInt jrow = 0 ; jrow < rows.nelements() ; jrow++ ) {
1104      atime[jrow] = (Double)timeCol.asdouble( rows[jrow] ) ;
1105      ainterval[jrow] = (Double)intervalCol.asdouble( rows[jrow] ) ;
1106      atcalid[jrow] = tcalidCol.asuInt( rows[jrow] ) ;
1107    }
1108    Vector<Float> dummy = tsysCol( rows[0] ) ;
1109    Matrix<Float> tsys( npol,dummy.nelements() ) ;
1110    tsys.row( 0 ) = dummy ;
1111    for ( uInt jrow = 1 ; jrow < npol ; jrow++ )
1112      tsys.row( jrow ) = tsysCol( rows[jrow] ) ;
1113
1114    // FEED_ID
1115    *feedRF = beamnoCol.asuInt( rows[0] ) ;
1116
1117    // SPECTRAL_WINDOW_ID
1118    *spwRF = ifnoCol.asuInt( rows[0] ) ;
1119
1120    // TIME and INTERVAL
1121    getValidTimeRange( midTime, interval, atime, ainterval ) ;
1122    *timeRF = midTime ;
1123    *intervalRF = interval ;
1124
1125    // TCAL and TSYS
1126    Matrix<Float> tcal ;
1127    Table t ;
1128    if ( idCol.asuInt( ids[0] ) == ids[0] ) {
1129      os_ << "sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
1130      Vector<Float> dummyC = tcalCol( ids[0] ) ;
1131      tcal.resize( npol, dummyC.size() ) ;
1132      tcal.row( 0 ) = dummyC ;
1133    }
1134    else {
1135      os_ << "NOT sorted at irow=" << irow << " ids[0]=" << ids[0] << LogIO::POST ;
1136      t = stt( stt.col("ID") == ids[0] ) ;
1137      Vector<Float> dummyC = tcalCol( 0 ) ;
1138      tcal.resize( npol, dummyC.size() ) ;
1139      tcal.row( 0 ) = dummyC ;
1140    }
1141    if ( npol == 2 ) {
1142      if ( idCol.asuInt( ids[1] ) == ids[1] ) {
1143        os_ << "sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
1144        tcal.row( 1 ) = tcalCol( ids[1] ) ;
1145      }
1146      else {
1147        os_ << "NOT sorted at irow=" << irow << " ids[1]=" << ids[1] << LogIO::POST ;
1148        t = stt( stt.col("ID") == ids[1] ) ;
1149        tcalCol.attach( t, "TCAL" ) ;
1150        tcal.row( 1 ) = tcalCol( 1 ) ;
1151      }
1152    }
1153    else if ( npol == 3 ) {
1154      if ( idCol.asuInt( ids[2] ) == ids[2] )
1155        tcal.row( 1 ) = tcalCol( ids[2] ) ;
1156      else {
1157        t = stt( stt.col("ID") == ids[2] ) ;
1158        tcalCol.attach( t, "TCAL" ) ;
1159        tcal.row( 1 ) = tcalCol( 0 ) ;
1160      }
1161      if ( idCol.asuInt( ids[1] ) == ids[1] )
1162        tcal.row( 2 ) = tcalCol( ids[1] ) ;
1163      else {
1164        t = stt( stt.col("ID") == ids[1] ) ;
1165        tcalCol.attach( t, "TCAL" ) ;
1166        tcal.row( 2 ) = tcalCol( 0 ) ;
1167      }
1168    }
1169    else if ( npol == 4 ) {
1170      if ( idCol.asuInt( ids[2] ) == ids[2] )
1171        tcal.row( 1 ) = tcalCol( ids[2] ) ;
1172      else {
1173        t = stt( stt.col("ID") == ids[2] ) ;
1174        tcalCol.attach( t, "TCAL" ) ;
1175        tcal.row( 1 ) = tcalCol( 0 ) ;
1176      }
1177      if ( idCol.asuInt( ids[3] ) == ids[3] )
1178        tcal.row( 2 ) = tcalCol( ids[3] ) ;
1179      else {
1180        t = stt( stt.col("ID") == ids[3] ) ;
1181        tcalCol.attach( t, "TCAL" ) ;
1182        tcal.row( 2 ) = tcalCol( 0 ) ;
1183      }
1184      if ( idCol.asuInt( ids[1] ) == ids[1] )
1185        tcal.row( 2 ) = tcalCol( ids[1] ) ;
1186      else {
1187        t = stt( stt.col("ID") == ids[1] ) ;
1188        tcalCol.attach( t, "TCAL" ) ;
1189        tcal.row( 3 ) = tcalCol( 0 ) ;
1190      }
1191    }
1192    if ( tcalSpec_ ) {
1193      // put TCAL_SPECTRUM
1194      //*tcalspRF = tcal ;
1195      tcalspRF.define( tcal ) ;
1196      // set TCAL (mean of TCAL_SPECTRUM)
1197      Matrix<Float> tcalMean( npol, 1 ) ;
1198      for ( uInt iid = 0 ; iid < npol ; iid++ ) {
1199        tcalMean( iid, 0 ) = mean( tcal.row(iid) ) ;
1200      }
1201      // put TCAL
1202      *tcalRF = tcalMean ;
1203    }
1204    else {
1205      // put TCAL
1206      *tcalRF = tcal ;
1207    }
1208   
1209    if ( tsysSpec_ ) {
1210      // put TSYS_SPECTRUM
1211      //*tsysspRF = tsys ;
1212      tsysspRF.define( tsys ) ;
1213      // set TSYS (mean of TSYS_SPECTRUM)
1214      Matrix<Float> tsysMean( npol, 1 ) ;
1215      for ( uInt iid = 0 ; iid < npol ; iid++ ) {
1216        tsysMean( iid, 0 ) = mean( tsys.row(iid) ) ;
1217      }
1218      // put TSYS
1219      *tsysRF = tsysMean ;
1220    }
1221    else {
1222      // put TSYS
1223      *tsysRF = tsys ;
1224    }
1225
1226    // add row
1227    mssc.addRow( 1, True ) ;
1228    row.put( mssc.nrow()-1 ) ;
1229
1230    double t2 = gettimeofday_sec() ;
1231    os_ << irow << "th loop elapsed time = " << t2-t1 << "sec" << LogIO::POST ;
1232  }
1233 
1234  double endSec = gettimeofday_sec() ;
1235  os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1236}
1237
1238void MSWriter::addFeed( Int id )
1239{
1240  double startSec = gettimeofday_sec() ;
1241  os_ << "start MSWriter::addFeed() startSec=" << startSec << LogIO::POST ;
1242
1243  // add row
1244  MSFeed msFeed = mstable_->feed() ;
1245  msFeed.addRow( 1, True ) ;
1246  Int nrow = msFeed.nrow() ;
1247  Int numReceptors = 2 ;
1248  Vector<String> polType( numReceptors ) ;
1249  Matrix<Double> beamOffset( 2, numReceptors ) ;
1250  beamOffset = 0.0 ;
1251  Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
1252  if ( polType_ == "linear" ) {
1253    polType[0] = "X" ;
1254    polType[1] = "Y" ;
1255  }
1256  else if ( polType_ == "circular" ) {
1257    polType[0] = "R" ;
1258    polType[1] = "L" ;
1259  }
1260  else {
1261    polType[0] = "X" ;
1262    polType[1] = "Y" ;
1263  }
1264  Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
1265  for ( Int i = 0 ; i < numReceptors ; i++ )
1266    polResponse( i, i ) = 0.0 ;
1267
1268  MSFeedColumns msFeedCols( mstable_->feed() ) ;
1269
1270  msFeedCols.feedId().put( nrow-1, id ) ;
1271  msFeedCols.antennaId().put( nrow-1, 0 ) ;
1272  msFeedCols.numReceptors().put( nrow-1, numReceptors ) ;
1273  msFeedCols.polarizationType().put( nrow-1, polType ) ;
1274  msFeedCols.beamOffset().put( nrow-1, beamOffset ) ;
1275  msFeedCols.receptorAngle().put( nrow-1, receptorAngle ) ;
1276  msFeedCols.polResponse().put( nrow-1, polResponse ) ;
1277
1278  double endSec = gettimeofday_sec() ;
1279  os_ << "end MSWriter::addFeed() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1280}
1281
1282void MSWriter::addSpectralWindow( Int spwid, Int freqid )
1283{
1284  double startSec = gettimeofday_sec() ;
1285  os_ << "start MSWriter::addSpectralWindow() startSec=" << startSec << LogIO::POST ;
1286 
1287  // add row
1288  MSSpectralWindow msSpw = mstable_->spectralWindow() ;
1289  while( (Int)msSpw.nrow() <= spwid ) {
1290    msSpw.addRow( 1, True ) ;
1291  }
1292 
1293  MSSpWindowColumns msSpwCols( msSpw ) ;
1294
1295  STFrequencies stf = table_->frequencies() ;
1296
1297  // MEAS_FREQ_REF
1298  msSpwCols.measFreqRef().put( spwid, stf.getFrame( True ) ) ;
1299
1300  Double refpix ;
1301  Double refval ;
1302  Double inc ;
1303  stf.getEntry( refpix, refval, inc, (uInt)freqid ) ;
1304
1305  // NUM_CHAN
1306  Int nchan = (Int)(refpix * 2) ;
1307  msSpwCols.numChan().put( spwid, nchan ) ;
1308
1309  // TOTAL_BANDWIDTH
1310  Double bw = nchan * inc ;
1311  msSpwCols.totalBandwidth().put( spwid, bw ) ;
1312
1313  // REF_FREQUENCY
1314  Double refFreq = refval - refpix * inc ;
1315  msSpwCols.refFrequency().put( spwid, refFreq ) ;
1316
1317  // NET_SIDEBAND
1318  // tentative: USB->0, LSB->1
1319  Int netSideband = 0 ;
1320  if ( inc < 0 )
1321    netSideband = 1 ;
1322  msSpwCols.netSideband().put( spwid, netSideband ) ;
1323
1324  // RESOLUTION, CHAN_WIDTH, EFFECTIVE_BW
1325  Vector<Double> sharedDoubleArr( nchan, inc ) ;
1326  msSpwCols.resolution().put( spwid, sharedDoubleArr ) ;
1327  msSpwCols.chanWidth().put( spwid, sharedDoubleArr ) ;
1328  msSpwCols.effectiveBW().put( spwid, sharedDoubleArr ) ;
1329
1330  // CHAN_FREQ
1331  indgen( sharedDoubleArr, refFreq, inc ) ;
1332  msSpwCols.chanFreq().put( spwid, sharedDoubleArr ) ;
1333
1334  double endSec = gettimeofday_sec() ;
1335  os_ << "end MSWriter::addSpectralWindow() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1336}
1337
1338void MSWriter::addField( Int fid, String fieldname, String srcname, Double t, Vector<Double> rate )
1339{
1340  double startSec = gettimeofday_sec() ;
1341  os_ << "start MSWriter::addField() startSec=" << startSec << LogIO::POST ;
1342 
1343  MSField msField = mstable_->field() ;
1344  while( (Int)msField.nrow() <= fid ) {
1345    msField.addRow( 1, True ) ;
1346  }
1347  MSFieldColumns msFieldCols( msField ) ;
1348
1349  // Access to SOURCE table
1350  MSSource msSrc = mstable_->source() ;
1351
1352  // fill target row
1353  msFieldCols.name().put( fid, fieldname ) ;
1354  msFieldCols.time().put( fid, t ) ;
1355  Int numPoly = 0 ;
1356  if ( anyNE( rate, 0.0 ) )
1357    numPoly = 1 ;
1358  msFieldCols.numPoly().put( fid, numPoly ) ;
1359  MSSourceIndex msSrcIdx( msSrc ) ;
1360  Int srcId = -1 ;
1361  Vector<Int> srcIdArr = msSrcIdx.matchSourceName( srcname ) ;
1362  if ( srcIdArr.size() != 0 ) {
1363    srcId = srcIdArr[0] ;
1364    MSSource msSrcSel = msSrc( msSrc.col("SOURCE_ID") == srcId ) ;
1365    ROMSSourceColumns msSrcCols( msSrcSel ) ;
1366    Vector<Double> srcDir = msSrcCols.direction()( 0 ) ;
1367    Matrix<Double> srcDirA( IPosition( 2, 2, 1+numPoly ) ) ;
1368    os_ << "srcDirA = " << srcDirA << LogIO::POST ;
1369    os_ << "sliced srcDirA = " << srcDirA.column( 0 ) << LogIO::POST ;
1370    srcDirA.column( 0 ) = srcDir ;
1371    os_ << "srcDirA = " << srcDirA << LogIO::POST ;
1372    if ( numPoly != 0 )
1373      srcDirA.column( 1 ) = rate ;
1374    msFieldCols.phaseDir().put( fid, srcDirA ) ;
1375    msFieldCols.referenceDir().put( fid, srcDirA ) ;
1376    msFieldCols.delayDir().put( fid, srcDirA ) ;
1377  }
1378  msFieldCols.sourceId().put( fid, srcId ) ;
1379
1380  double endSec = gettimeofday_sec() ;
1381  os_ << "end MSWriter::addField() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1382}
1383
1384void MSWriter::addPointing( String &name, Double &me, Double &interval, Matrix<Double> &dir )
1385{
1386  double startSec = gettimeofday_sec() ;
1387  os_ << "start MSWriter::addPointing() startSec=" << startSec << LogIO::POST ;
1388 
1389  // access to POINTING subtable
1390  MSPointing msp = mstable_->pointing() ;
1391  uInt nrow = msp.nrow() ;
1392
1393  // add row
1394  msp.addRow( 1, True ) ;
1395
1396  // fill row
1397  TableRow row( msp ) ;
1398  TableRecord &rec = row.record() ;
1399  RecordFieldPtr<Int> antennaRF( rec, "ANTENNA_ID" ) ;
1400  *antennaRF = 0 ;
1401  RecordFieldPtr<Int> numpolyRF( rec, "NUM_POLY" ) ;
1402  *numpolyRF = dir.ncolumn() ;
1403  RecordFieldPtr<Double> timeRF( rec, "TIME" ) ;
1404  *timeRF = me ;
1405  RecordFieldPtr<Double> toriginRF( rec, "TIME_ORIGIN" ) ;
1406  *toriginRF = me ;
1407  RecordFieldPtr<Double> intervalRF( rec, "INTERVAL" ) ;
1408  *intervalRF = interval ;
1409  RecordFieldPtr<String> nameRF( rec, "NAME" ) ;
1410  *nameRF = name ;
1411  RecordFieldPtr<Bool> trackRF( rec, "TRACKING" ) ;
1412  *trackRF = True ;
1413  RecordFieldPtr< Array<Double> > dirRF( rec, "DIRECTION" ) ;
1414  *dirRF = dir ;
1415  RecordFieldPtr< Array<Double> > targetRF( rec, "TARGET" ) ;
1416  *dirRF = dir ;
1417  row.put( nrow ) ;
1418
1419  double endSec = gettimeofday_sec() ;
1420  os_ << "end MSWriter::addPointing() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1421}
1422
1423Int MSWriter::addPolarization( Vector<Int> polnos )
1424{
1425  double startSec = gettimeofday_sec() ;
1426  os_ << "start MSWriter::addPolarization() startSec=" << startSec << LogIO::POST ;
1427
1428  os_ << "polnos = " << polnos << LogIO::POST ;
1429  MSPolarization msPol = mstable_->polarization() ;
1430  uInt nrow = msPol.nrow() ;
1431
1432  // only 1 POLARIZATION row for 1 scantable
1433  if ( nrow > 0 )
1434    return 0 ;
1435 
1436  Vector<Int> corrType = toCorrType( polnos ) ;
1437 
1438  ROArrayColumn<Int> corrtCol( msPol, "CORR_TYPE" ) ;
1439  Matrix<Int> corrTypeArr = corrtCol.getColumn() ;
1440  Int polid = -1 ;
1441  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1442    if ( allEQ( corrType, corrTypeArr.column( irow ) ) ) {
1443      polid = irow ;
1444      break ;
1445    }
1446  }
1447
1448  if ( polid == -1 ) {
1449    MSPolarizationColumns msPolCols( msPol ) ;
1450
1451    // add row
1452    msPol.addRow( 1, True ) ;
1453    polid = (Int)nrow ;
1454
1455    // CORR_TYPE
1456    msPolCols.corrType().put( nrow, corrType ) ;
1457
1458    // NUM_CORR
1459    uInt npol = corrType.size() ;
1460    msPolCols.numCorr().put( nrow, npol ) ;
1461
1462    // CORR_PRODUCT
1463    Matrix<Int> corrProd( 2, npol, -1 ) ;
1464    if ( npol == 1 ) {
1465      corrProd = 0 ;
1466    }
1467    else if ( npol == 2 ) {
1468      corrProd.column( 0 ) = 0 ;
1469      corrProd.column( 1 ) = 1 ;
1470    }
1471    else {
1472      corrProd.column( 0 ) = 0 ;
1473      corrProd.column( 3 ) = 1 ;
1474      corrProd( 0,1 ) = 0 ;
1475      corrProd( 1,1 ) = 1 ;
1476      corrProd( 0,2 ) = 1 ;
1477      corrProd( 1,2 ) = 0 ;
1478    }
1479    msPolCols.corrProduct().put( nrow, corrProd ) ;   
1480  }
1481
1482  double endSec = gettimeofday_sec() ;
1483  os_ << "end MSWriter::addPolarization() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1484
1485  return polid ;
1486}
1487
1488Int MSWriter::addDataDescription( Int polid, Int spwid )
1489{
1490  double startSec = gettimeofday_sec() ;
1491  os_ << "start MSWriter::addDataDescription() startSec=" << startSec << LogIO::POST ;
1492
1493  MSDataDescription msDataDesc = mstable_->dataDescription() ;
1494  uInt nrow = msDataDesc.nrow() ;
1495
1496  // only 1 POLARIZATION_ID for 1 scantable
1497  Int ddid = -1 ;
1498  ROScalarColumn<Int> spwCol( msDataDesc, "SPECTRAL_WINDOW_ID" ) ;
1499  Vector<Int> spwIds = spwCol.getColumn() ;
1500  //ROScalarColumn<Int> polCol( msDataDesc, "POLARIZATION_ID" ) ;
1501  //Vector<Int> polIds = polCol.getColumn() ;
1502  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1503    //if ( spwid == spwIds[irow] && polid == polIds[irow] ) {
1504    if ( spwid == spwIds[irow] ) {
1505      ddid = irow ;
1506      break ;
1507    }
1508  }
1509  os_ << "ddid = " << ddid << LogIO::POST ;
1510 
1511
1512  if ( ddid == -1 ) {
1513    msDataDesc.addRow( 1, True ) ;
1514    MSDataDescColumns msDataDescCols( msDataDesc ) ;
1515    msDataDescCols.polarizationId().put( nrow, polid ) ;
1516    msDataDescCols.spectralWindowId().put( nrow, spwid ) ;
1517    ddid = (Int)nrow ;
1518  }
1519
1520  double endSec = gettimeofday_sec() ;
1521  os_ << "end MSWriter::addDataDescription() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1522
1523  return ddid ;
1524}
1525
1526Int MSWriter::addState( Int st, Int &subscan )
1527{
1528  double startSec = gettimeofday_sec() ;
1529  os_ << "start MSWriter::addState() startSec=" << startSec << LogIO::POST ;
1530
1531  // access to STATE subtable
1532  MSState msState = mstable_->state() ;
1533  uInt nrow = msState.nrow() ;
1534
1535  String obsMode ;
1536  Bool isSignal ;
1537  Double tnoise ;
1538  Double tload ;
1539  queryType( st, obsMode, isSignal, tnoise, tload ) ;
1540  os_ << "obsMode = " << obsMode << " isSignal = " << isSignal << LogIO::POST ;
1541
1542  Int idx = -1 ;
1543  ROScalarColumn<String> obsModeCol( msState, "OBS_MODE" ) ;
1544  //ROScalarColumn<Bool> sigCol( msState, "SIG" ) ;
1545  //ROScalarColumn<Bool> refCol( msState, "REF" ) ;
1546  ROScalarColumn<Int> subscanCol( msState, "SUB_SCAN" ) ;
1547//   Vector<String> obsModeArr = obsModeCol.getColumn() ;
1548//   Vector<Bool> sigArr = sigCol.getColumn() ;
1549//   Vector<Bool> refArr = refCol.getColumn() ;
1550//   Vector<Int> subscanArr = subscanCol.getColumn() ;
1551  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1552    if ( obsModeCol(irow) == obsMode
1553         //&& sigCol(irow) == isSignal
1554         //&& refCol(irow) != isSignal
1555         && subscanCol(irow) == subscan ) {
1556//     if ( obsModeArr[irow] == obsMode
1557//          && sigArr[irow] == isSignal
1558//          && refArr[irow] != isSignal
1559//          && subscanArr[irow] == subscan ) {
1560      idx = irow ;
1561      break ;
1562    }
1563  }
1564  if ( idx == -1 ) {
1565    msState.addRow( 1, True ) ;
1566    TableRow row( msState ) ;
1567    TableRecord &rec = row.record() ;
1568    RecordFieldPtr<String> obsmodeRF( rec, "OBS_MODE" ) ;
1569    *obsmodeRF = obsMode ;
1570    RecordFieldPtr<Bool> sigRF( rec, "SIG" ) ;
1571    *sigRF = isSignal ;
1572    RecordFieldPtr<Bool> refRF( rec, "REF" ) ;
1573    *refRF = !isSignal ;
1574    RecordFieldPtr<Int> subscanRF( rec, "SUB_SCAN" ) ;
1575    *subscanRF = subscan ;
1576    RecordFieldPtr<Double> noiseRF( rec, "CAL" ) ;
1577    *noiseRF = tnoise ;
1578    RecordFieldPtr<Double> loadRF( rec, "LOAD" ) ;
1579    *loadRF = tload ;
1580    row.put( nrow ) ;
1581//     ScalarColumn<String> obsModeCol( msState, "OBS_MODE" ) ;
1582//     obsModeCol.put( nrow, obsMode ) ;
1583//     ScalarColumn<Bool> sharedBCol( msState, "SIG" ) ;
1584//     sharedBCol.put( nrow, isSignal ) ;
1585//     sharedBCol.attach( msState, "REF" ) ;
1586//     sharedBCol.put( nrow, !isSignal ) ;
1587//     ScalarColumn<Int> subscanCol( msState, "SUB_SCAN" ) ;
1588//     subscanCol.put( nrow, subscan ) ;
1589    idx = nrow ;
1590  }
1591  subscan++ ;
1592
1593  double endSec = gettimeofday_sec() ;
1594  os_ << "end MSWriter::addState() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1595
1596  return idx ;
1597}
1598
1599Vector<Int> MSWriter::toCorrType( Vector<Int> polnos )
1600{
1601  double startSec = gettimeofday_sec() ;
1602  os_ << "start MSWriter::toCorrType() startSec=" << startSec << LogIO::POST ;
1603
1604  uInt npol = polnos.size() ;
1605  Vector<Int> corrType( npol, Stokes::Undefined ) ;
1606 
1607  if ( npol == 4 ) {
1608    if ( polType_ == "linear" ) {
1609      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1610        if ( polnos[ipol] == 0 )
1611          corrType[ipol] = Stokes::XX ;
1612        else if ( polnos[ipol] == 1 )
1613          corrType[ipol] = Stokes::XY ;
1614        else if ( polnos[ipol] == 2 )
1615          corrType[ipol] = Stokes::YX ;
1616        else if ( polnos[ipol] == 3 )
1617          corrType[ipol] = Stokes::YY ;
1618      }
1619    }
1620    else if ( polType_ == "circular" ) {
1621      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1622        if ( polnos[ipol] == 0 )
1623          corrType[ipol] = Stokes::RR ;
1624        else if ( polnos[ipol] == 1 )
1625          corrType[ipol] = Stokes::RL ;
1626        else if ( polnos[ipol] == 2 )
1627          corrType[ipol] = Stokes::LR ;
1628        else if ( polnos[ipol] == 3 )
1629          corrType[ipol] = Stokes::LL ;
1630      }
1631    }
1632    else if ( polType_ == "stokes" ) {
1633      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1634        if ( polnos[ipol] == 0 )
1635          corrType[ipol] = Stokes::I ;
1636        else if ( polnos[ipol] == 1 )
1637          corrType[ipol] = Stokes::Q ;
1638        else if ( polnos[ipol] == 2 )
1639          corrType[ipol] = Stokes::U ;
1640        else if ( polnos[ipol] == 3 )
1641          corrType[ipol] = Stokes::V ;
1642      }
1643    }
1644  }
1645  else if ( npol == 2 ) {
1646    if ( polType_ == "linear" ) {
1647      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1648        if ( polnos[ipol] == 0 )
1649          corrType[ipol] = Stokes::XX ;
1650        else if ( polnos[ipol] == 1 )
1651          corrType[ipol] = Stokes::YY ;
1652      }
1653    }
1654    else if ( polType_ == "circular" ) {
1655      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1656        if ( polnos[ipol] == 0 )
1657          corrType[ipol] = Stokes::RR ;
1658        else if ( polnos[ipol] == 1 )
1659          corrType[ipol] = Stokes::LL ;
1660      }
1661    }
1662    else if ( polType_ == "stokes" ) {
1663      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1664        if ( polnos[ipol] == 0 )
1665          corrType[ipol] = Stokes::I ;
1666        else if ( polnos[ipol] == 1 )
1667          corrType[ipol] = Stokes::V ;
1668      }
1669    }
1670    else if ( polType_ == "linpol" ) {
1671      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1672        if ( polnos[ipol] == 1 )
1673          corrType[ipol] = Stokes::Plinear ;
1674        else if ( polnos[ipol] == 2 )
1675          corrType[ipol] = Stokes::Pangle ;
1676      }
1677    }
1678  }     
1679  else if ( npol == 1 ) {
1680    if ( polType_ == "linear" )
1681      corrType[0] = Stokes::XX ;
1682    else if ( polType_ == "circular" )
1683      corrType[0] = Stokes::RR ;
1684    else if ( polType_ == "stokes" )
1685      corrType[0] = Stokes::I ;
1686  }
1687
1688  double endSec = gettimeofday_sec() ;
1689  os_ << "end MSWriter::toCorrType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1690
1691  return corrType ;
1692}
1693
1694void MSWriter::getValidTimeRange( Double &me, Double &interval, Table &tab )
1695{
1696  double startSec = gettimeofday_sec() ;
1697  os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
1698
1699  // sort table
1700  //Table stab = tab.sort( "TIME" ) ;
1701
1702  ROScalarColumn<Double> timeCol( tab, "TIME" ) ;
1703  Vector<Double> timeArr = timeCol.getColumn() ;
1704  Double minTime ;
1705  Double maxTime ;
1706  minMax( minTime, maxTime, timeArr ) ;
1707  Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
1708  // unit for TIME
1709  // Scantable: "d"
1710  // MS: "s"
1711  me = midTime ;
1712  interval = ( maxTime - minTime ) * 86400.0 ;
1713
1714  double endSec = gettimeofday_sec() ;
1715  os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1716}
1717
1718void MSWriter::getValidTimeRange( Double &me, Double &interval, Vector<Double> &atime, Vector<Double> &ainterval )
1719{
1720  double startSec = gettimeofday_sec() ;
1721  os_ << "start MSWriter::getVaridTimeRange() startSec=" << startSec << LogIO::POST ;
1722
1723  // sort table
1724  //Table stab = tab.sort( "TIME" ) ;
1725
1726  Double minTime ;
1727  Double maxTime ;
1728  minMax( minTime, maxTime, atime ) ;
1729  Double midTime = 0.5 * ( minTime + maxTime ) * 86400.0 ;
1730  // unit for TIME
1731  // Scantable: "d"
1732  // MS: "s"
1733  me = midTime ;
1734  interval = ( maxTime - minTime ) * 86400.0 + mean( ainterval ) ;
1735
1736  double endSec = gettimeofday_sec() ;
1737  os_ << "end MSWriter::getValidTimeRange() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1738}
1739
1740//void MSWriter::queryType( Int type, String &stype, Bool &b )
1741void MSWriter::queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
1742{
1743  double startSec = gettimeofday_sec() ;
1744  os_ << "start MSWriter::queryType() startSec=" << startSec << LogIO::POST ;
1745
1746  switch ( type ) {
1747  case SrcType::PSON:
1748    stype = "POSITION_SWITCH.OBSERVE_TARGET.ON_SOURCE" ;
1749    b = True ;
1750    t = 0.0 ;
1751    l = 0.0 ;
1752    break ;
1753  case SrcType::PSOFF:
1754    stype = "POSITION_SWITCH.OBSERVE_TARGET.OFF_SOURCE" ;
1755    b = False ;
1756    t = 0.0 ;
1757    l = 0.0 ;
1758    break ;
1759  case SrcType::NOD:
1760    stype = "NOD.OBSERVE_TARGET.ON_SOURCE" ;
1761    b = True ;
1762    t = 0.0 ;
1763    l = 0.0 ;
1764    break ;
1765  case SrcType::FSON:
1766    stype = "FREQUENCY_SWITCH.OBSERVE_TARGET.ON_SOURCE" ;
1767    b = True ;
1768    t = 0.0 ;
1769    l = 0.0 ;
1770    break ;
1771  case SrcType::FSOFF:
1772    stype = "FREQUENCY_SWITCH.OBSERVE_TARGET.ON_SOURCE" ;
1773    b = False ;
1774    t = 0.0 ;
1775    l = 0.0 ;
1776    break ;
1777  case SrcType::SKY:
1778    stype = "UNSPECIFIED.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1779    b = False ;
1780    t = 0.0 ;
1781    l = 1.0 ;
1782    break ;
1783  case SrcType::HOT:
1784    stype = "UNSPECIFIED.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1785    b = False ;
1786    t = 0.0 ;
1787    l = 1.0 ;
1788    break ;
1789  case SrcType::WARM:
1790    stype = "UNSPECIFIED.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1791    t = 0.0 ;
1792    b = False ;
1793    l = 1.0 ;
1794    break ;
1795  case SrcType::COLD:
1796    stype = "UNSPECIFIED.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1797    b = False ;
1798    t = 0.0 ;
1799    l = 1.0 ;
1800    break ;
1801  case SrcType::PONCAL:
1802    stype = "POSITION_SWITCH.CALIBRATE_TEMPERATURE.ON_SOURCE" ;
1803    b = True ;
1804    t = 1.0 ;
1805    l = 0.0 ;
1806    break ;
1807  case SrcType::POFFCAL:
1808    stype = "POSITION_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1809    b = False ;
1810    t = 1.0 ;
1811    l = 0.0 ;
1812    break ;
1813  case SrcType::NODCAL:
1814    stype = "NOD.CALIBRATE_TEMPERATURE.ON_SOURCE" ;
1815    b = True ;
1816    t = 1.0 ;
1817    l = 0.0 ;
1818    break ;
1819  case SrcType::FONCAL:
1820    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.ON_SOURCE" ;
1821    b = True ;
1822    t = 1.0 ;
1823    l = 0.0 ;
1824    break ;
1825  case SrcType::FOFFCAL:
1826    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1827    b = False ;
1828    t = 1.0 ;
1829    l = 0.0 ;
1830    break ;
1831  case SrcType::FSLO:
1832    stype = "FREQUENCY_SWITCH.OBSERVE_TARGET.ON_SOURCE" ;
1833    b = True ;
1834    t = 0.0 ;
1835    l = 0.0 ;
1836    break ;
1837  case SrcType::FLOOFF:
1838    stype = "FREQUENCY_SWITCH.OBSERVE_TARGET.OFF_SOURCE" ;
1839    b = False ;
1840    t = 0.0 ;
1841    l = 0.0 ;
1842    break ;
1843  case SrcType::FLOSKY:
1844    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1845    b = False ;
1846    t = 0.0 ;
1847    l = 1.0 ;
1848    break ;
1849  case SrcType::FLOHOT:
1850    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1851    b = False ;
1852    t = 0.0 ;
1853    l = 1.0 ;
1854    break ;
1855  case SrcType::FLOWARM:
1856    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1857    b = False ;
1858    t = 0.0 ;
1859    l = 1.0 ;
1860    break ;
1861  case SrcType::FLOCOLD:
1862    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1863    b = False ;
1864    t = 0.0 ;
1865    l = 1.0 ;
1866    break ;
1867  case SrcType::FSHI:
1868    stype = "FREQUENCY_SWITCH.OBSERVE_TARGET.ON_SOURCE" ;
1869    b = True ;
1870    t = 0.0 ;
1871    l = 0.0 ;
1872    break ;
1873  case SrcType::FHIOFF:
1874    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1875    b = False ;
1876    t = 0.0 ;
1877    l = 0.0 ;
1878    break ;
1879  case SrcType::FHISKY:
1880    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1881    b = False ;
1882    t = 0.0 ;
1883    l = 1.0 ;
1884    break ;
1885  case SrcType::FHIHOT:
1886    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1887    b = False ;
1888    t = 0.0 ;
1889    l = 1.0 ;
1890    break ;
1891  case SrcType::FHIWARM:
1892    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1893    b = False ;
1894    t = 0.0 ;
1895    l = 1.0 ;
1896    break ;
1897  case SrcType::FHICOLD:
1898    stype = "FREQUENCY_SWITCH.CALIBRATE_TEMPERATURE.OFF_SOURCE" ;
1899    b = False ;
1900    t = 0.0 ;
1901    l = 1.0 ;
1902    break ;
1903  case SrcType::SIG:
1904    stype = "UNSPECIFIED.OBSERVE_TARGET.ON_SOURCE" ;
1905    b = True ;
1906    t = 0.0 ;
1907    l = 0.0 ;
1908    break ;
1909  case SrcType::REF:
1910    stype = "UNSPECIFIED.OBSERVE_TARGET.ON_SOURCE" ;
1911    b = False ;
1912    t = 0.0 ;
1913    l = 0.0 ;
1914    break ;
1915  default:
1916    stype = "UNSPECIFIED" ;
1917    b = True ;
1918    t = 0.0 ;
1919    l = 0.0 ;
1920    break ;
1921  }
1922
1923  double endSec = gettimeofday_sec() ;
1924  os_ << "end MSWriter::queryType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1925}
1926}
Note: See TracBrowser for help on using the repository browser.