source: branches/casa-prerelease/pre-asap/src/MSFiller.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.2 KB
Line 
1//
2// C++ Interface: MSFiller
3//
4// Description:
5//
6// This class is specific filler 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 <iostream>
15#include <map>
16
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/TableIter.h>
19#include <tables/Tables/TableColumn.h>
20#include <tables/Tables/ScalarColumn.h>
21#include <tables/Tables/ArrayColumn.h>
22#include <tables/Tables/TableParse.h>
23#include <tables/Tables/TableRow.h>
24
25#include <casa/Containers/RecordField.h>
26#include <casa/Logging/LogIO.h>
27#include <casa/Arrays/Slicer.h>
28#include <casa/Quanta/MVTime.h>
29#include <casa/OS/Path.h>
30
31#include <measures/Measures/Stokes.h>
32#include <measures/Measures/MEpoch.h>
33#include <measures/Measures/MCEpoch.h>
34#include <measures/Measures/MFrequency.h>
35#include <measures/Measures/MCFrequency.h>
36#include <measures/Measures/MPosition.h>
37#include <measures/Measures/MCPosition.h>
38#include <measures/Measures/MDirection.h>
39#include <measures/Measures/MCDirection.h>
40#include <measures/Measures/MeasConvert.h>
41#include <measures/TableMeasures/ScalarMeasColumn.h>
42#include <measures/TableMeasures/ArrayMeasColumn.h>
43#include <measures/TableMeasures/ScalarQuantColumn.h>
44#include <measures/TableMeasures/ArrayQuantColumn.h>
45
46#include <ms/MeasurementSets/MSAntennaIndex.h>
47
48#include <atnf/PKSIO/SrcType.h>
49
50#include "MSFiller.h"
51#include "STHeader.h"
52
53#include <ctime>
54#include <sys/time.h>
55
56using namespace casa ;
57using namespace std ;
58
59namespace asap {
60double MSFiller::gettimeofday_sec()
61{
62  struct timeval tv ;
63  gettimeofday( &tv, NULL ) ;
64  return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
65}
66
67MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
68  : table_( stable ),
69    tablename_( "" ),
70    antenna_( -1 ),
71    antennaStr_(""),
72    getPt_( False ),
73    isFloatData_( False ),
74    isData_( False ),
75    isDoppler_( False ),
76    isFlagCmd_( False ),
77    isFreqOffset_( False ),
78    isHistory_( False ),
79    isProcessor_( False ),
80    isSysCal_( False ),
81    isWeather_( False ),
82    colTsys_( "TSYS_SPECTRUM" ),
83    colTcal_( "TCAL_SPECTRUM" )
84{
85  os_ = LogIO() ;
86  os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
87}
88
89MSFiller::~MSFiller()
90{
91  os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
92}
93
94bool MSFiller::open( const std::string &filename, const casa::Record &rec )
95{
96  os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
97//   double startSec = gettimeofday_sec() ;
98//   os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
99  //os_ << "   filename = " << filename << endl ;
100
101  // parsing MS options
102  if ( rec.isDefined( "ms" ) ) {
103    Record msrec = rec.asRecord( "ms" ) ;
104    if ( msrec.isDefined( "getpt" ) ) {
105      getPt_ = msrec.asBool( "getpt" ) ;
106    }
107    if ( msrec.isDefined( "antenna" ) ) {
108      if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
109        antenna_ = msrec.asInt( "antenna" ) ;
110      }
111      else {
112        //antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
113        antennaStr_ = msrec.asString( "antenna" ) ;
114      }
115    }
116    else {
117      antenna_ = 0 ;
118    }
119  }
120
121  MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
122  //mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
123  //                     && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;
124  tablename_ = tmpMS->tableName() ;
125  if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
126    MSAntennaIndex msAntIdx( tmpMS->antenna() ) ;
127    Vector<Int> id = msAntIdx.matchAntennaName( antennaStr_ ) ;
128    if ( id.size() > 0 )
129      antenna_ = id[0] ;
130  }
131
132  os_ << "Parsing MS options" << endl ;
133  os_ << "   getPt = " << getPt_ << endl ;
134  os_ << "   antenna = " << antenna_ << endl ;
135  os_ << "   antennaStr = " << antennaStr_ << LogIO::POST ;
136
137  mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
138                                       && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
139//   stringstream ss ;
140//   ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 == " << antenna_ ;
141//   String taql( ss.str() ) ;
142//   mstable_ = MeasurementSet( tableCommand( taql, *tmpMS ) ) ;
143  delete tmpMS ;
144
145  // check which data column exists
146  isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
147  isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
148
149//   double endSec = gettimeofday_sec() ;
150//   os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
151  return true ;
152}
153
154void MSFiller::fill()
155{
156  os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
157//   double startSec = gettimeofday_sec() ;
158//   os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
159
160//   double time0 = gettimeofday_sec() ;
161//   os_ << "start init fill: " << time0 << LogIO::POST ;
162
163  // Initialize header
164  STHeader sdh ; 
165  sdh.nchan = 0 ;
166  sdh.npol = 0 ;
167  sdh.nif = 0 ;
168  sdh.nbeam = 0 ;
169  sdh.observer = "" ;
170  sdh.project = "" ;
171  sdh.obstype = "" ;
172  sdh.antennaname = "" ;
173  sdh.antennaposition.resize( 0 ) ;
174  sdh.equinox = 0.0 ;
175  sdh.freqref = "" ;
176  sdh.reffreq = -1.0 ;
177  sdh.bandwidth = 0.0 ;
178  sdh.utc = 0.0 ;
179  sdh.fluxunit = "" ;
180  sdh.epoch = "" ;
181  sdh.poltype = "" ;
182 
183  // check if optional table exists
184  //const TableRecord msrec = tablesel_.keywordSet() ;
185  const TableRecord msrec = mstable_.keywordSet() ;
186  isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
187  if ( isDoppler_ )
188    if ( mstable_.doppler().nrow() == 0 )
189      isDoppler_ = False ;
190  isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
191  if ( isFlagCmd_ )
192    if ( mstable_.flagCmd().nrow() == 0 )
193      isFlagCmd_ = False ;
194  isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
195  if ( isFreqOffset_ )
196    if ( mstable_.freqOffset().nrow() == 0 )
197      isFreqOffset_ = False ;
198  isHistory_ = msrec.isDefined( "HISTORY" ) ;
199  if ( isHistory_ )
200    if ( mstable_.history().nrow() == 0 )
201      isHistory_ = False ;
202  isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
203  if ( isProcessor_ )
204    if ( mstable_.processor().nrow() == 0 )
205      isProcessor_ = False ;
206  isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
207  if ( isSysCal_ )
208    if ( mstable_.sysCal().nrow() == 0 )
209      isSysCal_ = False ;
210  isWeather_ = msrec.isDefined( "WEATHER" ) ;
211  if ( isWeather_ )
212    if ( mstable_.weather().nrow() == 0 )
213      isWeather_ = False ;
214
215  // Access to MS subtables
216  MSField fieldtab = mstable_.field() ;
217  MSPolarization poltab = mstable_.polarization() ;
218  MSDataDescription ddtab = mstable_.dataDescription() ;
219  MSObservation obstab = mstable_.observation() ;
220  MSSource srctab = mstable_.source() ;
221  MSSpectralWindow spwtab = mstable_.spectralWindow() ;
222  MSSysCal caltab = mstable_.sysCal() ;
223  if ( caltab.nrow() == 0 )
224    isSysCal_ = False ;
225  else {
226    if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
227      colTcal_ = "TCAL" ;
228      if ( !caltab.tableDesc().isColumn( colTcal_ ) )
229        colTcal_ = "NONE" ;
230    }
231    if ( !caltab.tableDesc().isColumn( colTsys_ ) ) {
232      colTsys_ = "TSYS" ;
233      if ( !caltab.tableDesc().isColumn( colTcal_ ) )
234        colTsys_ = "NONE" ;
235    }
236  }
237//   colTcal_ = "TCAL" ;
238//   colTsys_ = "TSYS" ;
239  MSPointing pointtab = mstable_.pointing() ;
240  if ( mstable_.weather().nrow() == 0 )
241    isWeather_ = False ;
242  MSState stattab = mstable_.state() ;
243  MSAntenna anttab = mstable_.antenna() ;
244
245  // TEST
246  // memory allocation by boost::object_pool
247  boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ;
248  //
249
250  // SUBTABLES: FREQUENCIES
251  table_->frequencies().setFrame( "LSRK" ) ;
252  table_->frequencies().setFrame( "LSRK", True ) ;
253
254  // SUBTABLES: WEATHER
255  fillWeather() ;
256
257  // SUBTABLES: FOCUS
258  fillFocus() ;
259
260  // SUBTABLES: TCAL
261  fillTcal( tpoolr ) ;
262
263  // SUBTABLES: FIT
264  //fillFit() ;
265
266  // SUBTABLES: HISTORY
267  //fillHistory() ;
268
269  // shared pointers
270  ROTableColumn *tcolr ;
271
272  // MAIN
273  // Iterate over several ids
274  map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair
275  ROArrayQuantColumn<Double> *sharedQDArrCol = new ROArrayQuantColumn<Double>( anttab, "POSITION" ) ;
276  Vector< Quantum<Double> > antpos = (*sharedQDArrCol)( antenna_ ) ;
277  delete sharedQDArrCol ;
278  MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
279  if ( getPt_ ) {
280    //pointtab = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;
281    pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
282  }
283  tcolr = tpoolr->construct( anttab, "STATION" ) ;
284  String stationName = tcolr->asString( antenna_ ) ;
285  tpoolr->destroy( tcolr ) ;
286  tcolr = tpoolr->construct( anttab, "NAME" ) ;
287  String antennaName = tcolr->asString( antenna_ ) ;
288  tpoolr->destroy( tcolr ) ;
289  sdh.antennaposition.resize( 3 ) ;
290  for ( int i = 0 ; i < 3 ; i++ )
291    sdh.antennaposition[i] = antpos[i].getValue( "m" ) ;
292  String telescopeName = "" ;
293
294//   double time1 = gettimeofday_sec() ;
295//   os_ << "end fill init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
296
297  // row based
298  Table &stab = table_->table() ;
299  TableRow row( stab ) ;
300  TableRecord &trec = row.record() ;
301  RecordFieldPtr< Array<Float> > spRF( trec, "SPECTRA" ) ;
302  RecordFieldPtr< Array<uChar> > ucarrRF( trec, "FLAGTRA" ) ;
303  RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
304  RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
305  RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
306  RecordFieldPtr< Array<Double> > dirRF( trec, "DIRECTION" ) ;
307  RecordFieldPtr<Float> azRF( trec, "AZIMUTH" ) ;
308  RecordFieldPtr<Float> elRF( trec, "ELEVATION" ) ;
309  RecordFieldPtr< Array<Double> > scrRF( trec, "SCANRATE" ) ;
310  RecordFieldPtr<uInt> cycleRF( trec, "CYCLENO" ) ;
311  RecordFieldPtr<uInt> flrRF( trec, "FLAGROW" ) ;
312  RecordFieldPtr<uInt> tcalidRF( trec, "TCAL_ID" ) ;
313  RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ;
314  RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ;
315
316
317  // REFBEAMNO
318  RecordFieldPtr<Int> intRF( trec, "REFBEAMNO" ) ;
319  *intRF = 0 ;
320
321  // FIT_ID
322  intRF.attachToRecord( trec, "FIT_ID" ) ;
323  *intRF = -1 ;
324
325  // OPACITY
326  RecordFieldPtr<Float> floatRF( trec, "OPACITY" ) ;
327  *floatRF = 0.0 ;
328
329  //
330  // ITERATION: OBSERVATION_ID
331  //
332  TableIterator iter0( mstable_, "OBSERVATION_ID" ) ;
333  while( !iter0.pastEnd() ) {
334//     time0 = gettimeofday_sec() ;
335//     os_ << "start 0th iteration: " << time0 << LogIO::POST ;
336    Table t0 = iter0.table() ;
337    tcolr = tpoolr->construct( t0, "OBSERVATION_ID" ) ;
338    Int obsId = tcolr->asInt( 0 ) ;
339    tpoolr->destroy( tcolr ) ;
340    if ( sdh.observer == "" ) {
341      tcolr = tpoolr->construct( obstab, "OBSERVER" ) ;
342      sdh.observer = tcolr->asString( obsId ) ;
343      tpoolr->destroy( tcolr ) ;
344    }
345    if ( sdh.project == "" ) {
346      tcolr = tpoolr->construct( obstab, "PROJECT" ) ;
347      sdh.observer = tcolr->asString( obsId ) ;
348      tpoolr->destroy( tcolr ) ;
349    }
350    ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ;
351    MEpoch me = (*tmpMeasCol)( obsId )( IPosition(1,0) ) ;
352    delete tmpMeasCol ;
353    if ( sdh.utc == 0.0 ) {
354      sdh.utc = me.get( "d" ).getValue() ;
355    }
356    if ( telescopeName == "" ) {
357      tcolr = tpoolr->construct( obstab, "TELESCOPE_NAME" ) ;
358      telescopeName = tcolr->asString( obsId ) ;
359      tpoolr->destroy( tcolr ) ;
360    }
361    Int nbeam = 0 ;
362//     time1 = gettimeofday_sec() ;
363//     os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
364    //
365    // ITERATION: FEED1
366    //
367    TableIterator iter1( t0, "FEED1" ) ;
368    while( !iter1.pastEnd() ) {
369//       time0 = gettimeofday_sec() ;
370//       os_ << "start 1st iteration: " << time0 << LogIO::POST ;
371      Table t1 = iter1.table() ;
372      // assume FEED1 == FEED2
373      tcolr = tpoolr->construct( t1, "FEED1" ) ;
374      Int feedId = tcolr->asInt( 0 ) ;
375      tpoolr->destroy( tcolr ) ;
376      nbeam++ ;
377
378      // BEAMNO
379      RecordFieldPtr<uInt> uintRF( trec, "BEAMNO" ) ;
380      *uintRF = feedId ;
381
382      // FOCUS_ID
383      uintRF.attachToRecord( trec, "FOCUS_ID" ) ;
384      *uintRF = 0 ;
385
386//       time1 = gettimeofday_sec() ;
387//       os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
388      //
389      // ITERATION: FIELD_ID
390      //
391      TableIterator iter2( t1, "FIELD_ID" ) ;
392      while( !iter2.pastEnd() ) {
393//         time0 = gettimeofday_sec() ;
394//         os_ << "start 2nd iteration: " << time0 << LogIO::POST ;
395        Table t2 = iter2.table() ;
396        tcolr = tpoolr->construct( t2, "FIELD_ID" ) ;
397        Int fieldId = tcolr->asInt( 0 ) ;
398        tpoolr->destroy( tcolr ) ;
399        tcolr = tpoolr->construct( fieldtab, "SOURCE_ID" ) ;
400        Int srcId = tcolr->asInt( fieldId ) ;
401        tpoolr->destroy( tcolr ) ;
402        tcolr = tpoolr->construct( fieldtab, "NAME" ) ;
403        String fieldName = tcolr->asString( fieldId ) + "__" + String::toString(fieldId) ;
404        tpoolr->destroy( tcolr ) ;
405        ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ;
406        Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ;
407        delete delayDirCol ;
408        Vector<Double> defaultScanrate( 2, 0.0 ) ;
409        Vector<Double> defaultDir = delayDir[0].getAngle( "rad" ).getValue() ;
410        if ( delayDir.nelements() > 1 )
411          defaultScanrate = delayDir[1].getAngle( "rad" ).getValue() ;
412         
413
414        // FIELDNAME
415        RecordFieldPtr<String> strRF( trec, "FIELDNAME" ) ;
416        *strRF = fieldName ;
417
418
419//         time1 = gettimeofday_sec() ;
420//         os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
421        //
422        // ITERATION: DATA_DESC_ID
423        //
424        TableIterator iter3( t2, "DATA_DESC_ID" ) ;
425        while( !iter3.pastEnd() ) {
426//           time0 = gettimeofday_sec() ;
427//           os_ << "start 3rd iteration: " << time0 << LogIO::POST ;
428          Table t3 = iter3.table() ;
429          tcolr = tpoolr->construct( t3, "DATA_DESC_ID" ) ;
430          Int ddId = tcolr->asInt( 0 ) ;
431          tpoolr->destroy( tcolr ) ;
432          tcolr = tpoolr->construct( ddtab, "POLARIZATION_ID" ) ;
433          Int polId = tcolr->asInt( ddId ) ;
434          tpoolr->destroy( tcolr ) ;
435          tcolr = tpoolr->construct( ddtab, "SPECTRAL_WINDOW_ID" ) ;
436          Int spwId = tcolr->asInt( ddId ) ;
437          tpoolr->destroy( tcolr ) ;
438
439          // IFNO
440          uintRF.attachToRecord( trec, "IFNO" ) ;
441          *uintRF = (uInt)spwId ;
442
443          // polarization information
444          tcolr = tpoolr->construct( poltab, "NUM_CORR" ) ;
445          Int npol = tcolr->asInt( polId ) ;
446          tpoolr->destroy( tcolr ) ;
447          ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ;
448          Vector<Int> corrtype = (*roArrICol)( polId ) ;
449          delete roArrICol ;
450//           os_ << "npol = " << npol << LogIO::POST ;
451//           os_ << "corrtype = " << corrtype << LogIO::POST ;
452          // source information
453//           os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;
454          MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;
455          if ( srctabSel.nrow() == 0 ) {
456            srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;
457          }
458          String srcName( "" ) ;
459          Vector<Double> srcPM( 2, 0.0 ) ;
460          Vector<Double> srcDir( 2, 0.0 ) ;
461          MDirection md ;
462          Int numLines = 0 ;
463          ROArrayColumn<Double> *roArrDCol = 0 ;
464          if ( srctabSel.nrow() > 0 ) {
465            // source name
466            tcolr = tpoolr->construct( srctabSel, "NAME" ) ;
467            srcName = tcolr->asString( 0 ) ;
468            tpoolr->destroy( tcolr ) ;
469
470            // source proper motion
471            roArrDCol = new ROArrayColumn<Double>( srctabSel, "PROPER_MOTION" ) ;
472            srcPM = (*roArrDCol)( 0 ) ;
473            delete roArrDCol ;
474           
475            // source direction
476            roArrDCol = new ROArrayColumn<Double>( srctabSel, "DIRECTION" ) ;
477            srcDir = (*roArrDCol)( 0 ) ;
478            delete roArrDCol ;
479
480            // source direction as MDirection object
481            ROScalarMeasColumn<MDirection> *tmpMeasCol = new ROScalarMeasColumn<MDirection>( srctabSel, "DIRECTION" ) ;
482            md = (*tmpMeasCol)( 0 ) ;
483            delete tmpMeasCol ;
484
485            // number of lines
486            tcolr = tpoolr->construct( srctabSel, "NUM_LINES" ) ;
487            numLines = tcolr->asInt( 0 ) ;
488            tpoolr->destroy( tcolr ) ;
489
490          }
491          else {
492            md = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
493          }
494
495          // SRCNAME
496          strRF.attachToRecord( trec, "SRCNAME" ) ;
497          *strRF = srcName ;
498
499//           os_ << "srcName = " << srcName << LogIO::POST ;
500
501          // SRCPROPERMOTION
502          RecordFieldPtr< Array<Double> > darrRF( trec, "SRCPROPERMOTION" ) ;
503          *darrRF = srcPM ;
504
505          //os_ << "srcPM = " << srcPM << LogIO::POST ;
506
507          // SRCDIRECTION
508          darrRF.attachToRecord( trec, "SRCDIRECTION" ) ;
509          *darrRF = srcDir ;
510
511          //os_ << "srcDir = " << srcDir << LogIO::POST ;
512
513          // for MOLECULES subtable
514//           os_ << "numLines = " << numLines << LogIO::POST ;
515
516          Vector<Double> restFreqs( numLines, 0.0 ) ;
517          Vector<String> transitionName( numLines, "" ) ;
518          Vector<Double> sysVels ;
519          Double sysVel = 0.0 ;
520          if ( numLines != 0 ) {
521            if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) {
522              sharedQDArrCol = new ROArrayQuantColumn<Double>( srctabSel, "REST_FREQUENCY" ) ;
523              Array< Quantum<Double> > qRestFreqs = (*sharedQDArrCol)( 0 ) ;
524              delete sharedQDArrCol ;
525              for ( int i = 0 ; i < numLines ; i++ ) {
526                restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
527              }
528            }
529//             os_ << "restFreqs = " << restFreqs << LogIO::POST ;
530            if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
531              ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ;
532              if ( transitionCol.isDefined( 0 ) )
533                transitionName = transitionCol( 0 ) ;
534              //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ;
535            }
536            if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) {
537              roArrDCol = new ROArrayColumn<Double>( srctabSel, "SYSVEL" ) ;
538              sysVels = (*roArrDCol)( 0 ) ;
539              delete roArrDCol ;
540            }
541            if ( !sysVels.empty() ) {
542              //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ;
543              // NB: assume all SYSVEL values are the same
544              sysVel = sysVels( IPosition(1,0) ) ;
545            }
546          }
547
548          // SRCVELOCITY
549          RecordFieldPtr<Double> doubleRF( trec, "SRCVELOCITY" ) ;
550          *doubleRF = sysVel ;
551
552//           os_ << "sysVel = " << sysVel << LogIO::POST ;
553
554          uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ;
555
556          // MOLECULE_ID
557          uintRF.attachToRecord( trec, "MOLECULE_ID" ) ;
558          *uintRF = molId ;
559
560          // spectral setup
561          MeasFrame mf( me, mp, md ) ;
562          tcolr = tpoolr->construct( spwtab, "MEAS_FREQ_REF" ) ;
563          MFrequency::Types freqRef = MFrequency::castType((uInt)(tcolr->asInt(spwId))) ;
564          tpoolr->destroy( tcolr ) ;
565          tcolr = tpoolr->construct( spwtab, "NUM_CHAN" ) ;
566          Int nchan = tcolr->asInt( spwId ) ;
567          Bool iswvr = False ;
568          if ( nchan == 4 ) iswvr = True ;
569          tpoolr->destroy( tcolr ) ;
570          Bool even = False ;
571          if ( (nchan/2)*2 == nchan ) even = True ;
572          sdh.nchan = max( sdh.nchan, nchan ) ;
573          ROScalarQuantColumn<Double> *tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "TOTAL_BANDWIDTH" ) ;
574          Double totbw = (*tmpQuantCol)( spwId ).getValue( "Hz" ) ;
575          delete tmpQuantCol ;
576          sdh.bandwidth = max( sdh.bandwidth, totbw ) ;
577          if ( sdh.freqref == "" )
578            //sdh.freqref = MFrequency::showType( freqRef ) ;
579            sdh.freqref = "LSRK" ;
580          if ( sdh.reffreq == -1.0 ) {
581            tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "REF_FREQUENCY" ) ;
582            Quantum<Double> qreffreq = (*tmpQuantCol)( spwId ) ;
583            delete tmpQuantCol ;
584            if ( freqRef == MFrequency::LSRK ) {
585              sdh.reffreq = qreffreq.getValue("Hz") ;
586            }
587            else {
588              MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
589              sdh.reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
590            }
591          }
592          Int refchan = nchan / 2 ;
593          IPosition refip( 1, refchan ) ;
594          Double refpix = 0.5*(nchan-1) ;
595          Double refval = 0.0 ;
596          sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_WIDTH" ) ;
597          Double increment = (*sharedQDArrCol)( spwId )( refip ).getValue( "Hz" ) ;
598          delete sharedQDArrCol ;
599//           os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
600          sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_FREQ" ) ;
601          Vector< Quantum<Double> > chanFreqs = (*sharedQDArrCol)( spwId ) ;
602          delete sharedQDArrCol ;
603          if ( freqRef == MFrequency::LSRK ) {
604            if ( even ) {
605              IPosition refip0( 1, refchan-1 ) ;
606              Double refval0 = chanFreqs(refip0).getValue("Hz") ;
607              Double refval1 = chanFreqs(refip).getValue("Hz") ;
608              refval = 0.5 * ( refval0 + refval1 ) ;
609            }
610            else {
611              refval = chanFreqs(refip).getValue("Hz") ;
612            }
613          }
614          else {
615            MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
616            if ( even ) {
617              IPosition refip0( 1, refchan-1 ) ;
618              Double refval0 = chanFreqs(refip0).getValue("Hz") ;
619              Double refval1 = chanFreqs(refip).getValue("Hz") ;
620              refval = 0.5 * ( refval0 + refval1 ) ;
621              refval = tolsr( refval ).get("Hz").getValue() ;
622            }
623            else {
624              refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
625            }
626          }
627          uInt freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
628          if ( ifmap.find( spwId ) == ifmap.end() ) {
629            ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
630            //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
631          }
632
633          // FREQ_ID
634          uintRF.attachToRecord( trec, "FREQ_ID" ) ;
635          *uintRF = freqId ;
636
637          // for TSYS and TCAL
638          Vector<MEpoch> scTime ;
639          Vector<Double> scInterval ;
640          ROArrayColumn<Float> scTsysCol ;
641          MSSysCal caltabsel ;
642          if ( isSysCal_ ) {
643            caltabsel = caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ;
644            ROScalarMeasColumn<MEpoch> scTimeCol( caltabsel, "TIME" ) ;
645            scTime.resize( caltabsel.nrow() ) ;
646            for ( uInt irow = 0 ; irow < caltabsel.nrow() ; irow++ )
647              scTime[irow] = scTimeCol( irow ) ;
648            ROScalarColumn<Double> *scIntervalCol = new ROScalarColumn<Double>( caltabsel, "INTERVAL" ) ;
649            scInterval = scIntervalCol->getColumn() ;
650            delete scIntervalCol ;
651            if ( colTsys_ != "NONE" )
652              scTsysCol.attach( caltabsel, colTsys_ ) ;
653          }
654
655          sdh.npol = max( sdh.npol, npol ) ;
656          if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
657
658//           time1 = gettimeofday_sec() ;
659//           os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
660          //
661          // ITERATION: SCAN_NUMBER
662          //
663          TableIterator iter4( t3, "SCAN_NUMBER" ) ;
664          while( !iter4.pastEnd() ) {
665//             time0 = gettimeofday_sec() ;
666//             os_ << "start 4th iteration: " << time0 << LogIO::POST ;
667            Table t4 = iter4.table() ;
668            tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ;
669            Int scanNum = tcolr->asInt( 0 ) ;
670            tpoolr->destroy( tcolr ) ;
671
672            // SCANNO
673            uintRF.attachToRecord( trec, "SCANNO" ) ;
674            *uintRF = scanNum - 1 ;
675
676//             time1 = gettimeofday_sec() ;
677//             os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
678            //
679            // ITERATION: STATE_ID
680            //
681            TableIterator iter5( t4, "STATE_ID" ) ;
682            while( !iter5.pastEnd() ) {
683//               time0 = gettimeofday_sec() ;
684//               os_ << "start 5th iteration: " << time0 << LogIO::POST ;
685              Table t5 = iter5.table() ;
686              tcolr = tpoolr->construct( t5, "STATE_ID" ) ;
687              Int stateId = tcolr->asInt( 0 ) ;
688              tpoolr->destroy( tcolr ) ;
689              tcolr = tpoolr->construct( stattab, "OBS_MODE" ) ;
690              String obstype = tcolr->asString( stateId ) ;
691              tpoolr->destroy( tcolr ) ;
692              if ( sdh.obstype == "" ) sdh.obstype = obstype ;
693
694              Int nrow = t5.nrow() ;
695//               time1 = gettimeofday_sec() ;
696//               os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
697
698              uInt cycle = 0 ;
699
700              Cube<Float> spArr ;
701              Cube<Bool> flArr ;
702              if ( isFloatData_ ) {
703                ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
704                ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ;
705                spArr = mFloatDataCol.getColumn() ;
706                flArr = mFlagCol.getColumn() ;
707                if ( sdh.fluxunit == "" ) {
708                  const TableRecord &dataColKeys = mFloatDataCol.keywordSet() ;
709                  if ( dataColKeys.isDefined( "UNIT" ) )
710                    sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
711                }
712              }
713              else if ( isData_ ) {
714                spArr.resize( npol, nchan, nrow ) ;
715                flArr.resize( npol, nchan, nrow ) ;
716                ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
717                ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ;
718                for ( Int irow = 0 ; irow < nrow ; irow++ ) {
719                  Bool crossOK = False ;
720                  Matrix<Complex> sp = mDataCol( irow ) ;
721                  Matrix<Bool> fl = mFlagCol( irow ) ;
722                  Matrix<Float> spxy = spArr.xyPlane( irow ) ;
723                  Matrix<Bool> flxy = flArr.xyPlane( irow ) ;
724                  for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
725                    if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
726                         || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
727                      if ( !crossOK ) {
728                        spxy.row( ipol ) = real( sp.row( ipol ) ) ;
729                        flxy.row( ipol ) = fl.row( ipol ) ;
730                        if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) {
731                          spxy.row( ipol+1 ) = imag( sp.row( ipol ) ) ;
732                          flxy.row( ipol+1 ) = fl.row( ipol ) ;
733                        }                       
734                        else {
735                          spxy.row( ipol+1 ) = imag( conj( sp.row( ipol ) ) ) ;
736                          flxy.row( ipol+1 ) = fl.row( ipol ) ;
737                        }
738                        crossOK = True ;
739                      }
740                    }
741                    else {
742                      spxy.row( ipol ) = real( sp.row( ipol ) ) ;
743                      flxy.row( ipol ) = fl.row( ipol ) ;
744                    }
745                  }
746                }
747                if ( sdh.fluxunit == "" ) {
748                  const TableRecord &dataColKeys = mDataCol.keywordSet() ;
749                  if ( dataColKeys.isDefined( "UNIT" ) )
750                    sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
751                }
752              }
753              ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
754              Block<MEpoch> mTimeB( nrow ) ;
755              for ( Int irow = 0 ; irow < nrow ; irow++ )
756                mTimeB[irow] = (*mTimeCol)( irow ) ;
757              ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;
758              ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;
759              Block<Int> sysCalIdx( nrow, -1 ) ;
760              if ( isSysCal_ ) {
761                getSysCalTime( scTime, scInterval, mTimeB, sysCalIdx ) ;
762              }
763              delete mTimeCol ;
764              Matrix<Float> defaulttsys( npol, 1, 1.0 ) ;
765              Int srcType = getSrcType( stateId, tpoolr ) ;
766              uInt diridx = 0 ;
767              MDirection::Types dirType ;
768              uInt wid = 0 ;
769              Int pidx = 0 ;
770              Bool crossOK = False ;
771              Block<uInt> polnos( npol, 99 ) ;
772              for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
773                Block<uInt> p = getPolNo( corrtype[ipol] ) ;
774                if ( p.size() > 1 ) {
775                  if ( crossOK ) continue ;
776                  else {
777                    polnos[pidx] = p[0] ;
778                    pidx++ ;
779                    polnos[pidx] = p[1] ;
780                    pidx++ ;
781                    crossOK = True ;
782                  }
783                }
784                else {
785                  polnos[pidx] = p[0] ;
786                  pidx++ ;
787                }
788              }
789             
790              // SRCTYPE
791              intRF.attachToRecord( trec, "SRCTYPE" ) ;
792              *intRF = srcType ;
793
794              for ( Int irow = 0 ; irow < nrow ; irow++ ) {
795                // CYCLENO
796                *cycleRF = cycle ;
797
798                // FLAGROW
799                *flrRF = (uInt)mFlagRowCol->asBool( irow ) ;
800
801                // SPECTRA, FLAG
802                Matrix<Float> sp = spArr.xyPlane( irow ) ;
803                Matrix<Bool> flb = flArr.xyPlane( irow ) ;
804                Matrix<uChar> fl( flb.shape() ) ;
805                convertArray( fl, flb ) ;
806
807                // TIME
808                *timeRF = mTimeB[irow].get("d").getValue() ;
809
810                // INTERVAL
811                *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ;
812
813                // TSYS
814                Matrix<Float> tsys ;
815                if ( sysCalIdx[irow] != -1 && colTsys_ != "NONE" )
816                  tsys = scTsysCol( irow ) ;
817                else
818                  tsys = defaulttsys ;
819
820                // TCAL_ID
821                Block<uInt> tcalids( npol, 0 ) ;
822                if ( sysCalIdx[irow] != -1 && colTcal_ != "NONE" ) {
823                  tcalids = getTcalId( feedId, spwId, scTime[sysCalIdx[irow]] ) ;
824                }
825
826                // WEATHER_ID
827                if ( isWeather_ )
828                  wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ;
829                *widRF = wid ;
830                 
831
832                // DIRECTION, AZEL, SCANRATE
833                if ( getPt_ ) {
834                  Vector<Double> dir ;
835                  Vector<Double> scanrate ;
836                  String refString ;
837                  diridx = getDirection( diridx, dir, scanrate, refString, pointtab, mTimeB[irow].get("s").getValue() ) ;
838                  MDirection::getType( dirType, refString ) ;
839                  mf.resetEpoch( mTimeB[irow] ) ;
840                  mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ;
841                  if ( refString == "J2000" ) {
842                    *dirRF = dir ;
843                    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
844                    Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
845                    *azRF = (Float)azel(0) ;
846                    *elRF = (Float)azel(1) ;
847                  }
848                  else if ( refString(0,4) == "AZEL" ) {
849                    *azRF = (Float)dir(0) ;
850                    *elRF = (Float)dir(1) ;
851                    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
852                    Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
853                    *dirRF = newdir ;
854                  }
855                  else {
856                    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
857                    Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
858                    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
859                    Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
860                    *dirRF = newdir ;
861                    *azRF = (Float)azel(0) ;
862                    *elRF = (Float)azel(1) ;
863                  }
864                  if ( scanrate.size() != 0 ) {
865                    *scrRF = scanrate ;
866                  }
867                  else {
868                    *scrRF = defaultScanrate ;
869                  }
870                }
871                else {
872                  String ref = md.getRefString() ;
873                  //Vector<Double> defaultDir = srcDir ;
874                  MDirection::getType( dirType, "J2000" ) ;
875                  if ( ref != "J2000" ) {
876                    ROScalarMeasColumn<MEpoch> tmCol( pointtab, "TIME" ) ;
877                    mf.resetEpoch( tmCol( 0 ) ) ;
878                    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
879                    defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ;
880                  }
881                  mf.resetEpoch( mTimeB[irow] ) ;
882                  MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
883                  Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ;
884                  *azRF = (Float)azel(0) ;
885                  *elRF = (Float)azel(1) ;
886                  *dirRF = defaultDir ;
887                  *scrRF = defaultScanrate ;
888                }
889
890                // Polarization dependent things
891                for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
892                  // POLNO
893                  *polnoRF = polnos[ipol] ;
894
895                  //*spRF = sp.row( ipol ) ;
896                  //*ucarrRF = fl.row( ipol ) ;
897                  //*tsysRF = tsys.row( ipol ) ;
898                  spRF.define( sp.row( ipol ) ) ;
899                  ucarrRF.define( fl.row( ipol ) ) ;
900                  tsysRF.define( tsys.row( ipol ) ) ;
901                  *tcalidRF = tcalids[ipol] ;
902
903                  // Commit row
904                  stab.addRow() ;
905                  row.put( stab.nrow()-1 ) ;
906                }
907
908                cycle++ ;
909              }
910              tpoolr->destroy( mIntervalCol ) ;
911              tpoolr->destroy( mFlagRowCol ) ;
912
913//               time1 = gettimeofday_sec() ;
914//               os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
915
916              iter5.next() ;
917            }
918
919
920            iter4.next() ;
921          }
922
923
924          iter3.next() ;
925        }
926
927             
928        iter2.next() ;
929      }
930
931
932      iter1.next() ;
933    }
934    if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
935
936    iter0.next() ;
937  }
938
939
940  delete tpoolr ;
941
942
943  // Table Keywords
944  sdh.nif = ifmap.size() ;
945  if ( ( telescopeName == "" ) || ( antennaName == telescopeName ) ) {
946    sdh.antennaname = antennaName ;
947  }
948  else {
949    sdh.antennaname = telescopeName + "//" + antennaName ;
950  }
951  if ( stationName != "" ) {
952    sdh.antennaname += "@" + stationName ;
953  }
954  ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ;
955  String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ;
956  if ( dirref == "AZELGEO" || dirref == "AZEL" ) {
957    dirref = "J2000" ;
958  }
959  sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
960  sdh.epoch = "UTC" ;
961  if (sdh.freqref == "TOPO") {
962    sdh.freqref = "TOPOCENT";
963  } else if (sdh.freqref == "GEO") {
964    sdh.freqref = "GEOCENTR";
965  } else if (sdh.freqref == "BARY") {
966    sdh.freqref = "BARYCENT";
967  } else if (sdh.freqref == "GALACTO") {
968    sdh.freqref = "GALACTOC";
969  } else if (sdh.freqref == "LGROUP") {
970    sdh.freqref = "LOCALGRP";
971  } else if (sdh.freqref == "CMB") {
972    sdh.freqref = "CMBDIPOL";
973  } else if (sdh.freqref == "REST") {
974    sdh.freqref = "SOURCE";
975  }
976  table_->setHeader( sdh ) ;
977
978  // save path to POINTING table
979  // 2011/3/2 TN
980  // So far, path to original POINTING table is always stored
981  // since sd tasks and regressions don't support getpt control
982  //if ( !getPt_ ) {
983  Path datapath( tablename_ ) ;
984  String pTabName = datapath.absoluteName() + "/POINTING" ;
985  stab.rwKeywordSet().define( "POINTING", pTabName ) ;
986  //}
987
988  // for GBT
989  if ( antennaName.contains( "GBT" ) ) {
990    String goTabName = datapath.absoluteName() + "/GBT_GO" ;
991    stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
992  }
993
994  // for MS created from ASDM
995  //mstable_.keywordSet().print(cout) ;
996  const TableRecord &msKeys = mstable_.keywordSet() ;
997  uInt nfields = msKeys.nfields() ;
998  for ( uInt ifield = 0 ; ifield < nfields ; ifield++ ) {
999    String name = msKeys.name( ifield ) ;
1000    //os_ << "name = " << name << LogIO::POST ;
1001    if ( name.find( "ASDM" ) != String::npos ) {
1002      String asdmpath = msKeys.asTable( ifield ).tableName() ;
1003      os_ << "ASDM table: " << asdmpath << LogIO::POST ;
1004      stab.rwKeywordSet().define( name, asdmpath ) ;
1005    }
1006  }
1007
1008//   double endSec = gettimeofday_sec() ;
1009//   os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1010}
1011
1012void MSFiller::close()
1013{
1014  //tablesel_.closeSubTables() ;
1015  mstable_.closeSubTables() ;
1016  //tablesel_.unlock() ;
1017  mstable_.unlock() ;
1018}
1019
1020Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
1021{
1022//   double startSec = gettimeofday_sec() ;
1023//   os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
1024
1025  MSState statetab = mstable_.state() ;
1026  ROTableColumn *sharedCol ;
1027  sharedCol = tpool->construct( statetab, "OBS_MODE" ) ;
1028  String obsMode = sharedCol->asString( stateId ) ;
1029  tpool->destroy( sharedCol ) ;
1030  sharedCol = tpool->construct( statetab, "SIG" ) ;
1031  Bool sig = sharedCol->asBool( stateId ) ;
1032  tpool->destroy( sharedCol ) ;
1033  sharedCol = tpool->construct( statetab, "REF" ) ;
1034  Bool ref = sharedCol->asBool( stateId ) ;
1035  tpool->destroy( sharedCol ) ;
1036  sharedCol = tpool->construct( statetab, "CAL" ) ;
1037  Double cal = (Double)(sharedCol->asdouble( stateId )) ;
1038  tpool->destroy( sharedCol ) ;
1039  //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
1040
1041  // determine separator
1042  String sep = "" ;
1043  String tmpStr = obsMode.substr( 0, obsMode.find_first_of( "," ) ) ;
1044  //os_ << "tmpStr = " << tmpStr << LogIO::POST ;
1045  //if ( obsMode.find( ":" ) != String::npos ) {
1046  if ( tmpStr.find( ":" ) != String::npos ) {
1047    sep = ":" ;
1048  }
1049  //else if ( obsMode.find( "." ) != String::npos ) {
1050  else if ( tmpStr.find( "." ) != String::npos ) {
1051    sep = "." ;
1052  }
1053  //else if ( obsMode.find( "_" ) != String::npos ) {
1054  else if ( tmpStr.find( "_" ) != String::npos ) {
1055    sep = "_" ;
1056  }
1057  //os_ << "separator = " << sep << LogIO::POST ;
1058
1059  // determine SRCTYPE
1060  Int srcType = SrcType::NOTYPE ;
1061  if ( sep == ":" ) {
1062    // sep == ":"
1063    //
1064    // GBT case
1065    //
1066    // obsMode1=Nod
1067    //    NOD
1068    // obsMode1=OffOn
1069    //    obsMode2=PSWITCHON:  PSON
1070    //    obsMode2=PSWITCHOFF: PSOFF
1071    // obsMode1=??
1072    //    obsMode2=FSWITCH:
1073    //       SIG=1: FSON
1074    //       REF=1: FSOFF
1075    // Calibration scan if CAL != 0
1076    Int epos = obsMode.find_first_of( sep ) ;
1077    Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
1078    String obsMode1 = obsMode.substr( 0, epos ) ;
1079    String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
1080    if ( obsMode1 == "Nod" ) {
1081      srcType = SrcType::NOD ;
1082    }
1083    else if ( obsMode1 == "OffOn" ) {
1084      if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
1085      if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
1086    }
1087    else {
1088      if ( obsMode2 == "FSWITCH" ) {
1089        if ( sig ) srcType = SrcType::FSON ;
1090        if ( ref ) srcType = SrcType::FSOFF ;
1091      }
1092    }
1093    if ( cal > 0.0 ) {
1094      if ( srcType == SrcType::NOD )
1095        srcType = SrcType::NODCAL ;
1096      else if ( srcType == SrcType::PSON )
1097        srcType = SrcType::PONCAL ;
1098      else if ( srcType == SrcType::PSOFF )
1099        srcType = SrcType::POFFCAL ;
1100      else if ( srcType == SrcType::FSON )
1101        srcType = SrcType::FONCAL ;
1102      else if ( srcType == SrcType::FSOFF )
1103        srcType = SrcType::FOFFCAL ;
1104      else
1105        srcType = SrcType::CAL ;
1106    }
1107  }
1108  else if ( sep == "." ) {
1109    // sep == "."
1110    //
1111    // ALMA & EVLA case (MS via ASDM) before3.1
1112    //
1113    // obsMode1=CALIBRATE_*
1114    //    obsMode2=ON_SOURCE: PONCAL
1115    //    obsMode2=OFF_SOURCE: POFFCAL
1116    // obsMode1=OBSERVE_TARGET
1117    //    obsMode2=ON_SOURCE: PON
1118    //    obsMode2=OFF_SOURCE: POFF
1119    string substr[2] ;
1120    int numSubstr = split( obsMode, substr, 2, "," ) ;
1121    //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
1122    //for ( int i = 0 ; i < numSubstr ; i++ )
1123    //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
1124    String obsType( substr[0] ) ;
1125    //os_ << "obsType = " << obsType << LogIO::POST ;
1126    Int epos = obsType.find_first_of( sep ) ;
1127    Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
1128    String obsMode1 = obsType.substr( 0, epos ) ;
1129    String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
1130    //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
1131    //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
1132    if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
1133      if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
1134      if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
1135    }
1136    else if ( obsMode1 == "OBSERVE_TARGET" ) {
1137      if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
1138      if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
1139    }
1140  }
1141  else if ( sep == "_" ) {
1142    // sep == "_"
1143    //
1144    // ALMA & EVLA case (MS via ASDM) after 3.2
1145    //
1146    // obsMode1=CALIBRATE_*
1147    //    obsMode2=ON_SOURCE: PONCAL
1148    //    obsMode2=OFF_SOURCE: POFFCAL
1149    // obsMode1=OBSERVE_TARGET
1150    //    obsMode2=ON_SOURCE: PON
1151    //    obsMode2=OFF_SOURCE: POFF
1152    string substr[2] ;
1153    int numSubstr = split( obsMode, substr, 2, "," ) ;
1154    //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
1155    //for ( int i = 0 ; i < numSubstr ; i++ )
1156    //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
1157    String obsType( substr[0] ) ;
1158    //os_ << "obsType = " << obsType << LogIO::POST ;
1159    string substr2[4] ;
1160    int numSubstr2 = split( obsType, substr2, 4, sep ) ;
1161    //Int epos = obsType.find_first_of( sep ) ;
1162    //Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
1163    //String obsMode1 = obsType.substr( 0, epos ) ;
1164    //String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
1165    String obsMode1( substr2[0] ) ;
1166    String obsMode2( substr2[2] ) ;
1167    //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
1168    //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
1169    if ( obsMode1.find( "CALIBRATE" ) == 0 ) {
1170      if ( obsMode2 == "ON" ) srcType = SrcType::PONCAL ;
1171      if ( obsMode2 == "OFF" ) srcType = SrcType::POFFCAL ;
1172    }
1173    else if ( obsMode1 == "OBSERVE" ) {
1174      if ( obsMode2 == "ON" ) srcType = SrcType::PSON ;
1175      if ( obsMode2 == "OFF" ) srcType = SrcType::PSOFF ;
1176    }
1177  }
1178  else {
1179    if ( sig ) srcType = SrcType::SIG ;
1180    if ( ref ) srcType = SrcType::REF ;
1181  }
1182   
1183  //os_ << "srcType = " << srcType << LogIO::POST ;
1184//   double endSec = gettimeofday_sec() ;
1185//   os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1186  return srcType ;
1187}
1188
1189//Vector<uInt> MSFiller::getPolNo( Int corrType )
1190Block<uInt> MSFiller::getPolNo( Int corrType )
1191{
1192//   double startSec = gettimeofday_sec() ;
1193//   os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
1194  Block<uInt> polno( 1 ) ;
1195
1196  if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {
1197    polno = 0 ;
1198  }
1199  else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {
1200    polno = 1 ;
1201  }
1202  else if ( corrType == Stokes::U ) {
1203    polno = 2 ;
1204  }
1205  else if ( corrType == Stokes::V ) {
1206    polno = 3 ;
1207  }
1208  else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {
1209    polno.resize( 2 ) ;
1210    polno[0] = 2 ;
1211    polno[1] = 3 ;
1212  }
1213  else if ( corrType == Stokes::Plinear ) {
1214    polno[0] = 1 ;
1215  }
1216  else if ( corrType == Stokes::Pangle ) {
1217    polno[0] = 2 ;
1218  }
1219  else {
1220    polno = 99 ;
1221  }
1222  //os_ << "polno = " << polno << LogIO::POST ;
1223//   double endSec = gettimeofday_sec() ;
1224//   os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1225 
1226  return polno ;
1227}
1228
1229String MSFiller::getPolType( Int corrType )
1230{
1231//   double startSec = gettimeofday_sec() ;
1232//   os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
1233  String poltype = "" ;
1234
1235  if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )
1236    poltype = "stokes" ;
1237  else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )
1238    poltype = "linear" ;
1239  else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )
1240    poltype = "circular" ;
1241  else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )
1242    poltype = "linpol" ;
1243
1244//   double endSec = gettimeofday_sec() ;
1245//   os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1246  return poltype ;
1247}
1248
1249void MSFiller::fillWeather()
1250{
1251//   double startSec = gettimeofday_sec() ;
1252//   os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1253
1254  if ( !isWeather_ ) {
1255    // add dummy row
1256    table_->weather().table().addRow(1,True) ;
1257    return ;
1258  }
1259
1260  Table mWeather = mstable_.weather()  ;
1261  //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1262  Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1263  //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1264  if ( mWeatherSel.nrow() == 0 ) {
1265    os_ << "No rows with ANTENNA_ID = " << antenna_ << " in WEATHER table, Try -1..." << LogIO::POST ;
1266    mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1267    if ( mWeatherSel.nrow() == 0 ) {
1268      os_ << "No rows in WEATHER table" << LogIO::POST ;
1269    }
1270  }
1271  uInt wnrow = mWeatherSel.nrow() ;
1272  //os_ << "wnrow = " << wnrow << LogIO::POST ;
1273
1274  if ( wnrow == 0 )
1275    return ;
1276
1277  Table wtab = table_->weather().table() ;
1278  wtab.addRow( wnrow ) ;
1279
1280  ScalarColumn<Float> *fCol ;
1281  ROScalarColumn<Float> *sharedFloatCol ;
1282  if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1283    fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1284    sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1285    fCol->putColumn( *sharedFloatCol ) ;
1286    delete sharedFloatCol ;
1287    delete fCol ;
1288  }
1289  if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1290    fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1291    sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1292    fCol->putColumn( *sharedFloatCol ) ;
1293    delete sharedFloatCol ;
1294    delete fCol ;
1295  }
1296  if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1297    fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1298    sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1299    fCol->putColumn( *sharedFloatCol ) ;
1300    delete sharedFloatCol ;
1301    delete fCol ;
1302  }
1303  if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) { 
1304    fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1305    sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1306    fCol->putColumn( *sharedFloatCol ) ;
1307    delete sharedFloatCol ;
1308    delete fCol ;
1309  }
1310  if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1311    fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1312    sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1313    fCol->putColumn( *sharedFloatCol ) ;
1314    delete sharedFloatCol ;
1315    delete fCol ;
1316  }
1317  ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1318  for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1319    idCol.put( irow, irow ) ;
1320
1321  ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1322  ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1323  String tUnit = tqCol.getUnits() ;
1324  mwTime_ = tCol.getColumn() ;
1325  if ( tUnit == "d" )
1326    mwTime_ *= 86400.0 ;
1327  tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1328  tCol.attach( mWeatherSel, "INTERVAL" ) ;
1329  String iUnit = tqCol.getUnits() ;
1330  mwInterval_ = tCol.getColumn() ;
1331  if ( iUnit == "d" )
1332    mwInterval_ *= 86400.0 ;
1333  //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1334//   double endSec = gettimeofday_sec() ;
1335//   os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1336}
1337
1338void MSFiller::fillFocus()
1339{
1340//   double startSec = gettimeofday_sec() ;
1341//   os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
1342  // tentative
1343  Table tab = table_->focus().table() ;
1344  tab.addRow( 1 ) ;
1345  ScalarColumn<uInt> idCol( tab, "ID" ) ;
1346  idCol.put( 0, 0 ) ;
1347//   double endSec = gettimeofday_sec() ;
1348//   os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1349}
1350
1351void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr )
1352{
1353//   double startSec = gettimeofday_sec() ;
1354//   os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
1355
1356  if ( !isSysCal_ ) {
1357    // add dummy row
1358    os_ << "No SYSCAL rows" << LogIO::POST ;
1359    table_->tcal().table().addRow(1,True) ;
1360    Vector<Float> defaultTcal( 1, 1.0 ) ;
1361    ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1362    tcalCol.put( 0, defaultTcal ) ;
1363    return ;
1364  }
1365
1366  if ( colTcal_ == "NONE" ) {
1367    // add dummy row
1368    os_ << "No TCAL column" << LogIO::POST ;
1369    table_->tcal().table().addRow(1,True) ;
1370    Vector<Float> defaultTcal( 1, 1.0 ) ;
1371    ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1372    tcalCol.put( 0, defaultTcal ) ;
1373    return ;
1374  }
1375
1376  Table sctab = mstable_.sysCal() ;
1377  if ( sctab.nrow() == 0 ) {
1378    os_ << "No SYSCAL rows" << LogIO::POST ;
1379    return ;
1380  }
1381  Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
1382  if ( sctabsel.nrow() == 0 ) {
1383    os_ << "No SYSCAL rows" << LogIO::POST ;
1384    return ;
1385  }
1386  ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, colTcal_ ) ;
1387  // return if any rows without Tcal value exists
1388  Bool notDefined = False ;
1389  for ( uInt irow = 0 ; irow < sctabsel.nrow() ; irow++ ) {
1390    if ( !tmpTcalCol->isDefined( irow ) ) {
1391      notDefined = True ;
1392      break ;
1393    }
1394  }
1395  if ( notDefined ) {
1396    os_ << "No TCAL value" << LogIO::POST ;
1397    delete tmpTcalCol ;
1398    table_->tcal().table().addRow(1,True) ;
1399    Vector<Float> defaultTcal( 1, 1.0 ) ;
1400    ArrayColumn<Float> tcalCol( table_->tcal().table(), "TCAL" ) ;
1401    tcalCol.put( 0, defaultTcal ) ;
1402    return ;
1403  }   
1404  uInt npol = tmpTcalCol->shape( 0 )(0) ;
1405  delete tmpTcalCol ;
1406  //os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
1407  Table tab = table_->tcal().table() ;
1408  ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
1409  ROTableColumn *sharedCol ;
1410  uInt oldnr = 0 ;
1411  uInt newnr = 0 ;
1412  TableRow row( tab ) ;
1413  TableRecord &trec = row.record() ;
1414  RecordFieldPtr<uInt> idRF( trec, "ID" ) ;
1415  RecordFieldPtr<String> timeRF( trec, "TIME" ) ;
1416  RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
1417  TableIterator iter0( sctabsel, "FEED_ID" ) ;
1418  while( !iter0.pastEnd() ) {
1419    Table t0 = iter0.table() ;
1420    sharedCol = tpoolr->construct( t0, "FEED_ID" ) ;
1421    Int feedId = sharedCol->asInt( 0 ) ;
1422    tpoolr->destroy( sharedCol ) ;
1423    TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
1424    while( !iter1.pastEnd() ) {
1425      Table t1 = iter1.table() ;
1426      sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ;
1427      Int spwId = sharedCol->asInt( 0 ) ;
1428      tpoolr->destroy( sharedCol ) ;
1429      tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ;
1430      ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ;
1431      Vector<uInt> idminmax( 2, oldnr ) ;
1432      for ( uInt irow = 0 ; irow < t1.nrow() ; irow++ ) {
1433        String sTime = MVTime( scTimeCol(irow) ).string( MVTime::YMD ) ;
1434        *timeRF = sTime ;
1435        uInt idx = oldnr ;
1436        Matrix<Float> subtcal = (*tmpTcalCol)( irow ) ;
1437        for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1438          *idRF = idx++ ;
1439          //*tcalRF = subtcal.row( ipol ) ;
1440          tcalRF.define( subtcal.row( ipol ) ) ;
1441
1442          // commit row
1443          tab.addRow() ;
1444          row.put( tab.nrow()-1 ) ;
1445
1446          newnr++ ;
1447        }
1448        idminmax[0] = oldnr ;
1449        idminmax[1] = newnr - 1 ;
1450        oldnr = newnr ;
1451
1452        String key = keyTcal( feedId, spwId, sTime ) ;
1453        tcalrec_.define( key, idminmax ) ;
1454      }
1455      delete tmpTcalCol ;
1456      iter1++ ;
1457    }
1458    iter0++ ;
1459  }
1460
1461  //tcalrec_.print( std::cout ) ;
1462//   double endSec = gettimeofday_sec() ;
1463//   os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1464}
1465
1466uInt MSFiller::getWeatherId( uInt idx, Double wtime )
1467{
1468//   double startSec = gettimeofday_sec() ;
1469//   os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
1470  uInt nrow = mwTime_.size() ;
1471  if ( nrow < 2 )
1472    return 0 ;
1473  uInt wid = nrow ;
1474  for ( uInt i = idx ; i < nrow-1 ; i++ ) {
1475    Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
1476    // use of INTERVAL column is problematic
1477    // since there are "blank" time of weather monitoring
1478    //Double tEnd = tStart + mwInterval_[i] ;
1479    Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
1480    //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1481    if ( wtime >= tStart && wtime <= tEnd ) {
1482      wid = i ;
1483      break ;
1484    }
1485  }
1486  if ( wid == nrow ) {
1487    uInt i = nrow - 1 ;
1488    Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
1489    Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
1490    //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1491    if ( wtime >= tStart && wtime <= tEnd )
1492      wid = i-1 ;
1493    else
1494      wid = i ;
1495  }
1496
1497  //if ( wid == nrow )
1498  //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
1499
1500//   double endSec = gettimeofday_sec() ;
1501//   os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1502  return wid ;
1503}
1504
1505void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )
1506{
1507//   double startSec = gettimeofday_sec() ;
1508//   os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
1509
1510  if ( !isSysCal_ )
1511    return ;
1512
1513  uInt nrow = tidx.nelements() ;
1514  if ( scTime.nelements() == 0 )
1515    return ;
1516  else if ( scTime.nelements() == 1 ) {
1517    tidx[0] = 0 ;
1518    return ;
1519  }
1520  uInt scnrow = scTime.nelements() ;
1521  uInt idx = 0 ;
1522  const Double half = 0.5e0 ;
1523  // execute  binary search
1524  idx = binarySearch( scTime, tcol[0].get( "s" ).getValue() ) ;
1525  if ( idx != 0 )
1526    idx -= 1 ;
1527  for ( uInt i = 0 ; i < nrow ; i++ ) {
1528    Double t = tcol[i].get( "s" ).getValue() ;
1529    Double tsc = scTime[0].get( "s" ).getValue() ;
1530    if ( t < tsc ) {
1531      tidx[i] = 0 ;
1532      continue ;
1533    }
1534    for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
1535      Double tsc1 = scTime[j].get( "s" ).getValue() ;
1536      Double dt1 = scInterval[j] ;
1537      Double tsc2 = scTime[j+1].get( "s" ).getValue() ;
1538      Double dt2 = scInterval[j+1] ;
1539      if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {
1540        tidx[i] = j ;
1541        idx = j ;
1542        break ;
1543      }
1544    }
1545    if ( tidx[i] == -1 ) {
1546//       Double tsc = scTime[scnrow-1].get( "s" ).getValue() ;
1547//       Double dt = scInterval[scnrow-1] ;
1548//       if ( t <= tsc+0.5*dt ) {
1549//         tidx[i] = scnrow-1 ;
1550//       }
1551      tidx[i] = scnrow-1 ;
1552    }
1553  }
1554//   double endSec = gettimeofday_sec() ;
1555//   os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
1556  return ;
1557}
1558
1559Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )
1560{
1561//   double startSec = gettimeofday_sec() ;
1562//   os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
1563  //if ( table_->tcal().table().nrow() == 0 ) {
1564  if ( !isSysCal_ ) {
1565    os_ << "No TCAL rows" << LogIO::POST ;
1566    Block<uInt> tcalids( 0 ) ;
1567    return  tcalids ;
1568  }   
1569  //String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
1570  String sctime = MVTime( t.getValue() ).string(MVTime::YMD) ;
1571  String key = keyTcal( fid, spwid, sctime ) ;
1572  if ( !tcalrec_.isDefined( key ) ) {
1573    os_ << "No TCAL rows" << LogIO::POST ;
1574    Block<uInt> tcalids( 0 ) ;
1575    return tcalids ;
1576  }
1577  Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;
1578  uInt npol = ids[1] - ids[0] + 1 ;
1579  Block<uInt> tcalids( npol ) ;
1580  tcalids[0] = ids[0] ;
1581  tcalids[1] = ids[1] ;
1582  for ( uInt ipol = 2 ; ipol < npol ; ipol++ )
1583    tcalids[ipol] = ids[0] + ipol - 1 ;
1584
1585//   double endSec = gettimeofday_sec() ;
1586//   os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1587  return tcalids ;
1588}
1589
1590uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t )
1591{
1592//   double startSec = gettimeofday_sec() ;
1593//   os_ << "start MSFiller::getDirection() startSec=" << startSec << LogIO::POST ;
1594  // assume that cols is sorted by TIME
1595  Bool doInterp = False ;
1596  //uInt nrow = cols.nrow() ;
1597  uInt nrow = tab.nrow() ;
1598  if ( nrow == 0 )
1599    return 0 ;
1600  ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ;
1601  ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ;
1602  ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ;
1603  // binary search if idx == 0
1604  if ( idx == 0 ) {
1605    uInt nrowb = 75000 ;
1606    if ( nrow > nrowb ) {
1607      uInt nblock = nrow / nrowb + 1 ;
1608      for ( uInt iblock = 0 ; iblock < nblock ; iblock++ ) {
1609        uInt high = min( nrowb, nrow-iblock*nrowb ) ;
1610
1611        if ( tcol( high-1 ).get( "s" ).getValue() < t ) {
1612          idx = iblock * nrowb ;
1613          continue ;
1614        }
1615
1616        Vector<MEpoch> tarr( high ) ;
1617        for ( uInt irow = 0 ; irow < high ; irow++ ) {
1618          tarr[irow] = tcol( iblock*nrowb+irow ) ;
1619        }
1620
1621        uInt bidx = binarySearch( tarr, t ) ;
1622
1623        idx = iblock * nrowb + bidx ;
1624        break ;
1625      }
1626    }
1627    else {
1628      Vector<MEpoch> tarr( nrow ) ;
1629      for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1630        tarr[irow] = tcol( irow ) ;
1631      }
1632      idx = binarySearch( tarr, t ) ;
1633    }
1634  }
1635  // ensure that tcol(idx) < t
1636  //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
1637  while ( tcol(idx).get("s").getValue() > t && idx > 0 )
1638    idx-- ;
1639  //os_ << "idx = " << idx << LogIO::POST ;
1640
1641  // index search
1642  for ( uInt i = idx ; i < nrow ; i++ ) {
1643    Double tref = tcol( i ).get( "s" ).getValue() ;
1644    if ( tref == t ) {
1645      idx = i ;
1646      break ;
1647    }
1648    else if ( tref > t ) {
1649      if ( i == 0 ) {
1650        idx = i ;
1651      }
1652      else {
1653        idx = i-1 ;
1654        doInterp = True ;
1655      }
1656      break ;
1657    }
1658    else {
1659      idx = nrow - 1 ;
1660    }
1661  }
1662  //os_ << "searched idx = " << idx << LogIO::POST ;
1663
1664  //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
1665  IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
1666  //os_ << "ip = " << ip << LogIO::POST ;
1667  ref = dmcol(idx)(ip).getRefString() ;
1668  //os_ << "ref = " << ref << LogIO::POST ;
1669  if ( doInterp ) {
1670    //os_ << "do interpolation" << LogIO::POST ;
1671    //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
1672    Double tref0 = tcol(idx).get("s").getValue() ;
1673    Double tref1 = tcol(idx+1).get("s").getValue() ;
1674    Matrix<Double> mdir0 = dcol( idx ) ;
1675    Matrix<Double> mdir1 = dcol( idx+1 ) ;
1676    Vector<Double> dir0 = mdir0.column( 0 ) ;
1677    //os_ << "dir0 = " << dir0 << LogIO::POST ;
1678    Vector<Double> dir1 = mdir1.column( 0 ) ;
1679    //os_ << "dir1 = " << dir1 << LogIO::POST ;
1680    Double dt0 = t - tref0 ;
1681    Double dt1 = tref1 - t ;
1682    dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;
1683    if ( mdir0.ncolumn() > 1 ) {
1684      if ( dt0 >= dt1 )
1685        srate.reference( mdir0.column( 1 ) ) ;
1686      else
1687        srate.reference( mdir1.column( 1 ) ) ;
1688    }
1689    //os_ << "dir = " << dir << LogIO::POST ;
1690  }
1691  else {
1692    //os_ << "no interpolation" << LogIO::POST ;
1693    Matrix<Double> mdir0 = dcol( idx ) ;
1694    dir.reference( mdir0.column( 0 ) ) ;
1695    if ( mdir0.ncolumn() > 1 )
1696      srate.reference( mdir0.column( 1 ) ) ;
1697  }
1698
1699//   double endSec = gettimeofday_sec() ;
1700//   os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1701  return idx ;
1702}
1703
1704String MSFiller::keyTcal( Int feedid, Int spwid, String stime )
1705{
1706  String sfeed = "FEED" + String::toString( feedid ) ;
1707  String sspw = "SPW" + String::toString( spwid ) ;
1708  return sfeed+":"+sspw+":"+stime ;
1709}
1710
1711uInt MSFiller::binarySearch( Vector<MEpoch> &timeList, Double target )
1712{
1713  Int low = 0 ;
1714  Int high = timeList.nelements() ;
1715  uInt idx = 0 ;
1716
1717  while ( low <= high ) {
1718    idx = (Int)( 0.5 * ( low + high ) ) ;
1719    Double t = timeList[idx].get( "s" ).getValue() ;
1720    if ( t < target )
1721      low = idx + 1 ;
1722    else if ( t > target )
1723      high = idx - 1 ;
1724    else
1725      return idx ;
1726  }
1727
1728  idx = max( 0, min( low, high ) ) ;
1729
1730  return idx ;
1731 
1732}
1733
1734} ;
1735
Note: See TracBrowser for help on using the repository browser.