source: trunk/src/MSFiller.cpp @ 2022

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

New Development: No

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

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