source: branches/casa-prerelease/pre-asap/src/MSWriter.cpp @ 2191

Last change on this file since 2191 was 2191, checked in by Kana Sugimoto, 13 years ago

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