Changeset 2240 for branches


Ignore:
Timestamp:
07/22/11 12:26:10 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

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...

Merge changes for speed up of MSFiller.
Revisions to be merged are r2217, r2219, r2221-2224, r2226, r2228,
r2231-2232, r2234, and r2237-2239.


Location:
branches/parallel
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • branches/parallel

  • branches/parallel/Makefile

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/SConstruct

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/apps

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/external-alma

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/external-alma/atnf/pks/pks_maths.cc

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/getsvnrev.sh

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/python

    • Property svn:mergeinfo changed (with no actual effect on merging)
  • branches/parallel/python/flagtoolbar.py

  • branches/parallel/src

  • branches/parallel/src/MSFiller.cpp

    r2191 r2240  
    9595{
    9696  os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
    97 //   double startSec = gettimeofday_sec() ;
    98 //   os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
     97  //double startSec = gettimeofday_sec() ;
     98  //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
    9999  //os_ << "   filename = " << filename << endl ;
    100100
     
    147147  isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
    148148
    149 //   double endSec = gettimeofday_sec() ;
    150 //   os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     149  //double endSec = gettimeofday_sec() ;
     150  //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    151151  return true ;
    152152}
     
    155155{
    156156  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 ;
     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 ;
    162162
    163163  // Initialize header
    164164  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 = "" ;
     165  initHeader( sdh ) ;
    182166 
    183167  // check if optional table exists
     
    249233
    250234  // SUBTABLES: FREQUENCIES
    251   table_->frequencies().setFrame( "LSRK" ) ;
    252   table_->frequencies().setFrame( "LSRK", True ) ;
     235  //string freqFrame = getFrame() ;
     236  string freqFrame = "LSRK" ;
     237  table_->frequencies().setFrame( freqFrame ) ;
     238  table_->frequencies().setFrame( freqFrame, True ) ;
    253239
    254240  // SUBTABLES: WEATHER
     
    266252  // SUBTABLES: HISTORY
    267253  //fillHistory() ;
    268 
    269   // shared pointers
    270   ROTableColumn *tcolr ;
    271254
    272255  // MAIN
     
    278261  MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
    279262  if ( getPt_ ) {
    280     //pointtab = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;
    281263    pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
    282264  }
    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 ) ;
     265  ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ;
     266  ROArrayColumn<Double> pdcol( pointtab, "DIRECTION" ) ;
     267  String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ;
     268  String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ;
    289269  sdh.antennaposition.resize( 3 ) ;
    290270  for ( int i = 0 ; i < 3 ; i++ )
     
    292272  String telescopeName = "" ;
    293273
    294 //   double time1 = gettimeofday_sec() ;
    295 //   os_ << "end fill init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     274  //double time1 = gettimeofday_sec() ;
     275  //os_ << "end init fill: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    296276
    297277  // row based
     
    313293  RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ;
    314294  RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ;
    315 
     295  RecordFieldPtr<Int> refbRF( trec, "REFBEAMNO" ) ;
     296  RecordFieldPtr<Int> fitidRF( trec, "FIT_ID" ) ;
     297  RecordFieldPtr<Float> tauRF( trec, "OPACITY" ) ;
     298  RecordFieldPtr<uInt> beamRF( trec, "BEAMNO" ) ;
     299  RecordFieldPtr<uInt> focusidRF( trec, "FOCUS_ID" ) ;
     300  RecordFieldPtr<uInt> ifnoRF( trec, "IFNO" ) ;
     301  RecordFieldPtr<uInt> molidRF( trec, "MOLECULE_ID" ) ;
     302  RecordFieldPtr<uInt> freqidRF( trec, "FREQ_ID" ) ;
     303  RecordFieldPtr<uInt> scanRF( trec, "SCANNO" ) ;
     304  RecordFieldPtr<Int> srctypeRF( trec, "SRCTYPE" ) ;
     305  RecordFieldPtr<String> srcnameRF( trec, "SRCNAME" ) ;
     306  RecordFieldPtr<String> fieldnameRF( trec, "FIELDNAME" ) ;
     307  RecordFieldPtr< Array<Double> > srcpmRF( trec, "SRCPROPERMOTION" ) ;
     308  RecordFieldPtr< Array<Double> > srcdirRF( trec, "SRCDIRECTION" ) ;
     309  RecordFieldPtr<Double> sysvelRF( trec, "SRCVELOCITY" ) ;
    316310
    317311  // REFBEAMNO
    318   RecordFieldPtr<Int> intRF( trec, "REFBEAMNO" ) ;
    319   *intRF = 0 ;
     312  *refbRF = -1 ;
    320313
    321314  // FIT_ID
    322   intRF.attachToRecord( trec, "FIT_ID" ) ;
    323   *intRF = -1 ;
     315  *fitidRF = -1 ;
    324316
    325317  // OPACITY
    326   RecordFieldPtr<Float> floatRF( trec, "OPACITY" ) ;
    327   *floatRF = 0.0 ;
     318  *tauRF = 0.0 ;
    328319
    329320  //
     
    332323  TableIterator iter0( mstable_, "OBSERVATION_ID" ) ;
    333324  while( !iter0.pastEnd() ) {
    334 //     time0 = gettimeofday_sec() ;
    335 //     os_ << "start 0th iteration: " << time0 << LogIO::POST ;
     325    //time0 = gettimeofday_sec() ;
     326    //os_ << "start 0th iteration: " << time0 << LogIO::POST ;
    336327    Table t0 = iter0.table() ;
    337     tcolr = tpoolr->construct( t0, "OBSERVATION_ID" ) ;
    338     Int obsId = tcolr->asInt( 0 ) ;
    339     tpoolr->destroy( tcolr ) ;
     328    Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ;
    340329    if ( sdh.observer == "" ) {
    341       tcolr = tpoolr->construct( obstab, "OBSERVER" ) ;
    342       sdh.observer = tcolr->asString( obsId ) ;
    343       tpoolr->destroy( tcolr ) ;
     330      sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ;
    344331    }
    345332    if ( sdh.project == "" ) {
    346       tcolr = tpoolr->construct( obstab, "PROJECT" ) ;
    347       sdh.observer = tcolr->asString( obsId ) ;
    348       tpoolr->destroy( tcolr ) ;
     333      sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ;
    349334    }
    350335    ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ;
     
    355340    }
    356341    if ( telescopeName == "" ) {
    357       tcolr = tpoolr->construct( obstab, "TELESCOPE_NAME" ) ;
    358       telescopeName = tcolr->asString( obsId ) ;
    359       tpoolr->destroy( tcolr ) ;
     342      telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ;
    360343    }
    361344    Int nbeam = 0 ;
    362 //     time1 = gettimeofday_sec() ;
    363 //     os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     345    //time1 = gettimeofday_sec() ;
     346    //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    364347    //
    365348    // ITERATION: FEED1
     
    367350    TableIterator iter1( t0, "FEED1" ) ;
    368351    while( !iter1.pastEnd() ) {
    369 //       time0 = gettimeofday_sec() ;
    370 //       os_ << "start 1st iteration: " << time0 << LogIO::POST ;
     352      //time0 = gettimeofday_sec() ;
     353      //os_ << "start 1st iteration: " << time0 << LogIO::POST ;
    371354      Table t1 = iter1.table() ;
    372355      // assume FEED1 == FEED2
    373       tcolr = tpoolr->construct( t1, "FEED1" ) ;
    374       Int feedId = tcolr->asInt( 0 ) ;
    375       tpoolr->destroy( tcolr ) ;
     356      Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ;
    376357      nbeam++ ;
    377358
    378359      // BEAMNO
    379       RecordFieldPtr<uInt> uintRF( trec, "BEAMNO" ) ;
    380       *uintRF = feedId ;
     360      *beamRF = feedId ;
    381361
    382362      // 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 ;
     363      *focusidRF = 0 ;
     364
     365      //time1 = gettimeofday_sec() ;
     366      //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    388367      //
    389368      // ITERATION: FIELD_ID
     
    391370      TableIterator iter2( t1, "FIELD_ID" ) ;
    392371      while( !iter2.pastEnd() ) {
    393 //         time0 = gettimeofday_sec() ;
    394 //         os_ << "start 2nd iteration: " << time0 << LogIO::POST ;
     372        //time0 = gettimeofday_sec() ;
     373        //os_ << "start 2nd iteration: " << time0 << LogIO::POST ;
    395374        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 ) ;
     375        Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ;
     376        Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ;
     377        String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ;
     378        fieldName += "__" + String::toString(fieldId) ;
    405379        ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ;
    406380        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          
     381        delete delayDirCol ;         
    413382
    414383        // 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 ;
     384        *fieldnameRF = fieldName ;
     385
     386
     387        //time1 = gettimeofday_sec() ;
     388        //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    421389        //
    422390        // ITERATION: DATA_DESC_ID
     
    424392        TableIterator iter3( t2, "DATA_DESC_ID" ) ;
    425393        while( !iter3.pastEnd() ) {
    426 //           time0 = gettimeofday_sec() ;
    427 //           os_ << "start 3rd iteration: " << time0 << LogIO::POST ;
     394          //time0 = gettimeofday_sec() ;
     395          //os_ << "start 3rd iteration: " << time0 << LogIO::POST ;
    428396          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 ) ;
     397          Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ;
     398          Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ;
     399          Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ;
    438400
    439401          // IFNO
    440           uintRF.attachToRecord( trec, "IFNO" ) ;
    441           *uintRF = (uInt)spwId ;
     402          *ifnoRF = (uInt)spwId ;
    442403
    443404          // polarization information
    444           tcolr = tpoolr->construct( poltab, "NUM_CORR" ) ;
    445           Int npol = tcolr->asInt( polId ) ;
    446           tpoolr->destroy( tcolr ) ;
     405          Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ;
    447406          ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ;
    448407          Vector<Int> corrtype = (*roArrICol)( polId ) ;
    449408          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           }
    458409          String srcName( "" ) ;
    459410          Vector<Double> srcPM( 2, 0.0 ) ;
    460411          Vector<Double> srcDir( 2, 0.0 ) ;
    461412          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           }
     413          Vector<Double> restFreqs ;
     414          Vector<String> transitionName ;
     415          Vector<Double> sysVels ;
     416//           os_ << "npol = " << npol << LogIO::POST ;
     417//           os_ << "corrtype = " << corrtype << LogIO::POST ;
     418
     419          // source and molecular transition
     420          sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ;
     421//           os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;
    494422
    495423          // SRCNAME
    496           strRF.attachToRecord( trec, "SRCNAME" ) ;
    497           *strRF = srcName ;
     424          *srcnameRF = srcName ;
    498425
    499426//           os_ << "srcName = " << srcName << LogIO::POST ;
    500427
    501428          // SRCPROPERMOTION
    502           RecordFieldPtr< Array<Double> > darrRF( trec, "SRCPROPERMOTION" ) ;
    503           *darrRF = srcPM ;
     429          *srcpmRF = srcPM ;
    504430
    505431          //os_ << "srcPM = " << srcPM << LogIO::POST ;
    506432
    507433          // SRCDIRECTION
    508           darrRF.attachToRecord( trec, "SRCDIRECTION" ) ;
    509           *darrRF = srcDir ;
     434          *srcdirRF = md.getAngle().getValue( "rad" ) ;
    510435
    511436          //os_ << "srcDir = " << srcDir << LogIO::POST ;
    512437
    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 ;
     438          // SRCVELOCITY
    519439          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 ;
     440          if ( !sysVels.empty() )
     441            sysVel = sysVels[0] ;
     442          *sysvelRF = sysVel ;
    551443
    552444//           os_ << "sysVel = " << sysVel << LogIO::POST ;
     
    555447
    556448          // MOLECULE_ID
    557           uintRF.attachToRecord( trec, "MOLECULE_ID" ) ;
    558           *uintRF = molId ;
     449          *molidRF = molId ;
    559450
    560451          // 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 ) ;
     452          uInt freqId ;
     453          Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ;
    567454          Bool iswvr = False ;
    568455          if ( nchan == 4 ) iswvr = True ;
    569           tpoolr->destroy( tcolr ) ;
    570           Bool even = False ;
    571           if ( (nchan/2)*2 == nchan ) even = True ;
    572456          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() ) {
     457          map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
     458          if ( iter == ifmap.end() ) {
     459            ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ;
     460            me = (*tmpMeasCol)( 0 ) ;
     461            delete tmpMeasCol ;
     462            Double refpix ;
     463            Double refval ;
     464            Double increment ;
     465            spectralSetup( spwId,
     466                           me,
     467                           mp,
     468                           md,
     469                           refpix,
     470                           refval,
     471                           increment,
     472                           nchan,
     473                           sdh.freqref,
     474                           sdh.reffreq,
     475                           sdh.bandwidth,
     476                           tpoolr ) ;
     477            freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
    629478            ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
    630479            //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
    631480          }
     481          else {
     482            freqId = iter->second ;
     483          }
    632484
    633485          // FREQ_ID
    634           uintRF.attachToRecord( trec, "FREQ_ID" ) ;
    635           *uintRF = freqId ;
     486          *freqidRF = freqId ;
    636487
    637488          // for TSYS and TCAL
     
    656507          if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
    657508
    658 //           time1 = gettimeofday_sec() ;
    659 //           os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     509          //time1 = gettimeofday_sec() ;
     510          //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    660511          //
    661512          // ITERATION: SCAN_NUMBER
     
    663514          TableIterator iter4( t3, "SCAN_NUMBER" ) ;
    664515          while( !iter4.pastEnd() ) {
    665 //             time0 = gettimeofday_sec() ;
    666 //             os_ << "start 4th iteration: " << time0 << LogIO::POST ;
     516            //time0 = gettimeofday_sec() ;
     517            //os_ << "start 4th iteration: " << time0 << LogIO::POST ;
    667518            Table t4 = iter4.table() ;
    668             tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ;
    669             Int scanNum = tcolr->asInt( 0 ) ;
    670             tpoolr->destroy( tcolr ) ;
     519            Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ;
    671520
    672521            // 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 ;
     522            *scanRF = scanNum - 1 ;
     523
     524            //time1 = gettimeofday_sec() ;
     525            //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    678526            //
    679527            // ITERATION: STATE_ID
     
    681529            TableIterator iter5( t4, "STATE_ID" ) ;
    682530            while( !iter5.pastEnd() ) {
    683 //               time0 = gettimeofday_sec() ;
    684 //               os_ << "start 5th iteration: " << time0 << LogIO::POST ;
     531              //time0 = gettimeofday_sec() ;
     532              //os_ << "start 5th iteration: " << time0 << LogIO::POST ;
    685533              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 ) ;
     534              Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ;
     535              String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ;
    692536              if ( sdh.obstype == "" ) sdh.obstype = obstype ;
    693537
    694538              Int nrow = t5.nrow() ;
    695 //               time1 = gettimeofday_sec() ;
    696 //               os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     539              //time1 = gettimeofday_sec() ;
     540              //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    697541
    698542              uInt cycle = 0 ;
     
    700544              Cube<Float> spArr ;
    701545              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                 }
     546              reshapeSpectraAndFlagtra( spArr,
     547                                        flArr,
     548                                        t5,
     549                                        npol,
     550                                        nchan,
     551                                        nrow,
     552                                        corrtype ) ;
     553              if ( sdh.fluxunit == "" ) {
     554                String colName = "FLOAT_DATA" ;
     555                if ( isData_ ) colName = "DATA" ;
     556                ROTableColumn dataCol( t5, colName ) ;
     557                const TableRecord &dataColKeys = dataCol.keywordSet() ;
     558                if ( dataColKeys.isDefined( "UNIT" ) )
     559                  sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
     560                else if ( dataColKeys.isDefined( "QuantumUnits" ) )
     561                  sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ;
    712562              }
    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               }
     563
    753564              ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
    754565              Block<MEpoch> mTimeB( nrow ) ;
    755566              for ( Int irow = 0 ; irow < nrow ; irow++ )
    756567                mTimeB[irow] = (*mTimeCol)( irow ) ;
    757               ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;
    758               ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;
     568//               ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;
     569//               ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;
    759570              Block<Int> sysCalIdx( nrow, -1 ) ;
    760571              if ( isSysCal_ ) {
     
    765576              Int srcType = getSrcType( stateId, tpoolr ) ;
    766577              uInt diridx = 0 ;
    767               MDirection::Types dirType ;
    768578              uInt wid = 0 ;
    769579              Int pidx = 0 ;
     
    789599             
    790600              // SRCTYPE
    791               intRF.attachToRecord( trec, "SRCTYPE" ) ;
    792               *intRF = srcType ;
     601              *srctypeRF = srcType ;
    793602
    794603              for ( Int irow = 0 ; irow < nrow ; irow++ ) {
     
    797606
    798607                // FLAGROW
    799                 *flrRF = (uInt)mFlagRowCol->asBool( irow ) ;
     608//                 *flrRF = (uInt)mFlagRowCol->asBool( irow ) ;
     609                *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ;
    800610
    801611                // SPECTRA, FLAG
     
    809619
    810620                // INTERVAL
    811                 *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ;
     621//                 *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ;
     622                *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ;
    812623
    813624                // TSYS
     
    825636
    826637                // WEATHER_ID
    827                 if ( isWeather_ )
     638                if ( isWeather_ ) {
    828639                  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                   }
     640                  *widRF = mwIndex_[wid] ;
    870641                }
    871642                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 ;
     643                  *widRF = wid ;
    888644                }
     645                 
     646
     647                // DIRECTION, AZEL, SCANRATE
     648                Vector<Double> dir ;
     649                Vector<Double> azel ;
     650                Vector<Double> scanrate ;
     651                String refString ;
     652                if ( getPt_ )
     653                  //diridx = getDirection( diridx, dir, azel, scanrate, pointtab, mTimeB[irow], mp ) ;
     654                  diridx = getDirection( diridx, dir, azel, scanrate, ptcol, pdcol, mTimeB[irow], mp ) ;
     655                else
     656                  getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ;
     657                *dirRF = dir ;
     658                *azRF = azel[0] ;
     659                *elRF = azel[1] ;
     660                *scrRF = scanrate ;
    889661
    890662                // Polarization dependent things
     
    893665                  *polnoRF = polnos[ipol] ;
    894666
    895                   //*spRF = sp.row( ipol ) ;
    896                   //*ucarrRF = fl.row( ipol ) ;
    897                   //*tsysRF = tsys.row( ipol ) ;
    898667                  spRF.define( sp.row( ipol ) ) ;
    899668                  ucarrRF.define( fl.row( ipol ) ) ;
     
    908677                cycle++ ;
    909678              }
    910               tpoolr->destroy( mIntervalCol ) ;
    911               tpoolr->destroy( mFlagRowCol ) ;
    912 
    913 //               time1 = gettimeofday_sec() ;
    914 //               os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     679             
     680              //time1 = gettimeofday_sec() ;
     681              //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    915682
    916683              iter5.next() ;
    917684            }
    918 
    919 
    920685            iter4.next() ;
    921686          }
    922 
    923 
    924687          iter3.next() ;
    925688        }
    926 
    927              
    928689        iter2.next() ;
    929690      }
    930 
    931 
    932691      iter1.next() ;
    933692    }
     
    974733    sdh.freqref = "SOURCE";
    975734  }
     735
     736  if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" )
     737    sdh.fluxunit = "K" ;
    976738  table_->setHeader( sdh ) ;
    977739
     
    1006768  }
    1007769
    1008 //   double endSec = gettimeofday_sec() ;
    1009 //   os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     770  //double endSec = gettimeofday_sec() ;
     771  //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1010772}
    1011773
     
    1020782Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
    1021783{
    1022 //   double startSec = gettimeofday_sec() ;
    1023 //   os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
     784  //double startSec = gettimeofday_sec() ;
     785  //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
    1024786
    1025787  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 ) ;
     788  String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ;
     789  Bool sig = asBool( "SIG", stateId, statetab, tpool ) ;
     790  Bool ref = asBool( "REF", stateId, statetab, tpool ) ;
     791  Double cal = asDouble( "CAL", stateId, statetab, tpool ) ;
    1039792  //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
    1040793
     
    1182935   
    1183936  //os_ << "srcType = " << srcType << LogIO::POST ;
    1184 //   double endSec = gettimeofday_sec() ;
    1185 //   os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     937  //double endSec = gettimeofday_sec() ;
     938  //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1186939  return srcType ;
    1187940}
     
    1190943Block<uInt> MSFiller::getPolNo( Int corrType )
    1191944{
    1192 //   double startSec = gettimeofday_sec() ;
    1193 //   os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
     945  //double startSec = gettimeofday_sec() ;
     946  //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
    1194947  Block<uInt> polno( 1 ) ;
    1195948
     
    1221974  }
    1222975  //os_ << "polno = " << polno << LogIO::POST ;
    1223 //   double endSec = gettimeofday_sec() ;
    1224 //   os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     976  //double endSec = gettimeofday_sec() ;
     977  //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1225978 
    1226979  return polno ;
     
    1229982String MSFiller::getPolType( Int corrType )
    1230983{
    1231 //   double startSec = gettimeofday_sec() ;
    1232 //   os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
     984  //double startSec = gettimeofday_sec() ;
     985  //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
    1233986  String poltype = "" ;
    1234987
     
    1242995    poltype = "linpol" ;
    1243996
    1244 //   double endSec = gettimeofday_sec() ;
    1245 //   os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     997  //double endSec = gettimeofday_sec() ;
     998  //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1246999  return poltype ;
    12471000}
     
    12491002void MSFiller::fillWeather()
    12501003{
    1251 //   double startSec = gettimeofday_sec() ;
    1252 //   os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
     1004  //double startSec = gettimeofday_sec() ;
     1005  //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
    12531006
    12541007  if ( !isWeather_ ) {
     
    12781031  wtab.addRow( wnrow ) ;
    12791032
     1033  Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ;
     1034  Int stationId = -1 ;
     1035  if ( stationInfoExists ) {
     1036    // determine which station is closer
     1037    ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
     1038    ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ;
     1039    Vector<Int> stationIds = stationCol.getColumn() ;
     1040    Vector<Int> stationIdList( 0 ) ;
     1041    Matrix<Double> stationPosList( 0, 3, 0.0 ) ;
     1042    uInt numStation = 0 ;
     1043    for ( uInt i = 0 ; i < stationIds.size() ; i++ ) {
     1044      if ( !anyEQ( stationIdList, stationIds[i] ) ) {
     1045        numStation++ ;
     1046        stationIdList.resize( numStation, True ) ;
     1047        stationIdList[numStation-1] = stationIds[i] ;
     1048        stationPosList.resize( numStation, 3, True ) ;
     1049        stationPosList.row( numStation-1 ) = stationPosCol( i ) ;
     1050      }
     1051    }
     1052    //os_ << "staionIdList = " << stationIdList << endl ;
     1053    Table mAntenna = mstable_.antenna() ;
     1054    ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ;
     1055    Vector<Double> antpos = antposCol( antenna_ ) ;
     1056    Double minDiff = -1.0 ;
     1057    for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) {
     1058      Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ;
     1059      if ( minDiff < 0.0 || minDiff > diff ) {
     1060        minDiff = diff ;
     1061        stationId = stationIdList[i] ;
     1062      }
     1063    }
     1064  }
     1065  //os_ << "stationId = " << stationId << endl ;
     1066 
    12801067  ScalarColumn<Float> *fCol ;
    12811068  ROScalarColumn<Float> *sharedFloatCol ;
     
    13221109  ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
    13231110  String tUnit = tqCol.getUnits() ;
    1324   mwTime_ = tCol.getColumn() ;
     1111  Vector<Double> mwTime = tCol.getColumn() ;
    13251112  if ( tUnit == "d" )
    1326     mwTime_ *= 86400.0 ;
     1113    mwTime *= 86400.0 ;
    13271114  tqCol.attach( mWeatherSel, "INTERVAL" ) ;
    13281115  tCol.attach( mWeatherSel, "INTERVAL" ) ;
    13291116  String iUnit = tqCol.getUnits() ;
    1330   mwInterval_ = tCol.getColumn() ;
     1117  Vector<Double> mwInterval = tCol.getColumn() ;
    13311118  if ( iUnit == "d" )
    1332     mwInterval_ *= 86400.0 ;
     1119    mwInterval *= 86400.0 ;
     1120
     1121  if ( stationId > 0 ) {
     1122    ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ;
     1123    Vector<Int> stationVec = stationCol.getColumn() ;
     1124    uInt wsnrow = ntrue( stationVec == stationId ) ;
     1125    mwTime_.resize( wsnrow ) ;
     1126    mwInterval_.resize( wsnrow ) ;
     1127    mwIndex_.resize( wsnrow ) ;
     1128    uInt wsidx = 0 ;
     1129    for ( uInt irow = 0 ; irow < wnrow ; irow++ ) {
     1130      if ( stationId == stationVec[irow] ) {
     1131        mwTime_[wsidx] = mwTime[irow] ;
     1132        mwInterval_[wsidx] = mwInterval[irow] ;
     1133        mwIndex_[wsidx] = irow ;
     1134        wsidx++ ;
     1135      }
     1136    }
     1137  }
     1138  else {
     1139    mwTime_ = mwTime ;
     1140    mwInterval_ = mwInterval ;
     1141    mwIndex_.resize( mwTime_.size() ) ;
     1142    indgen( mwIndex_ ) ;
     1143  }
    13331144  //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 ;
     1145  //double endSec = gettimeofday_sec() ;
     1146  //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    13361147}
    13371148
    13381149void MSFiller::fillFocus()
    13391150{
    1340 //   double startSec = gettimeofday_sec() ;
    1341 //   os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
     1151  //double startSec = gettimeofday_sec() ;
     1152  //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
    13421153  // 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 ;
     1154  table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ;
     1155  //double endSec = gettimeofday_sec() ;
     1156  //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    13491157}
    13501158
    13511159void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr )
    13521160{
    1353 //   double startSec = gettimeofday_sec() ;
    1354 //   os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
     1161  //double startSec = gettimeofday_sec() ;
     1162  //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
    13551163
    13561164  if ( !isSysCal_ ) {
     
    14071215  Table tab = table_->tcal().table() ;
    14081216  ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
    1409   ROTableColumn *sharedCol ;
    14101217  uInt oldnr = 0 ;
    14111218  uInt newnr = 0 ;
     
    14181225  while( !iter0.pastEnd() ) {
    14191226    Table t0 = iter0.table() ;
    1420     sharedCol = tpoolr->construct( t0, "FEED_ID" ) ;
    1421     Int feedId = sharedCol->asInt( 0 ) ;
    1422     tpoolr->destroy( sharedCol ) ;
     1227    Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ;
    14231228    TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
    14241229    while( !iter1.pastEnd() ) {
    14251230      Table t1 = iter1.table() ;
    1426       sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ;
    1427       Int spwId = sharedCol->asInt( 0 ) ;
    1428       tpoolr->destroy( sharedCol ) ;
     1231      Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ;
    14291232      tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ;
    14301233      ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ;
     
    14601263
    14611264  //tcalrec_.print( std::cout ) ;
    1462 //   double endSec = gettimeofday_sec() ;
    1463 //   os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1265  //double endSec = gettimeofday_sec() ;
     1266  //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    14641267}
    14651268
    14661269uInt MSFiller::getWeatherId( uInt idx, Double wtime )
    14671270{
    1468 //   double startSec = gettimeofday_sec() ;
    1469 //   os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
     1271  //double startSec = gettimeofday_sec() ;
     1272  //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
    14701273  uInt nrow = mwTime_.size() ;
    14711274  if ( nrow < 2 )
    14721275    return 0 ;
    14731276  uInt wid = nrow ;
     1277  if ( idx == 0 ) {
     1278    wid = 0 ;
     1279    Double tStart = mwTime_[wid]-0.5*mwInterval_[wid] ;
     1280    if ( wtime < tStart )
     1281      return wid ;
     1282  }
    14741283  for ( uInt i = idx ; i < nrow-1 ; i++ ) {
    14751284    Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
     
    14981307  //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
    14991308
    1500 //   double endSec = gettimeofday_sec() ;
    1501 //   os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1309  //double endSec = gettimeofday_sec() ;
     1310  //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    15021311  return wid ;
    15031312}
     
    15051314void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )
    15061315{
    1507 //   double startSec = gettimeofday_sec() ;
    1508 //   os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
     1316  //double startSec = gettimeofday_sec() ;
     1317  //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
    15091318
    15101319  if ( !isSysCal_ )
     
    15521361    }
    15531362  }
    1554 //   double endSec = gettimeofday_sec() ;
    1555 //   os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
     1363  //double endSec = gettimeofday_sec() ;
     1364  //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
    15561365  return ;
    15571366}
     
    15591368Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )
    15601369{
    1561 //   double startSec = gettimeofday_sec() ;
    1562 //   os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
     1370  //double startSec = gettimeofday_sec() ;
     1371  //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
    15631372  //if ( table_->tcal().table().nrow() == 0 ) {
    15641373  if ( !isSysCal_ ) {
     
    15831392    tcalids[ipol] = ids[0] + ipol - 1 ;
    15841393
    1585 //   double endSec = gettimeofday_sec() ;
    1586 //   os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1394  //double endSec = gettimeofday_sec() ;
     1395  //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    15871396  return tcalids ;
    15881397}
    15891398
    1590 uInt 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 ;
     1399uInt MSFiller::getDirection( uInt idx,
     1400                             Vector<Double> &dir,
     1401                             Vector<Double> &srate,
     1402                             String &ref,
     1403                             ROScalarColumn<Double> &tcol,
     1404                             ROArrayColumn<Double> &dcol,
     1405                             Double t )
     1406{
     1407  //double startSec = gettimeofday_sec() ;
     1408  //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ;
     1409  //double time0 = gettimeofday_sec() ;
     1410  //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ;
    15941411  // assume that cols is sorted by TIME
    15951412  Bool doInterp = False ;
    1596   //uInt nrow = cols.nrow() ;
    1597   uInt nrow = tab.nrow() ;
     1413  uInt nrow = tcol.nrow() ;
    15981414  if ( nrow == 0 )
    15991415    return 0 ;
    1600   ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ;
    1601   ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ;
    1602   ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ;
     1416  TableRecord trec = tcol.keywordSet() ;
     1417  String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ;
     1418  Double factor = 1.0 ;
     1419  if ( tUnit == "d" )
     1420    factor = 86400.0 ;
    16031421  // binary search if idx == 0
    16041422  if ( idx == 0 ) {
    1605     uInt nrowb = 75000 ;
     1423    uInt nrowb = 1000 ;
    16061424    if ( nrow > nrowb ) {
    16071425      uInt nblock = nrow / nrowb + 1 ;
     
    16091427        uInt high = min( nrowb, nrow-iblock*nrowb ) ;
    16101428
    1611         if ( tcol( high-1 ).get( "s" ).getValue() < t ) {
     1429        if ( tcol( high-1 ) * factor < t ) {
    16121430          idx = iblock * nrowb ;
    16131431          continue ;
    16141432        }
    16151433
    1616         Vector<MEpoch> tarr( high ) ;
    1617         for ( uInt irow = 0 ; irow < high ; irow++ ) {
    1618           tarr[irow] = tcol( iblock*nrowb+irow ) ;
    1619         }
     1434        RefRows refrows( iblock*nrowb, iblock*nrowb+high, 1 ) ;
     1435        Vector<Double> tarr = tcol.getColumnCells( refrows ) ;
     1436        if ( tUnit == "d" )
     1437          tarr *= factor ;
    16201438
    16211439        uInt bidx = binarySearch( tarr, t ) ;
     
    16261444    }
    16271445    else {
    1628       Vector<MEpoch> tarr( nrow ) ;
    1629       for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
    1630         tarr[irow] = tcol( irow ) ;
    1631       }
     1446      Vector<Double> tarr = tcol.getColumn() ;
     1447      if ( tUnit == "d" )
     1448        tarr *= factor ;
    16321449      idx = binarySearch( tarr, t ) ;
    16331450    }
    16341451  }
     1452  //double time1 = gettimeofday_sec() ;
     1453  //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    16351454  // ensure that tcol(idx) < t
    16361455  //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 )
     1456  //time0 = gettimeofday_sec() ;
     1457  //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ;
     1458  while( tcol( idx ) * factor > t && idx > 0 )
    16381459    idx-- ;
    16391460  //os_ << "idx = " << idx << LogIO::POST ;
     
    16411462  // index search
    16421463  for ( uInt i = idx ; i < nrow ; i++ ) {
    1643     Double tref = tcol( i ).get( "s" ).getValue() ;
     1464    //Double tref = tcol( i ).get( "s" ).getValue() ;
     1465    Double tref = tcol( i ) * factor ;
    16441466    if ( tref == t ) {
    16451467      idx = i ;
     
    16601482    }
    16611483  }
     1484  //time1 = gettimeofday_sec() ;
     1485  //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    16621486  //os_ << "searched idx = " << idx << LogIO::POST ;
    16631487
     1488  //time0 = gettimeofday_sec() ;
     1489  //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ;
    16641490  //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
    1665   IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
     1491  //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
     1492  IPosition ip( dcol(idx).shape().nelements(), 0 ) ;
    16661493  //os_ << "ip = " << ip << LogIO::POST ;
    1667   ref = dmcol(idx)(ip).getRefString() ;
     1494  //ref = dmcol(idx)(ip).getRefString() ;
     1495  trec = dcol.keywordSet() ;
     1496  Record rec = trec.asRecord( "MEASINFO" ) ;
     1497  ref = rec.asString( "Ref" ) ;
    16681498  //os_ << "ref = " << ref << LogIO::POST ;
    16691499  if ( doInterp ) {
    16701500    //os_ << "do interpolation" << LogIO::POST ;
    16711501    //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() ;
     1502//     Double tref0 = tcol(idx).get("s").getValue() ;
     1503//     Double tref1 = tcol(idx+1).get("s").getValue() ;
     1504    Double tref0 = tcol(idx) * factor ;
     1505    Double tref1 = tcol(idx+1) * factor ;
    16741506    Matrix<Double> mdir0 = dcol( idx ) ;
    16751507    Matrix<Double> mdir1 = dcol( idx+1 ) ;
     
    16971529  }
    16981530
    1699 //   double endSec = gettimeofday_sec() ;
    1700 //   os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1531  //time1 = gettimeofday_sec() ;
     1532  //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
     1533  //double endSec = gettimeofday_sec() ;
     1534  //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    17011535  return idx ;
    17021536}
     
    17221556    else if ( t > target )
    17231557      high = idx - 1 ;
    1724     else
     1558    else {
    17251559      return idx ;
     1560    }
    17261561  }
    17271562
     
    17291564
    17301565  return idx ;
     1566}
     1567
     1568uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target )
     1569{
     1570  Int low = 0 ;
     1571  Int high = timeList.nelements() ;
     1572  uInt idx = 0 ;
     1573
     1574  while ( low <= high ) {
     1575    idx = (Int)( 0.5 * ( low + high ) ) ;
     1576    Double t = timeList[idx] ;
     1577    if ( t < target )
     1578      low = idx + 1 ;
     1579    else if ( t > target )
     1580      high = idx - 1 ;
     1581    else {
     1582      return idx ;
     1583    }
     1584  }
     1585
     1586  idx = max( 0, min( low, high ) ) ;
     1587
     1588  return idx ;
     1589}
     1590
     1591string MSFiller::getFrame()
     1592{
     1593  MFrequency::Types frame = MFrequency::DEFAULT ;
     1594  ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ;
     1595  ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ;
     1596  uInt nrow = numChanCol.nrow() ;
     1597  Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ;
     1598  uInt nref = 0 ;
     1599  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
     1600    if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR
     1601      measFreqRef[nref] = measFreqRefCol.asInt( irow ) ;
     1602      nref++ ;
     1603    }
     1604  }
     1605  if ( nref > 0 )
     1606    frame = (MFrequency::Types)measFreqRef[0] ;
     1607
     1608  return MFrequency::showType( frame ) ;
     1609}
     1610
     1611void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp,
     1612                                         Cube<Bool> &fl,
     1613                                         Table &tab,
     1614                                         Int &npol,
     1615                                         Int &nchan,
     1616                                         Int &nrow,
     1617                                         Vector<Int> &corrtype )
     1618{
     1619  //double startSec = gettimeofday_sec() ;
     1620  //os_ << "start MSFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ; 
     1621  if ( isFloatData_ ) {
     1622    ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
     1623    ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ;
     1624    sp = mFloatDataCol.getColumn() ;
     1625    fl = mFlagCol.getColumn() ;
     1626  }
     1627  else if ( isData_ ) {
     1628    sp.resize( npol, nchan, nrow ) ;
     1629    fl.resize( npol, nchan, nrow ) ;
     1630    ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
     1631    ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ;
     1632    if ( npol < 3 ) {
     1633//       Cube<Complex> tmp = mDataCol.getColumn() ;
     1634//       Float *sp_p = sp.data() ;
     1635//       Complex *tmp_p = tmp.data() ;
     1636//       for ( uInt i = 0 ; i < sp.nelements() ; i++ )
     1637//         sp_p[i] = tmp_p[i].real() ;
     1638      Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ;
     1639      IPosition start( 3, 0, 0, 0 ) ;
     1640      IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ;
     1641      IPosition inc( 3, 2, 1, 1 ) ;
     1642      sp = tmp( start, end, inc ) ;
     1643      fl = mFlagCol.getColumn() ;
     1644    }
     1645    else {
     1646      for ( Int irow = 0 ; irow < nrow ; irow++ ) {
     1647        Bool crossOK = False ;
     1648        Matrix<Complex> mSp = mDataCol( irow ) ;
     1649        Matrix<Bool> mFl = mFlagCol( irow ) ;
     1650        Matrix<Float> spxy = sp.xyPlane( irow ) ;
     1651        Matrix<Bool> flxy = fl.xyPlane( irow ) ;
     1652        for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     1653          if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
     1654               || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
     1655            if ( !crossOK ) {
     1656              Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
     1657              IPosition start( 1, 0 ) ;
     1658              IPosition end( 1, 2*nchan-1 ) ;
     1659              IPosition inc( 1, 2 ) ;
     1660              spxy.row( ipol ) = tmp( start, end, inc ) ;
     1661              flxy.row( ipol ) = mFl.row( ipol ) ;
     1662              start = IPosition( 1, 1 ) ;
     1663              spxy.row( ipol+1 ) = tmp( start, end, inc ) ;
     1664              flxy.row( ipol+1 ) = mFl.row( ipol ) ;
     1665              if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) {
     1666                spxy.row( ipol+1 ) = spxy.row( ipol+1 ) * (Float)-1.0 ;
     1667              }
     1668              crossOK = True ;
     1669            }
     1670          }
     1671          else {
     1672            Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
     1673            IPosition start( 1, 0 ) ;
     1674            IPosition end( 1, 2*nchan-1 ) ;
     1675            IPosition inc( 1, 2 ) ;
     1676            spxy.row( ipol ) = tmp( start, end, inc ) ;
     1677            flxy.row( ipol ) = mFl.row( ipol ) ;
     1678          }
     1679        }
     1680      }
     1681    }
     1682  }
     1683  //double endSec = gettimeofday_sec() ;
     1684  //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1685}
     1686
     1687uInt MSFiller::getDirection( uInt idx,
     1688                             Vector<Double> &dir,
     1689                             Vector<Double> &azel,
     1690                             Vector<Double> &srate,
     1691                             ROScalarColumn<Double> &ptcol,
     1692                             ROArrayColumn<Double> &pdcol,
     1693                             MEpoch &t,
     1694                             MPosition &antpos )
     1695{
     1696  //double startSec = gettimeofday_sec() ;
     1697  //os_ << "start MSFiller::getDirection2() startSec=" << startSec << LogIO::POST ; 
     1698  String refString ;
     1699  MDirection::Types dirType ;
     1700  uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ;
     1701  MDirection::getType( dirType, refString ) ;
     1702  MeasFrame mf( t, antpos ) ;
     1703  if ( refString == "J2000" ) {
     1704    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
     1705    azel = toazel( dir ).getAngle("rad").getValue() ;
     1706  }
     1707  else if ( refString(0,4) == "AZEL" ) {
     1708    azel = dir.copy() ;
     1709    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     1710    dir = toj2000( dir ).getAngle("rad").getValue() ;
     1711  }
     1712  else {
     1713    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
     1714    azel = toazel( dir ).getAngle("rad").getValue() ;
     1715    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     1716    dir = toj2000( dir ).getAngle("rad").getValue() ;
     1717  }
     1718  if ( srate.size() == 0 ) {
     1719    srate.resize( 2 ) ;
     1720    srate = 0.0 ;
     1721  }
     1722  //double endSec = gettimeofday_sec() ;
     1723  //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1724  return diridx ;
     1725}
     1726
     1727void MSFiller::getSourceDirection( Vector<Double> &dir,
     1728                                   Vector<Double> &azel,
     1729                                   Vector<Double> &srate,
     1730                                   MEpoch &t,
     1731                                   MPosition &antpos,
     1732                                   Vector<MDirection> &srcdir )
     1733{
     1734  //double startSec = gettimeofday_sec() ;
     1735  //os_ << "start MSFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ;
     1736  Vector<Double> defaultScanrate( 2, 0.0 ) ;
     1737  Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ;
     1738  if ( srcdir.nelements() > 1 )
     1739    defaultScanrate = srcdir[1].getAngle( "rad" ).getValue() ;
     1740  String ref = srcdir[0].getRefString() ;
     1741  MDirection::Types dirType ;
     1742  MDirection::getType( dirType, ref ) ;
     1743  MeasFrame mf( t, antpos ) ;
     1744  if ( ref != "J2000" ) {
     1745    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     1746    dir = toj2000( defaultDir ).getAngle("rad").getValue() ;
     1747  }
     1748  else
     1749    dir = defaultDir ;
     1750  if ( ref != "AZELGEO" ) {
     1751    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
     1752    azel = toazel( defaultDir ).getAngle("rad").getValue() ;
     1753  }
     1754  else
     1755    azel = defaultDir ;
     1756  srate = defaultScanrate ;
     1757  //double endSec = gettimeofday_sec() ;
     1758  //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1759}
     1760
     1761void MSFiller::initHeader( STHeader &header )
     1762{
     1763  header.nchan = 0 ;
     1764  header.npol = 0 ;
     1765  header.nif = 0 ;
     1766  header.nbeam = 0 ;
     1767  header.observer = "" ;
     1768  header.project = "" ;
     1769  header.obstype = "" ;
     1770  header.antennaname = "" ;
     1771  header.antennaposition.resize( 0 ) ;
     1772  header.equinox = 0.0 ;
     1773  header.freqref = "" ;
     1774  header.reffreq = -1.0 ;
     1775  header.bandwidth = 0.0 ;
     1776  header.utc = 0.0 ;
     1777  header.fluxunit = "" ;
     1778  header.epoch = "" ;
     1779  header.poltype = "" ;
     1780}
     1781
     1782String MSFiller::asString( String name,
     1783                           uInt idx,
     1784                           Table tab,
     1785                           boost::object_pool<ROTableColumn> *pool )
     1786{
     1787  ROTableColumn *col = pool->construct( tab, name ) ;
     1788  String v = col->asString( idx ) ;
     1789  pool->destroy( col ) ;
     1790  return v ;
     1791}
     1792
     1793Bool MSFiller::asBool( String name,
     1794                         uInt idx,
     1795                         Table &tab,
     1796                         boost::object_pool<ROTableColumn> *pool )
     1797{
     1798  ROTableColumn *col = pool->construct( tab, name ) ;
     1799  Bool v = col->asBool( idx ) ;
     1800  pool->destroy( col ) ;
     1801  return v ;
     1802}
     1803
     1804uInt MSFiller::asuInt( String name,
     1805                         uInt idx,
     1806                         Table &tab,
     1807                         boost::object_pool<ROTableColumn> *pool )
     1808{
     1809  ROTableColumn *col = pool->construct( tab, name ) ;
     1810  uInt v = col->asuInt( idx ) ;
     1811  pool->destroy( col ) ;
     1812  return v ;
     1813}
     1814
     1815Int MSFiller::asInt( String name,
     1816                        uInt idx,
     1817                        Table &tab,
     1818                        boost::object_pool<ROTableColumn> *pool )
     1819{
     1820  ROTableColumn *col = pool->construct( tab, name ) ;
     1821  Int v = col->asInt( idx ) ;
     1822  pool->destroy( col ) ;
     1823  return v ;
     1824}
     1825
     1826Float MSFiller::asFloat( String name,
     1827                          uInt idx,
     1828                          Table &tab,
     1829                          boost::object_pool<ROTableColumn> *pool )
     1830{
     1831  ROTableColumn *col = pool->construct( tab, name ) ;
     1832  Float v = col->asfloat( idx ) ;
     1833  pool->destroy( col ) ;
     1834  return v ;
     1835}
     1836
     1837Double MSFiller::asDouble( String name,
     1838                           uInt idx,
     1839                           Table &tab,
     1840                           boost::object_pool<ROTableColumn> *pool )
     1841{
     1842  ROTableColumn *col = pool->construct( tab, name ) ;
     1843  Double v = col->asdouble( idx ) ;
     1844  pool->destroy( col ) ;
     1845  return v ;
     1846}
     1847
     1848void MSFiller::sourceInfo( Int sourceId,
     1849                           Int spwId,
     1850                           String &name,
     1851                           MDirection &direction,
     1852                           Vector<casa::Double> &properMotion,
     1853                           Vector<casa::Double> &restFreqs,
     1854                           Vector<casa::String> &transitions,
     1855                           Vector<casa::Double> &sysVels,
     1856                           boost::object_pool<ROTableColumn> *tpoolr )
     1857{
     1858  //double startSec = gettimeofday_sec() ;
     1859  //os_ << "start MSFiller::sourceInfo() startSec=" << startSec << LogIO::POST ;
     1860
     1861  MSSource srctab = mstable_.source() ;
     1862  MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;
     1863  if ( srctabSel.nrow() == 0 ) {
     1864    srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;
     1865  }
     1866  Int numLines = 0 ;
     1867  if ( srctabSel.nrow() > 0 ) {
     1868    // source name
     1869    name = asString( "NAME", 0, srctabSel, tpoolr ) ;
     1870   
     1871    // source proper motion
     1872    ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ;
     1873    properMotion = roArrDCol( 0 ) ;
     1874   
     1875    // source direction as MDirection object
     1876    ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ;
     1877    direction = tmpMeasCol( 0 ) ;
     1878   
     1879    // number of lines
     1880    numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ;
     1881  }
     1882  else {
     1883    name = "" ;
     1884    properMotion = Vector<Double>( 2, 0.0 ) ;
     1885    direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
     1886  }
     1887
     1888  restFreqs.resize( numLines ) ;
     1889  transitions.resize( numLines ) ;
     1890  sysVels.resize( numLines ) ;
     1891  if ( numLines > 0 ) {
     1892    if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) {
     1893      ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ;
     1894      Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ;
     1895      for ( int i = 0 ; i < numLines ; i++ ) {
     1896        restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
     1897      }
     1898    }
     1899    //os_ << "restFreqs = " << restFreqs << LogIO::POST ;
     1900    if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
     1901      ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ;
     1902      if ( transitionCol.isDefined( 0 ) )
     1903        transitions = transitionCol( 0 ) ;
     1904      //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ;
     1905    }
     1906    if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) {
     1907      ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ;
     1908      sysVels = roArrDCol( 0 ) ;
     1909    }
     1910  }
    17311911 
     1912  //double endSec = gettimeofday_sec() ;
     1913  //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1914}
     1915
     1916void MSFiller::spectralSetup( Int spwId,
     1917                              MEpoch &me,
     1918                              MPosition &mp,
     1919                              MDirection &md,
     1920                              Double &refpix,
     1921                              Double &refval,
     1922                              Double &increment,
     1923                              Int &nchan,
     1924                              String &freqref,
     1925                              Double &reffreq,
     1926                              Double &bandwidth,
     1927                              boost::object_pool<ROTableColumn> *tpoolr )
     1928{
     1929  //double startSec = gettimeofday_sec() ;
     1930  //os_ << "start MSFiller::spectralSetup() startSec=" << startSec << LogIO::POST ;
     1931
     1932  MSSpectralWindow spwtab = mstable_.spectralWindow() ;
     1933  MeasFrame mf( me, mp, md ) ;
     1934  MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ;
     1935  Bool even = False ;
     1936  if ( (nchan/2)*2 == nchan ) even = True ;
     1937  ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ;
     1938  Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ;
     1939  if ( nchan != 4 )
     1940    bandwidth = max( bandwidth, totbw ) ;
     1941  if ( freqref == "" && nchan != 4)
     1942    //sdh.freqref = MFrequency::showType( freqRef ) ;
     1943    freqref = "LSRK" ;
     1944  if ( reffreq == -1.0 && nchan != 4 ) {
     1945    tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ;
     1946    Quantum<Double> qreffreq = tmpQuantCol( spwId ) ;
     1947    if ( freqRef == MFrequency::LSRK ) {
     1948      reffreq = qreffreq.getValue("Hz") ;
     1949    }
     1950    else {
     1951      MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
     1952      reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
     1953    }
     1954  }
     1955  Int refchan = nchan / 2 ;
     1956  IPosition refip( 1, refchan ) ;
     1957  refpix = 0.5*(nchan-1) ;
     1958  refval = 0.0 ;
     1959  ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ;
     1960  increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ;
     1961  //           os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
     1962  sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ;
     1963  Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ;
     1964  if ( nchan > 1 && chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") )
     1965    increment *= -1.0 ;
     1966  if ( freqRef == MFrequency::LSRK ) {
     1967    if ( even ) {
     1968      IPosition refip0( 1, refchan-1 ) ;
     1969      Double refval0 = chanFreqs(refip0).getValue("Hz") ;
     1970      Double refval1 = chanFreqs(refip).getValue("Hz") ;
     1971      refval = 0.5 * ( refval0 + refval1 ) ;
     1972    }
     1973    else {
     1974      refval = chanFreqs(refip).getValue("Hz") ;
     1975    }
     1976  }
     1977  else {
     1978    MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
     1979    if ( even ) {
     1980      IPosition refip0( 1, refchan-1 ) ;
     1981      Double refval0 = chanFreqs(refip0).getValue("Hz") ;
     1982      Double refval1 = chanFreqs(refip).getValue("Hz") ;
     1983      refval = 0.5 * ( refval0 + refval1 ) ;
     1984      refval = tolsr( refval ).get("Hz").getValue() ;
     1985    }
     1986    else {
     1987      refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
     1988    }
     1989  }
     1990 
     1991  //double endSec = gettimeofday_sec() ;
     1992  //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    17321993}
    17331994
  • branches/parallel/src/MSFiller.h

    r2021 r2240  
    9090
    9191  // get direction for DIRECTION, AZIMUTH, and ELEVATION columns
    92   casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &srate, casa::String &ref, casa::MSPointing &tab, casa::Double t ) ;
     92  casa::uInt getDirection( casa::uInt idx,
     93                           casa::Vector<casa::Double> &dir,
     94                           casa::Vector<casa::Double> &srate,
     95                           casa::String &ref,
     96                           casa::ROScalarColumn<casa::Double> &ptcol,
     97                           casa::ROArrayColumn<casa::Double> &pdcol,
     98                           casa::Double t ) ;
     99  casa::uInt getDirection( casa::uInt idx,
     100                           casa::Vector<casa::Double> &dir,
     101                           casa::Vector<casa::Double> &azel,
     102                           casa::Vector<casa::Double> &srate,
     103                           casa::ROScalarColumn<casa::Double> &ptcol,
     104                           casa::ROArrayColumn<casa::Double> &pdcol,
     105                           casa::MEpoch &t, casa::MPosition &antpos ) ;
     106  void getSourceDirection( casa::Vector<casa::Double> &dir,
     107                           casa::Vector<casa::Double> &azel,
     108                           casa::Vector<casa::Double> &srate,
     109                           casa::MEpoch &t,
     110                           casa::MPosition &antpos,
     111                           casa::Vector<casa::MDirection> &srcdir ) ;
    93112
    94113  // create key for TCAL table
     
    97116  // binary search
    98117  casa::uInt binarySearch( casa::Vector<casa::MEpoch> &timeList, casa::Double target ) ;
     118  casa::uInt binarySearch( casa::Vector<casa::Double> &timeList, casa::Double target ) ;
    99119
    100120  // tool for HPC
    101121  double gettimeofday_sec() ;
    102  
     122
     123  // get frequency frame
     124  std::string getFrame() ;
     125
     126  // reshape SPECTRA and FLAGTRA
     127  void reshapeSpectraAndFlagtra( casa::Cube<casa::Float> &sp,
     128                                 casa::Cube<casa::Bool> &fl,
     129                                 casa::Table &tab,
     130                                 casa::Int &npol,
     131                                 casa::Int &nchan,
     132                                 casa::Int &nrow,
     133                                 casa::Vector<casa::Int> &corrtype ) ;
     134 
     135  // initialize header
     136  void initHeader( STHeader &header ) ;
     137
     138  // get value from table column using object pool
     139  casa::String asString( casa::String name,
     140                         casa::uInt idx,
     141                         casa::Table tab,
     142                         boost::object_pool<casa::ROTableColumn> *pool ) ;
     143  casa::Bool asBool( casa::String name,
     144                     casa::uInt idx,
     145                     casa::Table &tab,
     146                     boost::object_pool<casa::ROTableColumn> *pool ) ;
     147  casa::Int asInt( casa::String name,
     148                   casa::uInt idx,
     149                   casa::Table &tab,
     150                   boost::object_pool<casa::ROTableColumn> *pool ) ;
     151  casa::uInt asuInt( casa::String name,
     152                     casa::uInt idx,
     153                     casa::Table &tab,
     154                     boost::object_pool<casa::ROTableColumn> *pool ) ;
     155  casa::Float asFloat( casa::String name,
     156                       casa::uInt idx,
     157                       casa::Table &tab,
     158                       boost::object_pool<casa::ROTableColumn> *pool ) ;
     159  casa::Double asDouble( casa::String name,
     160                         casa::uInt idx,
     161                         casa::Table &tab,
     162                         boost::object_pool<casa::ROTableColumn> *pool ) ;
     163
     164  void sourceInfo( casa::Int sourceId,
     165                   casa::Int spwId,
     166                   casa::String &name,
     167                   casa::MDirection &direction,
     168                   casa::Vector<casa::Double> &properMotion,
     169                   casa::Vector<casa::Double> &restFreqs,
     170                   casa::Vector<casa::String> &transitions,
     171                   casa::Vector<casa::Double> &sysVels,
     172                   boost::object_pool<casa::ROTableColumn> *pool ) ;
     173
     174  void spectralSetup( casa::Int spwId,
     175                      casa::MEpoch &me,
     176                      casa::MPosition &mp,
     177                      casa::MDirection &md,
     178                      casa::Double &refpix,
     179                      casa::Double &refval,
     180                      casa::Double &increment,
     181                      casa::Int &nchan,
     182                      casa::String &freqref,
     183                      casa::Double &reffreq,
     184                      casa::Double &bandwidth,
     185                      boost::object_pool<casa::ROTableColumn> *pool ) ;
     186
    103187  casa::CountedPtr<Scantable> table_ ;
    104188  casa::MeasurementSet mstable_ ;
     
    126210  casa::Vector<casa::Double> mwTime_ ;
    127211  casa::Vector<casa::Double> mwInterval_ ;
     212  casa::Vector<casa::uInt> mwIndex_ ;
    128213
    129214  // Record for TCAL_ID
  • branches/parallel/src/SConscript

    • Property svn:mergeinfo changed (with no actual effect on merging)
Note: See TracChangeset for help on using the changeset viewer.