Changeset 2291


Ignore:
Timestamp:
09/12/11 12:07:41 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: MSFillerUtils and MSWriterUtils added

Test Programs: sd regressions, test_sdsave

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote MSFiller and MSWriter based on TableTraverse.
TableTraverse is changed a bit since it doesn't work on plain table.


Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/MSFiller.cpp

    r2289 r2291  
     1
    12//
    23// C++ Interface: MSFiller
     
    45// Description:
    56//
    6 // This class is specific filler for MS format
     7// This class is specific filler for MS format
     8// New version that is implemented using TableVisitor instead of TableIterator
    79//
    8 // Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
     10// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2011
    911//
    1012// Copyright: See COPYING file that comes with this distribution
     
    1214//
    1315
     16#include <assert.h>
    1417#include <iostream>
    1518#include <map>
     
    5154#include "STHeader.h"
    5255
    53 // #include <ctime>
    54 // #include <sys/time.h>
    55 
    5656#include "MathUtils.h"
    5757
     
    6060
    6161namespace asap {
    62 // double MSFiller::gettimeofday_sec()
    63 // {
    64 //   struct timeval tv ;
    65 //   gettimeofday( &tv, NULL ) ;
    66 //   return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
    67 // }
     62
     63class BaseMSFillerVisitor: public TableVisitor {
     64  uInt lastRecordNo ;
     65  Int lastObservationId ;
     66  Int lastFeedId ;
     67  Int lastFieldId ;
     68  Int lastDataDescId ;
     69  Int lastScanNo ;
     70  Int lastStateId ;
     71  Double lastTime ;
     72protected:
     73  const Table &table;
     74  uInt count;
     75public:
     76  BaseMSFillerVisitor(const Table &table)
     77   : table(table)
     78  {
     79    count = 0;
     80  }
     81 
     82  virtual void enterObservationId(const uInt recordNo, Int columnValue) { }
     83  virtual void leaveObservationId(const uInt recordNo, Int columnValue) { }
     84  virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
     85  virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
     86  virtual void enterFieldId(const uInt recordNo, Int columnValue) { }
     87  virtual void leaveFieldId(const uInt recordNo, Int columnValue) { }
     88  virtual void enterDataDescId(const uInt recordNo, Int columnValue) { }
     89  virtual void leaveDataDescId(const uInt recordNo, Int columnValue) { }
     90  virtual void enterScanNo(const uInt recordNo, Int columnValue) { }
     91  virtual void leaveScanNo(const uInt recordNo, Int columnValue) { }
     92  virtual void enterStateId(const uInt recordNo, Int columnValue) { }
     93  virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
     94  virtual void enterTime(const uInt recordNo, Double columnValue) { }
     95  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     96
     97  virtual Bool visitRecord(const uInt recordNo,
     98                           const Int ObservationId,
     99                           const Int feedId,
     100                           const Int fieldId,
     101                           const Int dataDescId,
     102                           const Int scanNo,
     103                           const Int stateId,
     104                           const Double time) { return True ; }
     105
     106  virtual Bool visit(Bool isFirst, const uInt recordNo,
     107                     const uInt nCols, void const *const colValues[]) {
     108    Int observationId, feedId, fieldId, dataDescId, scanNo, stateId;
     109    Double time;
     110    { // prologue
     111      uInt i = 0;
     112      {
     113        const Int *col = (const Int *)colValues[i++];
     114        observationId = col[recordNo];
     115      }
     116      {
     117        const Int *col = (const Int *)colValues[i++];
     118        feedId = col[recordNo];
     119      }
     120      {
     121        const Int *col = (const Int *)colValues[i++];
     122        fieldId = col[recordNo];
     123      }
     124      {
     125        const Int *col = (const Int *)colValues[i++];
     126        dataDescId = col[recordNo];
     127      }
     128      {
     129        const Int *col = (const Int *)colValues[i++];
     130        scanNo = col[recordNo];
     131      }
     132      {
     133        const Int *col = (const Int *)colValues[i++];
     134        stateId = col[recordNo];
     135      }
     136      {
     137        const Double *col = (const Double *)colValues[i++];
     138        time = col[recordNo];
     139      }
     140      assert(nCols == i);
     141    }
     142
     143    if (isFirst) {
     144      enterObservationId(recordNo, observationId);
     145      enterFeedId(recordNo, feedId);
     146      enterFieldId(recordNo, fieldId);
     147      enterDataDescId(recordNo, dataDescId);
     148      enterScanNo(recordNo, scanNo);
     149      enterStateId(recordNo, stateId);
     150      enterTime(recordNo, time);
     151    } else {
     152      if (lastObservationId != observationId) {
     153        leaveTime(lastRecordNo, lastTime);
     154        leaveStateId(lastRecordNo, lastStateId);
     155        leaveScanNo(lastRecordNo, lastScanNo);
     156        leaveDataDescId(lastRecordNo, lastDataDescId);
     157        leaveFieldId(lastRecordNo, lastFieldId);
     158        leaveFeedId(lastRecordNo, lastFeedId);
     159        leaveObservationId(lastRecordNo, lastObservationId);
     160
     161        enterObservationId(recordNo, observationId);
     162        enterFeedId(recordNo, feedId);
     163        enterFieldId(recordNo, fieldId);
     164        enterDataDescId(recordNo, dataDescId);
     165        enterScanNo(recordNo, scanNo);
     166        enterStateId(recordNo, stateId);
     167        enterTime(recordNo, time);
     168      } else if (lastFeedId != feedId) {
     169        leaveTime(lastRecordNo, lastTime);
     170        leaveStateId(lastRecordNo, lastStateId);
     171        leaveScanNo(lastRecordNo, lastScanNo);
     172        leaveDataDescId(lastRecordNo, lastDataDescId);
     173        leaveFieldId(lastRecordNo, lastFieldId);
     174        leaveFeedId(lastRecordNo, lastFeedId);
     175
     176        enterFeedId(recordNo, feedId);
     177        enterFieldId(recordNo, fieldId);
     178        enterDataDescId(recordNo, dataDescId);
     179        enterScanNo(recordNo, scanNo);
     180        enterStateId(recordNo, stateId);
     181        enterTime(recordNo, time);
     182      } else if (lastFieldId != fieldId) {
     183        leaveTime(lastRecordNo, lastTime);
     184        leaveStateId(lastRecordNo, lastStateId);
     185        leaveScanNo(lastRecordNo, lastScanNo);
     186        leaveDataDescId(lastRecordNo, lastDataDescId);
     187        leaveFieldId(lastRecordNo, lastFieldId);
     188
     189        enterFieldId(recordNo, fieldId);
     190        enterDataDescId(recordNo, dataDescId);
     191        enterScanNo(recordNo, scanNo);
     192        enterStateId(recordNo, stateId);
     193        enterTime(recordNo, time);
     194      } else if (lastDataDescId != dataDescId) {
     195        leaveTime(lastRecordNo, lastTime);
     196        leaveStateId(lastRecordNo, lastStateId);
     197        leaveScanNo(lastRecordNo, lastScanNo);
     198        leaveDataDescId(lastRecordNo, lastDataDescId);
     199
     200        enterDataDescId(recordNo, dataDescId);
     201        enterScanNo(recordNo, scanNo);
     202        enterStateId(recordNo, stateId);
     203        enterTime(recordNo, time);
     204      } else if (lastScanNo != scanNo) {
     205        leaveTime(lastRecordNo, lastTime);
     206        leaveStateId(lastRecordNo, lastStateId);
     207        leaveScanNo(lastRecordNo, lastScanNo);
     208
     209        enterScanNo(recordNo, scanNo);
     210        enterStateId(recordNo, stateId);
     211        enterTime(recordNo, time);
     212      } else if (lastStateId != stateId) {
     213        leaveTime(lastRecordNo, lastTime);
     214        leaveStateId(lastRecordNo, lastStateId);
     215
     216        enterStateId(recordNo, stateId);
     217        enterTime(recordNo, time);
     218      } else if (lastTime != time) {
     219        leaveTime(lastRecordNo, lastTime);
     220        enterTime(recordNo, time);
     221      }
     222    }
     223    count++;
     224    Bool result = visitRecord(recordNo, observationId, feedId, fieldId, dataDescId,
     225                              scanNo, stateId, time);
     226
     227    { // epilogue
     228      lastRecordNo = recordNo;
     229
     230      lastObservationId = observationId;
     231      lastFeedId = feedId;
     232      lastFieldId = fieldId;
     233      lastDataDescId = dataDescId;
     234      lastScanNo = scanNo;
     235      lastStateId = stateId;
     236      lastTime = time;
     237    }
     238    return result ;
     239  }
     240
     241  virtual void finish() {
     242    if (count > 0) {
     243      leaveTime(lastRecordNo, lastTime);
     244      leaveStateId(lastRecordNo, lastStateId);
     245      leaveScanNo(lastRecordNo, lastScanNo);
     246      leaveDataDescId(lastRecordNo, lastDataDescId);
     247      leaveFieldId(lastRecordNo, lastFieldId);
     248      leaveFeedId(lastRecordNo, lastFeedId);
     249      leaveObservationId(lastRecordNo, lastObservationId);
     250    }
     251  }
     252};
     253
     254class MSFillerVisitor: public BaseMSFillerVisitor, public MSFillerUtils {
     255public:
     256  MSFillerVisitor(const Table &from, Scantable &to)
     257    : BaseMSFillerVisitor(from),
     258      scantable( to )
     259  {
     260    antennaId = 0 ;
     261    rowidx = 0 ;
     262    tablerow = TableRow( scantable.table() ) ;
     263    feedEntry = Vector<Int>( 64, -1 ) ;
     264    nbeam = 0 ;
     265    ifmap.clear() ;
     266    const TableDesc &desc = table.tableDesc() ;
     267    if ( desc.isColumn( "DATA" ) )
     268      dataColumnName = "DATA" ;
     269    else if ( desc.isColumn( "FLOAT_DATA" ) )
     270      dataColumnName = "FLOAT_DATA" ;
     271    getpt = False ;
     272    isWeather = False ;
     273    isSysCal = False ;
     274    cycleNo = 0 ;
     275    numSysCalRow = 0 ;
     276    header = scantable.getHeader() ;
     277    fluxUnit( header.fluxunit ) ;
     278
     279    // MS subtables
     280    const TableRecord &hdr = table.keywordSet();
     281    obstab = hdr.asTable( "OBSERVATION" ) ;
     282    sctab = hdr.asTable( "SYSCAL" ) ;
     283    spwtab = hdr.asTable( "SPECTRAL_WINDOW" ) ;
     284    statetab = hdr.asTable( "STATE" ) ;
     285    ddtab = hdr.asTable( "DATA_DESCRIPTION" ) ;
     286    poltab = hdr.asTable( "POLARIZATION" ) ;
     287    fieldtab = hdr.asTable( "FIELD" ) ;
     288    anttab = hdr.asTable( "ANTENNA" ) ;
     289    srctab = hdr.asTable( "SOURCE" ) ;
     290
     291    // attach to columns
     292    // MS MAIN
     293    intervalCol.attach( table, "INTERVAL" ) ;
     294    flagRowCol.attach( table, "FLAG_ROW" ) ;
     295    flagCol.attach( table, "FLAG" ) ;
     296    if ( dataColumnName.compare( "DATA" ) == 0 )
     297      dataCol.attach( table, dataColumnName ) ;
     298    else
     299      floatDataCol.attach( table, dataColumnName ) ;
     300
     301    // set dummy epoch
     302    mf.set( currentTime ) ;
     303
     304    // initialize dirType
     305    dirType = MDirection::N_Types ;
     306
     307    //
     308    // add rows to scantable
     309    //
     310    // number of polarization is up to 4
     311    uInt addrow = table.nrow() * maxNumPol() ;
     312    scantable.table().addRow( addrow ) ;
     313
     314    // attach to columns
     315    // Scantable MAIN
     316    TableRecord &r = tablerow.record() ;
     317    timeRF.attachToRecord( r, "TIME" ) ;
     318    intervalRF.attachToRecord( r, "INTERVAL" ) ;
     319    directionRF.attachToRecord( r, "DIRECTION" ) ;
     320    azimuthRF.attachToRecord( r, "AZIMUTH" ) ;
     321    elevationRF.attachToRecord( r, "ELEVATION" ) ;
     322    scanRateRF.attachToRecord( r, "SCANRATE" ) ;
     323    weatherIdRF.attachToRecord( r, "WEATHER_ID" ) ;
     324    cycleNoRF.attachToRecord( r, "CYCLENO" ) ;
     325    flagRowRF.attachToRecord( r, "FLAGROW" ) ;
     326    polNoRF.attachToRecord( r, "POLNO" ) ;
     327    tcalIdRF.attachToRecord( r, "TCAL_ID" ) ;
     328    spectraRF.attachToRecord( r, "SPECTRA" ) ;
     329    flagtraRF.attachToRecord( r, "FLAGTRA" ) ;
     330    tsysRF.attachToRecord( r, "TSYS" ) ;
     331    beamNoRF.attachToRecord( r, "BEAMNO" ) ;
     332    ifNoRF.attachToRecord( r, "IFNO" ) ;
     333    freqIdRF.attachToRecord( r, "FREQ_ID" ) ;
     334    moleculeIdRF.attachToRecord( r, "MOLECULE_ID" ) ;
     335    sourceNameRF.attachToRecord( r, "SRCNAME" ) ;
     336    sourceProperMotionRF.attachToRecord( r, "SRCPROPERMOTION" ) ;
     337    sourceDirectionRF.attachToRecord( r, "SRCDIRECTION" ) ;
     338    sourceVelocityRF.attachToRecord( r, "SRCVELOCITY" ) ;
     339    focusIdRF.attachToRecord( r, "FOCUS_ID" ) ;
     340    fieldNameRF.attachToRecord( r, "FIELDNAME" ) ;
     341    sourceTypeRF.attachToRecord( r, "SRCTYPE" ) ;
     342    scanNoRF.attachToRecord( r, "SCANNO" ) ;
     343
     344    // put values
     345    RecordFieldPtr<Int> refBeamNoRF( r, "REFBEAMNO" ) ;
     346    *refBeamNoRF = -1 ;
     347    RecordFieldPtr<Int> fitIdRF( r, "FIT_ID" ) ;
     348    *fitIdRF = -1 ;
     349    RecordFieldPtr<Float> opacityRF( r, "OPACITY" ) ;
     350    *opacityRF = 0.0 ;
     351  }
     352
     353  virtual void enterObservationId(const uInt recordNo, Int columnValue) {
     354    //printf("%u: ObservationId: %d\n", recordNo, columnValue);
     355    // update header
     356    if ( header.observer.empty() )
     357      getScalar( String("OBSERVER"), (uInt)columnValue, obstab, header.observer ) ;
     358    if ( header.project.empty() )
     359      getScalar( "PROJECT", (uInt)columnValue, obstab, header.project ) ;
     360    if ( header.utc == 0.0 ) {
     361      Vector<MEpoch> amp ;
     362      getArrayMeas( "TIME_RANGE", (uInt)columnValue, obstab, amp ) ;
     363      header.utc = amp[0].get( "d" ).getValue() ;
     364    }
     365    if ( header.antennaname.empty() )
     366      getScalar( "TELESCOPE_NAME", (uInt)columnValue, obstab, header.antennaname ) ;
     367  }
     368  virtual void leaveObservationId(const uInt recordNo, Int columnValue) {
     369    // update header
     370    header.nbeam = max( header.nbeam, (Int)nbeam ) ;
     371
     372    nbeam = 0 ;
     373    feedEntry = -1 ;
     374  }
     375  virtual void enterFeedId(const uInt recordNo, Int columnValue) {
     376    //printf("%u: FeedId: %d\n", recordNo, columnValue);
     377
     378    // update feed entry
     379    if ( allNE( feedEntry, columnValue ) ) {
     380      feedEntry[nbeam] = columnValue ;
     381      nbeam++ ;
     382    }
     383
     384    // put values
     385    *beamNoRF = (uInt)columnValue ;
     386    *focusIdRF = (uInt)0 ;
     387  }
     388  virtual void leaveFeedId(const uInt recordNo, Int columnValue) {
     389    Int nelem = (Int)feedEntry.nelements() ;
     390    if ( nbeam > nelem ) {
     391      feedEntry.resize( nelem+64, True ) ;
     392      Slicer slice( IPosition( 1, nelem ), IPosition( 1, feedEntry.nelements()-1 ) ) ;
     393      feedEntry( slice ) = -1 ;
     394    }
     395  }
     396  virtual void enterFieldId(const uInt recordNo, Int columnValue) {
     397    //printf("%u: FieldId: %d\n", recordNo, columnValue);
     398    // update sourceId and fieldName
     399    getScalar( "SOURCE_ID", (uInt)columnValue, fieldtab, sourceId ) ;
     400    String fieldName ;
     401    getScalar( "NAME", (uInt)columnValue, fieldtab, fieldName ) ;
     402    fieldName += "__" + String::toString( columnValue ) ;
     403
     404    // put values
     405    *fieldNameRF = fieldName ;
     406  }
     407  virtual void leaveFieldId(const uInt recordNo, Int columnValue) {
     408    sourceId = -1 ;
     409  }
     410  virtual void enterDataDescId(const uInt recordNo, Int columnValue) {
     411    //printf("%u: DataDescId: %d\n", recordNo, columnValue);
     412    // update polarization and spectral window ids
     413    getScalar( "POLARIZATION_ID", (uInt)columnValue, ddtab, polId ) ;
     414    getScalar( "SPECTRAL_WINDOW_ID", (uInt)columnValue, ddtab, spwId ) ;
     415
     416    // polarization setup
     417    getScalar( "NUM_CORR", (uInt)polId, poltab, npol ) ;
     418    Vector<Int> corrtype ;
     419    getArray( "CORR_TYPE", (uInt)polId, poltab, corrtype ) ;
     420    polnos = getPolNos( corrtype ) ;
     421   
     422    // process SOURCE table
     423    String sourceName ;
     424    Vector<Double> sourcePM, restFreqs, sysVels ;
     425    Vector<String> transition ;
     426    processSource( sourceId, spwId, sourceName, sourceDir, sourcePM,
     427                   restFreqs, transition, sysVels ) ;
     428
     429    // spectral setup
     430    uInt freqId ;
     431    Double reffreq, bandwidth ;
     432    String freqref ;
     433    getScalar( "NUM_CHAN", (uInt)spwId, spwtab, nchan ) ;
     434    Bool iswvr = (Bool)(nchan == 4) ;
     435    map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
     436    if ( iter == ifmap.end() ) {
     437      MEpoch me ;
     438      getScalarMeas( "TIME", recordNo, table, me ) ;
     439      spectralSetup( spwId, me, antpos, sourceDir,
     440                     freqId, nchan,
     441                     freqref, reffreq, bandwidth ) ;
     442      ifmap.insert( pair<Int,uInt>(spwId,freqId) ) ;
     443    }
     444    else {
     445      freqId = iter->second ;
     446    }
     447    sp.resize( npol, nchan ) ;
     448    fl.resize( npol, nchan ) ;
     449
     450
     451    // molecular setup
     452    STMolecules mtab = scantable.molecules() ;
     453    uInt molId = mtab.addEntry( restFreqs, transition, transition ) ;
     454
     455    // process SYSCAL table
     456    if ( isSysCal )
     457      processSysCal( spwId ) ;
     458
     459    // update header
     460    if ( !iswvr ) {
     461      header.nchan = max( header.nchan, nchan ) ;
     462      header.bandwidth = max( header.bandwidth, bandwidth ) ;
     463      if ( header.reffreq == -1.0 )
     464        header.reffreq = reffreq ;
     465      header.npol = max( header.npol, npol ) ;
     466      if ( header.poltype.empty() )
     467        header.poltype = getPolType( corrtype[0] ) ;
     468      if ( header.freqref.empty() )
     469        header.freqref = freqref ;
     470    }
     471   
     472    // put values
     473    *ifNoRF = (uInt)spwId ;
     474    *freqIdRF = freqId ;
     475    *moleculeIdRF = molId ;
     476    *sourceNameRF = sourceName ;
     477    sourceProperMotionRF.define( sourcePM ) ;
     478    Vector<Double> srcD = sourceDir.getAngle().getValue( "rad" ) ;
     479    sourceDirectionRF.define( srcD ) ;
     480    if ( !sysVels.empty() )
     481      *sourceVelocityRF = sysVels[0] ;
     482    else {
     483      *sourceVelocityRF = (Double)0.0 ;
     484    }
     485  }
     486  virtual void leaveDataDescId(const uInt recordNo, Int columnValue) {
     487    npol = 0 ;
     488    nchan = 0 ;
     489    numSysCalRow = 0 ;
     490  }
     491  virtual void enterScanNo(const uInt recordNo, Int columnValue) {
     492    //printf("%u: ScanNo: %d\n", recordNo, columnValue);
     493    // put value
     494    // scan number is 1-based in MS while 0-based in Scantable
     495    *scanNoRF = (uInt)columnValue - 1 ;
     496  }
     497  virtual void leaveScanNo(const uInt recordNo, Int columnValue) {
     498    cycleNo = 0 ;
     499  }
     500  virtual void enterStateId(const uInt recordNo, Int columnValue) {
     501    //printf("%u: StateId: %d\n", recordNo, columnValue);
     502    // SRCTYPE
     503    Int srcType = getSrcType( columnValue ) ;
     504   
     505    // update header
     506    if ( header.obstype.empty() )
     507      getScalar( "OBS_MODE", (uInt)columnValue, statetab, header.obstype ) ;
     508
     509    // put value
     510    *sourceTypeRF = srcType ;
     511  }
     512  virtual void leaveStateId(const uInt recordNo, Int columnValue) { }
     513  virtual void enterTime(const uInt recordNo, Double columnValue) {
     514    //printf("%u: Time: %f\n", recordNo, columnValue);
     515    currentTime = MEpoch( Quantity( columnValue, "s" ), MEpoch::UTC ) ;
     516
     517    // DIRECTION, AZEL, and SCANRATE
     518    Vector<Double> direction, azel ;
     519    Vector<Double> scanrate( 2, 0.0 ) ;
     520    if ( getpt )
     521      getDirection( direction, azel, scanrate ) ;
     522    else
     523      getSourceDirection( direction, azel, scanrate ) ;
     524
     525    // INTERVAL
     526    Double interval = intervalCol.asdouble( recordNo ) ;
     527
     528    // WEATHER_ID
     529    uInt wid = 0 ;
     530    if ( isWeather )
     531      wid = getWeatherId() ;
     532
     533    // put value
     534    Double t = currentTime.get( "d" ).getValue() ;
     535    *timeRF = t ;
     536    *intervalRF = interval ;
     537    directionRF.define( direction ) ;
     538    *azimuthRF = (Float)azel[0] ;
     539    *elevationRF = (Float)azel[1] ;
     540    scanRateRF.define( scanrate ) ;
     541    *weatherIdRF = wid ;
     542  }
     543  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     544  virtual Bool visitRecord(const uInt recordNo,
     545                           const Int observationId,
     546                           const Int feedId,
     547                           const Int fieldId,
     548                           const Int dataDescId,
     549                           const Int scanNo,
     550                           const Int stateId,
     551                           const Double time)
     552  {
     553    //printf("%u: %d, %d, %d, %d, %d, %d, %f\n", recordNo,
     554    //observationId, feedId, fieldId, dataDescId, scanNo, stateId, time);
     555
     556    // SPECTRA and FLAGTRA
     557    //Matrix<Float> sp;
     558    //Matrix<uChar> fl;
     559    spectraAndFlagtra( recordNo, sp, fl ) ;
     560
     561    // FLAGROW
     562    Bool flr = flagRowCol.asBool( recordNo ) ;
     563
     564    // TSYS
     565    Matrix<Float> tsys ;
     566    uInt scIdx = getSysCalIndex() ;
     567    if ( numSysCalRow > 0 ) {
     568      tsys = sysCalTsysCol( syscalRow[scIdx] ) ;
     569    }
     570    else {
     571      tsys.resize( npol, 1 ) ;
     572      tsys = 1.0 ;
     573    }
     574
     575    // TCAL_ID
     576    Block<uInt> tcalids( npol, 0 ) ;
     577    if ( numSysCalRow > 0 ) {
     578      tcalids = getTcalId( syscalTime[scIdx] ) ;
     579    }
     580
     581    // put value
     582    *cycleNoRF = cycleNo ;
     583    *flagRowRF = (uInt)flr ;
     584
     585    // for each polarization component
     586    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     587      // put value depending on polarization component
     588      *polNoRF = polnos[ipol] ;
     589      *tcalIdRF = tcalids[ipol] ;
     590      spectraRF.define( sp.row( ipol ) ) ;
     591      flagtraRF.define( fl.row( ipol ) ) ;
     592      tsysRF.define( tsys.row( ipol ) ) ;
     593     
     594      // commit row
     595      tablerow.put( rowidx ) ;
     596      rowidx++ ;
     597    }
     598
     599    // increment CYCLENO
     600    cycleNo++ ;
     601
     602    return True ;
     603  }
     604  virtual void finish()
     605  {
     606    BaseMSFillerVisitor::finish();
     607    //printf("Total: %u\n", count);
     608    // remove redundant rows
     609    //cout << "filled " << rowidx << " rows out of " << scantable.nrow() << " rows" << endl ;
     610    if ( scantable.nrow() > rowidx ) {
     611      uInt numRemove = scantable.nrow() - rowidx ;
     612      //cout << "numRemove = " << numRemove << endl ;
     613      Vector<uInt> rows( numRemove ) ;
     614      indgen( rows, rowidx ) ;
     615      scantable.table().removeRow( rows ) ;
     616    }
     617   
     618    // antenna name and station name
     619    String antennaName ;
     620    getScalar( "NAME", (uInt)antennaId, anttab, antennaName ) ;
     621    String stationName ;
     622    getScalar( "STATION", (uInt)antennaId, anttab, stationName ) ;
     623
     624    // update header
     625    header.nif = ifmap.size() ;
     626    header.antennaposition = antpos.get( "m" ).getValue() ;
     627    if ( header.antennaname.empty() || header.antennaname == antennaName )
     628      header.antennaname = antennaName ;
     629    else
     630      header.antennaname += "//" + antennaName ;
     631    if ( !stationName.empty() && stationName != antennaName )
     632      header.antennaname += "@" + stationName ;
     633    if ( header.fluxunit.empty() || header.fluxunit == "CNTS" )
     634      header.fluxunit = "K" ;
     635    header.epoch = "UTC" ;
     636    header.equinox = 2000.0 ;
     637    if (header.freqref == "TOPO") {
     638      header.freqref = "TOPOCENT";
     639    } else if (header.freqref == "GEO") {
     640      header.freqref = "GEOCENTR";
     641    } else if (header.freqref == "BARY") {
     642      header.freqref = "BARYCENT";
     643    } else if (header.freqref == "GALACTO") {
     644      header.freqref = "GALACTOC";
     645    } else if (header.freqref == "LGROUP") {
     646      header.freqref = "LOCALGRP";
     647    } else if (header.freqref == "CMB") {
     648      header.freqref = "CMBDIPOL";
     649    } else if (header.freqref == "REST") {
     650      header.freqref = "SOURCE";
     651    }
     652    scantable.setHeader( header ) ;
     653  }
     654  void setAntenna( Int id )
     655  {
     656    antennaId = id ;
     657
     658    Vector< Quantum<Double> > pos ;
     659    getArrayQuant( "POSITION", (uInt)antennaId, anttab, pos ) ;
     660    antpos = MPosition( MVPosition( pos ), MPosition::ITRF ) ;
     661    mf.set( antpos ) ;
     662  }
     663  void setPointingTable( const Table &tab, String columnToUse="DIRECTION" )
     664  {
     665    // input POINTING table must be
     666    //  1) selected by antenna
     667    //  2) sorted by TIME
     668    ROScalarColumn<Double> tcol( tab, "TIME" ) ;
     669    ROArrayColumn<Double> dcol( tab, columnToUse ) ;
     670    tcol.getColumn( pointingTime ) ;
     671    dcol.getColumn( pointingDirection ) ;
     672    const TableRecord &rec = dcol.keywordSet() ;
     673    String pointingRef = rec.asRecord( "MEASINFO" ).asString( "Ref" ) ;
     674    MDirection::getType( dirType, pointingRef ) ;
     675    getpt = True ;
     676  }
     677  void setWeatherTime( const Vector<Double> &t, const Vector<Double> &it )
     678  {
     679    isWeather = True ;
     680    weatherTime = t ;
     681    weatherInterval = it ;
     682  }
     683  void setSysCalRecord( const Record &r )
     684  //void setSysCalRecord( const map< String,Vector<uInt> > &r )
     685  {
     686    isSysCal = True ;
     687    syscalRecord = r ;
     688
     689    const TableDesc &desc = sctab.tableDesc() ;
     690    uInt nrow = sctab.nrow() ;
     691    syscalRow.resize( nrow ) ;
     692    syscalTime.resize( nrow ) ;
     693    syscalInterval.resize( nrow ) ;
     694    String tsysCol = "NONE" ;
     695    Vector<String> tsysCols = stringToVector( "TSYS_SPECTRUM,TSYS" ) ;
     696    for ( uInt i = 0 ; i < tsysCols.nelements() ; i++ ) {
     697      if ( tsysCol == "NONE" && desc.isColumn( tsysCols[i] ) )
     698        tsysCol = tsysCols[i] ;
     699    }
     700    sysCalTsysCol.attach( sctab, tsysCol ) ;
     701  }
     702  STHeader getHeader() { return header ; }
     703  uInt getNumBeam() { return nbeam ; }
     704  uInt getFilledRowNum() { return rowidx ; }
     705private:
     706  void fluxUnit( String &u )
     707  {
     708    ROTableColumn col( table, dataColumnName ) ;
     709    const TableRecord &rec = col.keywordSet() ;
     710    if ( rec.isDefined( "UNIT" ) )
     711      u = rec.asString( "UNIT" ) ;
     712    else if ( rec.isDefined( "QuantumUnits" ) )
     713      u = rec.asString( "QuantumUnits" ) ;
     714    if ( u.empty() )
     715      u = "K" ;
     716  }
     717  void processSource( Int sourceId, Int spwId,
     718                      String &name, MDirection &dir, Vector<Double> &pm,
     719                      Vector<Double> &rf, Vector<String> &trans, Vector<Double> &vel )
     720  {
     721    // find row
     722    uInt nrow = srctab.nrow() ;
     723    Int idx = -1 ;
     724    ROTableRow row( srctab ) ;
     725    for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
     726      const TableRecord &r = row.get( irow ) ;
     727      if ( r.asInt( "SOURCE_ID" ) == sourceId ) {
     728        Int tmpSpwId = r.asInt( "SPECTRAL_WINDOW_ID" ) ;
     729        if ( tmpSpwId == spwId || tmpSpwId == -1 ) {
     730          idx = (Int)irow ;
     731          break ;
     732        }
     733      }
     734    }
     735
     736    // fill
     737    Int numLines = 0 ;
     738    if ( idx != -1 ) {
     739      const TableRecord &r = row.get( idx ) ;
     740      name = r.asString( "NAME" ) ;
     741      getScalarMeas( "DIRECTION", idx, srctab, dir ) ;
     742      pm = r.toArrayDouble( "PROPER_MOTION" ) ;
     743      numLines = r.asInt( "NUM_LINES" ) ;
     744    }
     745    else {
     746      name = "" ;
     747      pm = Vector<Double>( 2, 0.0 ) ;
     748      dir = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
     749    }
     750    if ( !getpt ) {
     751      String ref = dir.getRefString() ;
     752      MDirection::getType( dirType, ref ) ;
     753    }
     754
     755    rf.resize( numLines ) ;
     756    trans.resize( numLines ) ;
     757    vel.resize( numLines ) ;
     758    if ( numLines > 0 ) {
     759      Block<Bool> isDefined = row.getDefined() ;
     760      Vector<String> colNames = row.columnNames() ;
     761      Vector<Int> indexes( 3, -1 ) ;
     762      Vector<String> cols = stringToVector( "REST_FREQUENCY,TRANSITION,SYSVEL" ) ;
     763      for ( uInt icol = 0 ; icol < colNames.nelements() ; icol++ ) {
     764        if ( anyEQ( indexes, -1 ) ) {
     765          for ( uInt jcol = 0 ; jcol < cols.nelements() ; jcol++ ) {
     766            if ( colNames[icol] == cols[jcol] )
     767              indexes[jcol] = icol ;
     768          }
     769        }
     770      }
     771      if ( indexes[0] != -1 && isDefined[indexes[0]] == True ) {
     772        Vector< Quantum<Double> > qrf ;
     773        getArrayQuant( "REST_FREQUENCY", idx, srctab, qrf ) ;
     774        for ( int i = 0 ; i < numLines ; i++ )
     775          rf[i] = qrf[i].getValue( "Hz" ) ;
     776      }
     777      if ( indexes[1] != -1 && isDefined[indexes[1]] == True ) {
     778        getArray( "TRANSITION", idx, srctab, trans ) ;
     779      }
     780      if ( indexes[2] != -1 && isDefined[indexes[2]] == True ) {
     781        Vector< Quantum<Double> > qsv ;
     782        getArrayQuant( "SYSVEL", idx, srctab, qsv ) ;
     783        for ( int i = 0 ; i < numLines ; i++ )
     784          vel[i] = qsv[i].getValue( "m/s" ) ;
     785      }
     786    }
     787  }
     788  void spectralSetup( Int &spwId, MEpoch &me, MPosition &mp, MDirection &md,
     789                      uInt &freqId, Int &nchan,
     790                      String &freqref, Double &reffreq, Double &bandwidth )
     791  {
     792    // fill
     793    Int measFreqRef ;
     794    getScalar( "MEAS_FREQ_REF", spwId, spwtab, measFreqRef ) ;
     795    MFrequency::Types freqRef = MFrequency::castType( measFreqRef ) ;
     796    //freqref = MFrequency::showType( freqRef ) ;
     797    freqref = "LSRK" ;
     798    Quantum<Double> q ;
     799    getScalarQuant( "TOTAL_BANDWIDTH", spwId, spwtab, q ) ;
     800    bandwidth = q.getValue( "Hz" ) ;
     801    getScalarQuant( "REF_FREQUENCY", spwId, spwtab, q ) ;
     802    reffreq = q.getValue( "Hz" ) ;
     803    Double refpix = 0.5 * ( (Double)nchan-1.0 ) ;
     804    Int refchan = ( nchan - 1 ) / 2 ;
     805    Bool even = (Bool)( nchan % 2 == 0 ) ;
     806    Vector< Quantum<Double> > qa ;
     807    getArrayQuant( "CHAN_WIDTH", spwId, spwtab, qa ) ;
     808    Double increment = qa[refchan].getValue( "Hz" ) ;
     809    getArrayQuant( "CHAN_FREQ", spwId, spwtab, qa ) ;
     810    if ( nchan == 1 ) {
     811      Int netSideband ;
     812      getScalar( "NET_SIDEBAND", spwId, spwtab, netSideband ) ;
     813      if ( netSideband == 1 ) increment *= -1.0 ;
     814    }
     815    else {
     816      if ( qa[0].getValue( "Hz" ) > qa[1].getValue( "Hz" ) )
     817        increment *= -1.0 ;
     818    }
     819    Double refval = qa[refchan].getValue( "Hz" ) ;
     820    if ( even )
     821      refval = 0.5 * ( refval + qa[refchan+1].getValue( "Hz" ) ) ;
     822    if ( freqRef != MFrequency::LSRK ) {
     823      MeasFrame mframe( me, mp, md ) ;
     824      MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mframe ) ) ;
     825      refval = tolsr( Quantum<Double>( refval, "Hz" ) ).get( "Hz" ).getValue() ;
     826    }
     827   
     828    // add new row to FREQUENCIES
     829    Table ftab = scantable.frequencies().table() ;
     830    freqId = ftab.nrow() ;
     831    ftab.addRow() ;
     832    TableRow row( ftab ) ;
     833    TableRecord &r = row.record() ;
     834    RecordFieldPtr<uInt> idRF( r, "ID" ) ;
     835    *idRF = freqId ;
     836    RecordFieldPtr<Double> refpixRF( r, "REFPIX" ) ;
     837    RecordFieldPtr<Double> refvalRF( r, "REFVAL" ) ;
     838    RecordFieldPtr<Double> incrRF( r, "INCREMENT" ) ;
     839    *refpixRF = refpix ;
     840    *refvalRF = refval ;
     841    *incrRF = increment ;
     842    row.put( freqId ) ;
     843  }
     844  void spectraAndFlagtra( uInt recordNo, Matrix<Float> &sp, Matrix<uChar> &fl )
     845  {
     846    Matrix<Bool> b = flagCol( recordNo ) ;
     847    if ( dataColumnName.compare( "FLOAT_DATA" ) == 0 ) {
     848      sp = floatDataCol( recordNo ) ;
     849      convertArray( fl, b ) ;
     850    }
     851    else {
     852      Bool notyet = True ;
     853      Matrix<Complex> c = dataCol( recordNo ) ;
     854      for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     855        if ( ( header.poltype == "linear" || header.poltype == "circular" )
     856             && ( polnos[ipol] == 2 || polnos[ipol] == 3 ) ) {
     857          if ( notyet ) {
     858            Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
     859            IPosition start( 1, 0 ) ;
     860            IPosition end( 1, 2*nchan-1 ) ;
     861            IPosition inc( 1, 2 ) ;
     862            if ( polnos[ipol] == 2 ) {
     863              sp.row( ipol ) = tmp( start, end, inc ) ;
     864              Vector<Bool> br = b.row( ipol ) ;
     865              Vector<uChar> flr = fl.row( ipol ) ;
     866              convertArray( flr, br ) ;
     867              start = IPosition( 1, 1 ) ;
     868              Int jpol = ipol+1 ;
     869              while( polnos[jpol] != 3 && jpol < npol )
     870                jpol++ ;
     871              sp.row( jpol ) = tmp( start, end, inc ) ;
     872              flr.reference( fl.row( jpol ) ) ;
     873              convertArray( flr, br ) ;
     874            }
     875            else if ( polnos[ipol] == 3 ) {
     876              sp.row( ipol ) = sp.row( ipol ) * (Float)(-1.0) ;
     877              Int jpol = ipol+1 ;
     878              while( polnos[jpol] != 2 && jpol < npol )
     879                jpol++ ;
     880              Vector<Bool> br = b.row( ipol ) ;
     881              Vector<uChar> flr = fl.row( jpol ) ;
     882              sp.row( jpol ) = tmp( start, end, inc ) ;
     883              convertArray( flr, br ) ;
     884              start = IPosition( 1, 1 ) ;
     885              sp.row( ipol ) = tmp( start, end, inc ) * (Float)(-1.0) ;
     886              flr.reference( fl.row( ipol ) ) ;
     887              convertArray( flr, br ) ;
     888            }
     889            notyet = False ;
     890          }
     891        }
     892        else {
     893          Vector<Float> tmp = ComplexToReal( c.row( ipol ) ) ;
     894          IPosition start( 1, 0 ) ;
     895          IPosition end( 1, 2*nchan-1 ) ;
     896          IPosition inc( 1, 2 ) ;
     897          sp.row( ipol ) = tmp( start, end, inc ) ;
     898          Vector<Bool> br = b.row( ipol ) ;
     899          Vector<uChar> flr = fl.row( ipol ) ;
     900          convertArray( flr, br ) ;
     901        }
     902      }
     903    }
     904  }
     905  uInt binarySearch( Vector<Double> &timeList, Double target )
     906  {
     907    Int low = 0 ;
     908    Int high = timeList.nelements() ;
     909    uInt idx = 0 ;
     910   
     911    while ( low <= high ) {
     912      idx = (Int)( 0.5 * ( low + high ) ) ;
     913      Double t = timeList[idx] ;
     914      if ( t < target )
     915        low = idx + 1 ;
     916      else if ( t > target )
     917        high = idx - 1 ;
     918      else {
     919        return idx ;
     920      }
     921    }
     922   
     923    idx = max( 0, min( low, high ) ) ;
     924    return idx ;
     925  }
     926  void getDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
     927  {
     928    // @todo At the moment, do binary search every time
     929    //       if this is bottleneck, frequency of binary search must be reduced
     930    Double t = currentTime.get( "s" ).getValue() ;
     931    uInt idx = min( binarySearch( pointingTime, t ), pointingTime.nelements()-1 ) ;
     932    Matrix<Double> d ;
     933    if ( pointingTime[idx] == t )
     934      d = pointingDirection.xyPlane( idx ) ;
     935    else if ( pointingTime[idx] < t ) {
     936      if ( idx == pointingTime.nelements()-1 )
     937        d = pointingDirection.xyPlane( idx ) ;
     938      else
     939        d = interp( pointingTime[idx], pointingTime[idx+1], t,
     940                    pointingDirection.xyPlane( idx ), pointingDirection.xyPlane( idx+1 ) ) ;
     941    }
     942    else {
     943      if ( idx == 0 )
     944        d = pointingDirection.xyPlane( idx ) ;
     945      else
     946        d = interp( pointingTime[idx-1], pointingTime[idx], t,
     947                    pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
     948    }
     949    mf.resetEpoch( currentTime ) ;
     950    Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
     951    if ( dirType != MDirection::J2000 ) {
     952      MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     953      dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
     954    }
     955    else {
     956      dir = d.column( 0 ) ;
     957    }
     958    if ( dirType != MDirection::AZELGEO ) {
     959      MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
     960      //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
     961      azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
     962    }
     963    else {
     964      azel = d.column( 0 ) ;
     965    }
     966    if ( d.ncolumn() > 1 )
     967      srate = d.column( 1 ) ;
     968  }
     969  void getSourceDirection( Vector<Double> &dir, Vector<Double> &azel, Vector<Double> &srate )
     970  {
     971    dir = sourceDir.getAngle( "rad" ).getValue() ;
     972    mf.resetEpoch( currentTime ) ;
     973    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
     974    azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
     975    if ( dirType != MDirection::J2000 ) {
     976      MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
     977      dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
     978    }
     979  }
     980  String detectSeparator( String &s )
     981  {
     982    String tmp = s.substr( 0, s.find_first_of( "," ) ) ;
     983    Char *separators[] = { ":", "#", ".", "_" } ;
     984    uInt nsep = 4 ;
     985    for ( uInt i = 0 ; i < nsep ; i++ ) {
     986      if ( tmp.find( separators[i] ) != String::npos )
     987        return separators[i] ;
     988    }
     989    return "" ;
     990  }
     991  Int getSrcType( Int stateId )
     992  {
     993    // get values
     994    Bool sig ;
     995    getScalar( "SIG", stateId, statetab, sig ) ;
     996    Bool ref ;
     997    getScalar( "REF", stateId, statetab, ref ) ;
     998    Double cal ;
     999    getScalar( "CAL", stateId, statetab, cal ) ;
     1000    String obsmode ;
     1001    getScalar( "OBS_MODE", stateId, statetab, obsmode ) ;
     1002    String sep = detectSeparator( obsmode ) ;
     1003   
     1004    Int srcType = SrcType::NOTYPE ;
     1005    if ( sep == ":" )
     1006      srcTypeGBT( srcType, sep, obsmode, sig, ref, cal ) ;
     1007    else if ( sep == "." || sep == "#" )
     1008      srcTypeALMA( srcType, sep, obsmode ) ;
     1009    else if ( sep == "_" )
     1010      srcTypeOldALMA( srcType, sep, obsmode, sig, ref ) ;
     1011    else
     1012      srcTypeDefault( srcType, sig, ref ) ;
     1013
     1014    return srcType ;
     1015  }
     1016  void srcTypeDefault( Int &st, Bool &sig, Bool &ref )
     1017  {
     1018    if ( sig ) st = SrcType::SIG ;
     1019    else if ( ref ) st = SrcType::REF ;
     1020  }
     1021  void srcTypeGBT( Int &st, String &sep, String &mode, Bool &sig, Bool &ref, Double &cal )
     1022  {
     1023    Int epos = mode.find_first_of( sep ) ;
     1024    Int nextpos = mode.find_first_of( sep, epos+1 ) ;
     1025    String m1 = mode.substr( 0, epos ) ;
     1026    String m2 = mode.substr( epos+1, nextpos-epos-1 ) ;
     1027    if ( m1 == "Nod" ) {
     1028      st = SrcType::NOD ;
     1029    }
     1030    else if ( m1 == "OffOn" ) {
     1031      if ( m2 == "PSWITCHON" ) st = SrcType::PSON ;
     1032      if ( m2 == "PSWITCHOFF" ) st = SrcType::PSOFF ;
     1033    }
     1034    else {
     1035      if ( m2 == "FSWITCH" ) {
     1036        if ( sig ) st = SrcType::FSON ;
     1037        else if ( ref ) st = SrcType::FSOFF ;
     1038      }
     1039    }
     1040    if ( cal > 0.0 ) {
     1041      if ( st == SrcType::NOD )
     1042        st = SrcType::NODCAL ;
     1043      else if ( st == SrcType::PSON )
     1044        st = SrcType::PONCAL ;
     1045      else if ( st == SrcType::PSOFF )
     1046        st = SrcType::POFFCAL ;
     1047      else if ( st == SrcType::FSON )
     1048        st = SrcType::FONCAL ;
     1049      else if ( st == SrcType::FSOFF )
     1050        st = SrcType::FOFFCAL ;
     1051      else
     1052        st = SrcType::CAL ;
     1053    }
     1054  }
     1055  void srcTypeALMA( Int &st, String &sep, String &mode )
     1056  {
     1057    Int epos = mode.find_first_of( "," ) ;
     1058    String first = mode.substr( 0, epos ) ;
     1059    epos = first.find_first_of( sep ) ;
     1060    Int nextpos = first.find_first_of( sep, epos+1 ) ;
     1061    String m1 = first.substr( 0, epos ) ;
     1062    String m2 = first.substr( epos+1, nextpos-epos-1 ) ;
     1063    if ( m1.find( "CALIBRATE_" ) == 0 ) {
     1064      if ( m2.find( "ON_SOURCE" ) == 0 )
     1065        st = SrcType::PONCAL ;
     1066      else if ( m2.find( "OFF_SOURCE" ) == 0 )
     1067        st = SrcType::POFFCAL ;
     1068    }
     1069    else if ( m1.find( "OBSERVE_TARGET" ) == 0 ) {
     1070      if ( m2.find( "ON_SOURCE" ) == 0 )
     1071        st = SrcType::PSON ;
     1072      else if ( m2.find( "OFF_SOURCE" ) == 0 )
     1073        st = SrcType::PSOFF ;
     1074    }
     1075  }
     1076  void srcTypeOldALMA( Int &st, String &sep, String &mode, Bool &sig, Bool &ref )
     1077  {
     1078    Int epos = mode.find_first_of( "," ) ;
     1079    String first = mode.substr( 0, epos ) ;
     1080    string substr[4] ;
     1081    int numSubstr = split( first, substr, 4, sep ) ;
     1082    String m1( substr[0] ) ;
     1083    String m2( substr[2] ) ;
     1084    if ( numSubstr == 4 ) {
     1085      if ( m1.find( "CALIBRATE" ) == 0 ) {
     1086        if ( m2.find( "ON" ) == 0 )
     1087          st = SrcType::PONCAL ;
     1088        else if ( m2.find( "OFF" ) == 0 )
     1089          st = SrcType::POFFCAL ;
     1090      }
     1091      else if ( m1.find( "OBSERVE" ) == 0 ) {
     1092        if ( m2.find( "ON" ) == 0 )
     1093          st = SrcType::PSON ;
     1094        else if ( m2.find( "OFF" ) == 0 )
     1095          st = SrcType::PSOFF ;
     1096      }
     1097    }
     1098    else {
     1099      if ( sig ) st = SrcType::SIG ;
     1100      else if ( ref ) st = SrcType::REF ;
     1101    }
     1102  }
     1103  Block<uInt> getPolNos( Vector<Int> &corr )
     1104  {
     1105    Block<uInt> polnos( npol ) ;
     1106    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     1107      if ( corr[ipol] == Stokes::I || corr[ipol] == Stokes::RR || corr[ipol] == Stokes::XX )
     1108        polnos[ipol] = 0 ;
     1109      else if ( corr[ipol] == Stokes::Q || corr[ipol] == Stokes::LL || corr[ipol] == Stokes::YY )
     1110        polnos[ipol] = 1 ;
     1111      else if ( corr[ipol] == Stokes::U || corr[ipol] == Stokes::RL || corr[ipol] == Stokes::XY )
     1112        polnos[ipol] = 2 ;
     1113      else if ( corr[ipol] == Stokes::V || corr[ipol] == Stokes::LR || corr[ipol] == Stokes::YX )
     1114        polnos[ipol] = 3 ;
     1115    }
     1116    return polnos ;
     1117  }
     1118  String getPolType( Int &corr )
     1119  {
     1120    String poltype = "" ;
     1121    if ( corr == Stokes::I || corr == Stokes::Q || corr == Stokes::U || corr == Stokes::V )
     1122      poltype = "stokes" ;
     1123    else if ( corr == Stokes::XX || corr == Stokes::YY || corr == Stokes::XY || corr == Stokes::YX )
     1124      poltype = "linear" ;
     1125    else if ( corr == Stokes::RR || corr == Stokes::LL || corr == Stokes::RL || corr == Stokes::LR )
     1126      poltype = "circular" ;
     1127    else if ( corr == Stokes::Plinear || corr == Stokes::Pangle )
     1128      poltype = "linpol" ;
     1129    return poltype ;   
     1130  }
     1131  uInt getWeatherId()
     1132  {
     1133    // if only one row, return 0
     1134    if ( weatherTime.nelements() == 1 )
     1135      return 0 ;
     1136
     1137    // @todo At the moment, do binary search every time
     1138    //       if this is bottleneck, frequency of binary search must be reduced
     1139    Double t = currentTime.get( "s" ).getValue() ;
     1140    uInt idx = min( binarySearch( weatherTime, t ), weatherTime.nelements()-1 ) ;
     1141    if ( weatherTime[idx] < t ) {
     1142      if ( idx != weatherTime.nelements()-1 ) {
     1143        if ( weatherTime[idx+1] - t < 0.5 * weatherInterval[idx+1] )
     1144          idx++ ;
     1145      }
     1146    }
     1147    else if ( weatherTime[idx] > t ) {
     1148      if ( idx != 0 ) {
     1149        if ( weatherTime[idx] - t > 0.5 * weatherInterval[idx] )
     1150          idx-- ;
     1151      }
     1152    }
     1153    return idx ;
     1154  }
     1155  void processSysCal( Int &spwId )
     1156  {
     1157    // get feedId from row
     1158    Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
     1159
     1160    uInt nrow = sctab.nrow() ;
     1161    ROScalarColumn<Int> col( sctab, "ANTENNA_ID" ) ;
     1162    Vector<Int> aids = col.getColumn() ;
     1163    col.attach( sctab, "FEED_ID" ) ;
     1164    Vector<Int> fids = col.getColumn() ;
     1165    col.attach( sctab, "SPECTRAL_WINDOW_ID" ) ;
     1166    Vector<Int> sids = col.getColumn() ;
     1167    ROScalarColumn<Double> timeCol( sctab, "TIME" ) ;
     1168    ROScalarColumn<Double> intCol( sctab, "INTERVAL" ) ;
     1169    for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
     1170      if ( aids[irow] == antennaId
     1171           && fids[irow] == feedId
     1172           && sids[irow] == spwId ) {
     1173        syscalRow[numSysCalRow] = irow ;
     1174        syscalTime[numSysCalRow] = timeCol( irow ) ;
     1175        syscalInterval[numSysCalRow] = intCol( irow ) ;
     1176        numSysCalRow++ ;
     1177      }
     1178    }
     1179  }
     1180  uInt getSysCalIndex()
     1181  {
     1182    // if only one row, return 0
     1183    if ( numSysCalRow == 1 || !isSysCal )
     1184      return 0 ;
     1185
     1186    // @todo At the moment, do binary search every time
     1187    //       if this is bottleneck, frequency of binary search must be reduced
     1188    Double t = currentTime.get( "s" ).getValue() ;
     1189    Vector<Double> tslice  = syscalTime( Slice(0, numSysCalRow) ) ;
     1190    uInt idx = min( binarySearch( tslice, t ), numSysCalRow-1 ) ;
     1191    if ( syscalTime[idx] < t ) {
     1192      if ( idx != numSysCalRow-1 ) {
     1193        if ( syscalTime[idx+1] - t < 0.5 * syscalInterval[idx+1] )
     1194          idx++ ;
     1195      }
     1196    }
     1197    else if ( syscalTime[idx] > t ) {
     1198      if ( idx != 0 ) {
     1199        if ( syscalTime[idx] - t > 0.5 * syscalInterval[idx] )
     1200          idx-- ;
     1201      }
     1202    }
     1203    return idx ;   
     1204  }
     1205  Block<uInt> getTcalId( Double &t )
     1206  {
     1207    // return 0 if no SysCal table
     1208    if ( !isSysCal ) {
     1209      return Block<uInt>( 4, 0 ) ;
     1210    }
     1211     
     1212    // get feedId from row
     1213    Int feedId = (Int)tablerow.record().asuInt( "BEAMNO" ) ;
     1214
     1215    // key
     1216    String key = keyTcal( feedId, spwId, t ) ;
     1217   
     1218    // retrieve ids
     1219    Vector<uInt> ids = syscalRecord.asArrayuInt( key ) ;
     1220    //Vector<uInt> ids = syscalRecord[key] ;
     1221    uInt np = ids[1] - ids[0] + 1 ;
     1222    Block<uInt> tcalids( np ) ;
     1223    if ( np > 0 ) {
     1224      tcalids[0] = ids[0] ;
     1225      if ( np > 1 ) {
     1226        tcalids[1] = ids[1] ;
     1227        for ( uInt ip = 2 ; ip < np ; ip++ )
     1228          tcalids[ip] = ids[0] + ip - 1 ;
     1229      }
     1230    }
     1231    return tcalids ;
     1232  }
     1233  uInt maxNumPol()
     1234  {
     1235    ROScalarColumn<Int> numCorrCol( poltab, "NUM_CORR" ) ;
     1236    return max( numCorrCol.getColumn() ) ;
     1237  }
     1238
     1239  Scantable &scantable;
     1240  Int antennaId;
     1241  uInt rowidx;
     1242  String dataColumnName;
     1243  TableRow tablerow;
     1244  STHeader header;
     1245  Vector<Int> feedEntry;
     1246  uInt nbeam;
     1247  Int npol;
     1248  Int nchan;
     1249  Int sourceId;
     1250  Int polId;
     1251  Int spwId;
     1252  uInt cycleNo;
     1253  MDirection sourceDir;
     1254  MPosition antpos;
     1255  MEpoch currentTime;
     1256  map<Int,uInt> ifmap;
     1257  Block<uInt> polnos;
     1258  Bool getpt;
     1259  Vector<Double> pointingTime;
     1260  Cube<Double> pointingDirection;
     1261  MDirection::Types dirType;
     1262  Bool isWeather;
     1263  Vector<Double> weatherTime;
     1264  Vector<Double> weatherInterval;
     1265  Bool isSysCal;
     1266  Record syscalRecord;
     1267  //map< String,Vector<uInt> > syscalRecord;
     1268  uInt numSysCalRow ;
     1269  Vector<uInt> syscalRow;
     1270  Vector<Double> syscalTime;
     1271  Vector<Double> syscalInterval;
     1272  //String tsysCol;
     1273  //String tcalCol;
     1274
     1275  // MS subtables
     1276  Table obstab;
     1277  Table sctab;
     1278  Table spwtab;
     1279  Table statetab;
     1280  Table ddtab;
     1281  Table poltab;
     1282  Table fieldtab;
     1283  Table anttab;
     1284  Table srctab;
     1285  Matrix<Float> sp;
     1286  Matrix<uChar> fl;
     1287  MeasFrame mf;
     1288
     1289  // MS MAIN columns
     1290  ROTableColumn intervalCol;
     1291  ROTableColumn flagRowCol;
     1292  ROArrayColumn<Float> floatDataCol;
     1293  ROArrayColumn<Complex> dataCol;
     1294  ROArrayColumn<Bool> flagCol;
     1295
     1296  // MS SYSCAL columns
     1297  ROArrayColumn<Float> sysCalTsysCol;
     1298
     1299  // Scantable MAIN columns
     1300  RecordFieldPtr<Double> timeRF,intervalRF,sourceVelocityRF;
     1301  RecordFieldPtr< Vector<Double> > directionRF,scanRateRF,
     1302    sourceProperMotionRF,sourceDirectionRF;
     1303  RecordFieldPtr<Float> azimuthRF,elevationRF;
     1304  RecordFieldPtr<uInt> weatherIdRF,cycleNoRF,flagRowRF,polNoRF,tcalIdRF,
     1305    ifNoRF,freqIdRF,moleculeIdRF,beamNoRF,focusIdRF,scanNoRF;
     1306  RecordFieldPtr< Vector<Float> > spectraRF,tsysRF;
     1307  RecordFieldPtr< Vector<uChar> > flagtraRF;
     1308  RecordFieldPtr<String> sourceNameRF,fieldNameRF;
     1309  RecordFieldPtr<Int> sourceTypeRF;
     1310};
     1311
     1312class BaseTcalVisitor: public TableVisitor {
     1313  uInt lastRecordNo ;
     1314  Int lastAntennaId ;
     1315  Int lastFeedId ;
     1316  Int lastSpwId ;
     1317  Double lastTime ;
     1318protected:
     1319  const Table &table;
     1320  uInt count;
     1321public:
     1322  BaseTcalVisitor(const Table &table)
     1323   : table(table)
     1324  {
     1325    count = 0;
     1326  }
     1327 
     1328  virtual void enterAntennaId(const uInt recordNo, Int columnValue) { }
     1329  virtual void leaveAntennaId(const uInt recordNo, Int columnValue) { }
     1330  virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
     1331  virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
     1332  virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
     1333  virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
     1334  virtual void enterTime(const uInt recordNo, Double columnValue) { }
     1335  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     1336
     1337  virtual Bool visitRecord(const uInt recordNo,
     1338                           const Int antennaId,
     1339                           const Int feedId,
     1340                           const Int spwId,
     1341                           const Double time) { return True ; }
     1342
     1343  virtual Bool visit(Bool isFirst, const uInt recordNo,
     1344                     const uInt nCols, void const *const colValues[]) {
     1345    Int antennaId, feedId, spwId;
     1346    Double time;
     1347    { // prologue
     1348      uInt i = 0;
     1349      {
     1350        const Int *col = (const Int *)colValues[i++];
     1351        antennaId = col[recordNo];
     1352      }
     1353      {
     1354        const Int *col = (const Int *)colValues[i++];
     1355        feedId = col[recordNo];
     1356      }
     1357      {
     1358        const Int *col = (const Int *)colValues[i++];
     1359        spwId = col[recordNo];
     1360      }
     1361      {
     1362        const Double *col = (const Double *)colValues[i++];
     1363        time = col[recordNo];
     1364      }
     1365      assert(nCols == i);
     1366    }
     1367
     1368    if (isFirst) {
     1369      enterAntennaId(recordNo, antennaId);
     1370      enterFeedId(recordNo, feedId);
     1371      enterSpwId(recordNo, spwId);
     1372      enterTime(recordNo, time);
     1373    } else {
     1374      if ( lastAntennaId != antennaId ) {
     1375        leaveTime(lastRecordNo, lastTime);
     1376        leaveSpwId(lastRecordNo, lastSpwId);
     1377        leaveFeedId(lastRecordNo, lastFeedId);
     1378        leaveAntennaId(lastRecordNo, lastAntennaId);
     1379
     1380        enterAntennaId(recordNo, antennaId);
     1381        enterFeedId(recordNo, feedId);
     1382        enterSpwId(recordNo, spwId);
     1383        enterTime(recordNo, time);
     1384      }       
     1385      else if (lastFeedId != feedId) {
     1386        leaveTime(lastRecordNo, lastTime);
     1387        leaveSpwId(lastRecordNo, lastSpwId);
     1388        leaveFeedId(lastRecordNo, lastFeedId);
     1389
     1390        enterFeedId(recordNo, feedId);
     1391        enterSpwId(recordNo, spwId);
     1392        enterTime(recordNo, time);
     1393      } else if (lastSpwId != spwId) {
     1394        leaveTime(lastRecordNo, lastTime);
     1395        leaveSpwId(lastRecordNo, lastSpwId);
     1396
     1397        enterSpwId(recordNo, spwId);
     1398        enterTime(recordNo, time);
     1399      } else if (lastTime != time) {
     1400        leaveTime(lastRecordNo, lastTime);
     1401        enterTime(recordNo, time);
     1402      }
     1403    }
     1404    count++;
     1405    Bool result = visitRecord(recordNo, antennaId, feedId, spwId, time);
     1406
     1407    { // epilogue
     1408      lastRecordNo = recordNo;
     1409
     1410      lastAntennaId = antennaId;
     1411      lastFeedId = feedId;
     1412      lastSpwId = spwId;
     1413      lastTime = time;
     1414    }
     1415    return result ;
     1416  }
     1417
     1418  virtual void finish() {
     1419    if (count > 0) {
     1420      leaveTime(lastRecordNo, lastTime);
     1421      leaveSpwId(lastRecordNo, lastSpwId);
     1422      leaveFeedId(lastRecordNo, lastFeedId);
     1423      leaveAntennaId(lastRecordNo, lastAntennaId);
     1424    }
     1425  }
     1426};
     1427
     1428class TcalVisitor: public BaseTcalVisitor, public MSFillerUtils {
     1429public:
     1430  TcalVisitor(const Table &table, Table &tcaltab, Record &r, Int aid )
     1431  //TcalVisitor(const Table &table, Table &tcaltab, map< String,Vector<uInt> > &r, Int aid )
     1432    : BaseTcalVisitor( table ),
     1433      tcal(tcaltab),
     1434      rec(r),
     1435      antenna(aid)
     1436  {
     1437    process = False ;
     1438    rowidx = 0 ;
     1439
     1440    // attach to SYSCAL columns
     1441    timeCol.attach( table, "TIME" ) ;
     1442
     1443    // add rows
     1444    uInt addrow = table.nrow() * 4 ;
     1445    tcal.addRow( addrow ) ;
     1446
     1447    // attach to TCAL columns
     1448    row = TableRow( tcal ) ;
     1449    TableRecord &trec = row.record() ;
     1450    idRF.attachToRecord( trec, "ID" ) ;
     1451    timeRF.attachToRecord( trec, "TIME" ) ;
     1452    tcalRF.attachToRecord( trec, "TCAL" ) ;
     1453  }
     1454
     1455  virtual void enterAntennaId(const uInt recordNo, Int columnValue) {
     1456    if ( columnValue == antenna )
     1457      process = True ;
     1458  }
     1459  virtual void leaveAntennaId(const uInt recordNo, Int columnValue) {
     1460    process = False ;
     1461  }
     1462  virtual void enterFeedId(const uInt recordNo, Int columnValue) { }
     1463  virtual void leaveFeedId(const uInt recordNo, Int columnValue) { }
     1464  virtual void enterSpwId(const uInt recordNo, Int columnValue) { }
     1465  virtual void leaveSpwId(const uInt recordNo, Int columnValue) { }
     1466  virtual void enterTime(const uInt recordNo, Double columnValue) {
     1467    qtime = timeCol( recordNo ) ;
     1468  }
     1469  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     1470  virtual Bool visitRecord(const uInt recordNo,
     1471                           const Int antennaId,
     1472                           const Int feedId,
     1473                           const Int spwId,
     1474                           const Double time)
     1475  {
     1476    //cout << "(" << recordNo << "," << antennaId << "," << feedId << "," << spwId << ")" << endl ;
     1477    if ( process ) {
     1478      String sTime = MVTime( qtime ).string( MVTime::YMD ) ;
     1479      *timeRF = sTime ;
     1480      uInt oldidx = rowidx ;
     1481      Matrix<Float> subtcal = tcalCol( recordNo ) ;
     1482      Vector<uInt> idminmax( 2 ) ;
     1483      for ( uInt ipol = 0 ; ipol < subtcal.nrow() ; ipol++ ) {
     1484        *idRF = rowidx ;
     1485        tcalRF.define( subtcal.row( ipol ) ) ;
     1486       
     1487        // commit row
     1488        row.put( rowidx ) ;
     1489        rowidx++ ;
     1490      }
     1491     
     1492      idminmax[0] = oldidx ;
     1493      idminmax[1] = rowidx - 1 ;
     1494     
     1495      String key = keyTcal( feedId, spwId, sTime ) ;
     1496      rec.define( key, idminmax ) ;
     1497      //rec[key] = idminmax ;
     1498    }
     1499    return True ;
     1500  }
     1501  virtual void finish()
     1502  {
     1503    BaseTcalVisitor::finish() ;
     1504
     1505    if ( tcal.nrow() > rowidx ) {
     1506      uInt numRemove = tcal.nrow() - rowidx ;
     1507      //cout << "numRemove = " << numRemove << endl ;
     1508      Vector<uInt> rows( numRemove ) ;
     1509      indgen( rows, rowidx ) ;
     1510      tcal.removeRow( rows ) ;
     1511    }
     1512
     1513  }
     1514  void setTcalColumn( String &col )
     1515  {
     1516    //colName = col ;
     1517    tcalCol.attach( table, col ) ;
     1518  }
     1519private:
     1520  Table &tcal;
     1521  Record &rec;
     1522  //map< String,Vector<uInt> > &rec;
     1523  Int antenna;
     1524  uInt rowidx;
     1525  Bool process;
     1526  Quantum<Double> qtime;
     1527  TableRow row;
     1528  String colName;
     1529
     1530  // MS SYSCAL columns
     1531  ROScalarQuantColumn<Double> timeCol;
     1532  ROArrayColumn<Float> tcalCol;
     1533
     1534  // TCAL columns
     1535  RecordFieldPtr<uInt> idRF;
     1536  RecordFieldPtr<String> timeRF;
     1537  RecordFieldPtr< Vector<Float> > tcalRF;
     1538};
    681539
    691540MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
     
    1221593
    1231594  MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
    124   //mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
    125   //                     && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;
    1261595  tablename_ = tmpMS->tableName() ;
    1271596  if ( antenna_ == -1 && antennaStr_.size() > 0 ) {
     
    1391608  mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
    1401609                                       && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
    141 //   stringstream ss ;
    142 //   ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 == " << antenna_ ;
    143 //   String taql( ss.str() ) ;
    144 //   mstable_ = MeasurementSet( tableCommand( taql, *tmpMS ) ) ;
     1610
    1451611  delete tmpMS ;
    1461612
     
    1561622void MSFiller::fill()
    1571623{
    158   os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
    1591624  //double startSec = mathutil::gettimeofday_sec() ;
    1601625  //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
    1611626
    162   //double time0 = mathutil::gettimeofday_sec() ;
    163   //os_ << "start init fill: " << time0 << LogIO::POST ;
     1627  os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
    1641628
    1651629  // Initialize header
    1661630  STHeader sdh ; 
    1671631  initHeader( sdh ) ;
     1632  table_->setHeader( sdh ) ;
    1681633 
    1691634  // check if optional table exists
    170   //const TableRecord msrec = tablesel_.keywordSet() ;
    171   const TableRecord msrec = mstable_.keywordSet() ;
     1635  const TableRecord &msrec = mstable_.keywordSet() ;
    1721636  isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
    1731637  if ( isDoppler_ )
     
    1991663      isWeather_ = False ;
    2001664
    201   // Access to MS subtables
    202   MSField fieldtab = mstable_.field() ;
    203   MSPolarization poltab = mstable_.polarization() ;
    204   MSDataDescription ddtab = mstable_.dataDescription() ;
    205   MSObservation obstab = mstable_.observation() ;
    206   MSSource srctab = mstable_.source() ;
    207   MSSpectralWindow spwtab = mstable_.spectralWindow() ;
    208   MSSysCal caltab = mstable_.sysCal() ;
    209   if ( caltab.nrow() == 0 )
    210     isSysCal_ = False ;
    211   else {
     1665  // column name for Tsys and Tcal
     1666  if ( isSysCal_ ) {
     1667    const MSSysCal &caltab = mstable_.sysCal() ;
    2121668    if ( !caltab.tableDesc().isColumn( colTcal_ ) ) {
    2131669      colTcal_ = "TCAL" ;
     
    2211677    }
    2221678  }
    223 //   colTcal_ = "TCAL" ;
    224 //   colTsys_ = "TSYS" ;
    225   MSPointing pointtab = mstable_.pointing() ;
    226   if ( mstable_.weather().nrow() == 0 )
    227     isWeather_ = False ;
    228   MSState stattab = mstable_.state() ;
    229   MSAntenna anttab = mstable_.antenna() ;
    230 
    231   // TEST
    232   // memory allocation by boost::object_pool
    233   boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ;
    234   //
     1679  else {
     1680    colTcal_ = "NONE" ;
     1681    colTsys_ = "NONE" ;
     1682  }
     1683
     1684  // Access to MS subtables
     1685  //MSField &fieldtab = mstable_.field() ;
     1686  //MSPolarization &poltab = mstable_.polarization() ;
     1687  //MSDataDescription &ddtab = mstable_.dataDescription() ;
     1688  //MSObservation &obstab = mstable_.observation() ;
     1689  //MSSource &srctab = mstable_.source() ;
     1690  //MSSpectralWindow &spwtab = mstable_.spectralWindow() ;
     1691  //MSSysCal &caltab = mstable_.sysCal() ;
     1692  MSPointing &pointtab = mstable_.pointing() ;
     1693  //MSState &stattab = mstable_.state() ;
     1694  //MSAntenna &anttab = mstable_.antenna() ;
    2351695
    2361696  // SUBTABLES: FREQUENCIES
     
    2471707
    2481708  // SUBTABLES: TCAL
    249   fillTcal( tpoolr ) ;
     1709  fillTcal() ;
    2501710
    2511711  // SUBTABLES: FIT
     
    2551715  //fillHistory() ;
    2561716
    257   // MAIN
    258   // Iterate over several ids
    259   map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair
    260   ROArrayQuantColumn<Double> *sharedQDArrCol = new ROArrayQuantColumn<Double>( anttab, "POSITION" ) ;
    261   Vector< Quantum<Double> > antpos = (*sharedQDArrCol)( antenna_ ) ;
    262   delete sharedQDArrCol ;
    263   MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
    264   Vector<Double> pt ;
    265   ROArrayColumn<Double> pdcol ;
    266   Vector<Double> defaultScanrate( 2, 0.0 ) ;
    267   if ( getPt_ ) {
    268     pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
    269     ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ;
    270     ptcol.getColumn( pt ) ;
    271     TableRecord trec = ptcol.keywordSet() ;
    272     String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ;
    273     if ( tUnit == "d" )
    274       pt *= 86400.0 ;
    275     pdcol.attach( pointtab, "DIRECTION" ) ;
    276   }
    277   String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ;
    278   String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ;
    279   sdh.antennaposition.resize( 3 ) ;
    280   for ( int i = 0 ; i < 3 ; i++ )
    281     sdh.antennaposition[i] = antpos[i].getValue( "m" ) ;
    282   String telescopeName = "" ;
    283 
    284   //double time1 = mathutil::gettimeofday_sec() ;
    285   //os_ << "end init fill: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    286 
    287   // row based
    288   Table &stab = table_->table() ;
    289   TableRow row( stab ) ;
    290   TableRecord &trec = row.record() ;
    291   RecordFieldPtr< Array<Float> > spRF( trec, "SPECTRA" ) ;
    292   RecordFieldPtr< Array<uChar> > ucarrRF( trec, "FLAGTRA" ) ;
    293   RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
    294   RecordFieldPtr< Array<Float> > tsysRF( trec, "TSYS" ) ;
    295   RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
    296   RecordFieldPtr< Array<Double> > dirRF( trec, "DIRECTION" ) ;
    297   RecordFieldPtr<Float> azRF( trec, "AZIMUTH" ) ;
    298   RecordFieldPtr<Float> elRF( trec, "ELEVATION" ) ;
    299   RecordFieldPtr< Array<Double> > scrRF( trec, "SCANRATE" ) ;
    300   RecordFieldPtr<uInt> cycleRF( trec, "CYCLENO" ) ;
    301   RecordFieldPtr<uInt> flrRF( trec, "FLAGROW" ) ;
    302   RecordFieldPtr<uInt> tcalidRF( trec, "TCAL_ID" ) ;
    303   RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ;
    304   RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ;
    305   RecordFieldPtr<Int> refbRF( trec, "REFBEAMNO" ) ;
    306   RecordFieldPtr<Int> fitidRF( trec, "FIT_ID" ) ;
    307   RecordFieldPtr<Float> tauRF( trec, "OPACITY" ) ;
    308   RecordFieldPtr<uInt> beamRF( trec, "BEAMNO" ) ;
    309   RecordFieldPtr<uInt> focusidRF( trec, "FOCUS_ID" ) ;
    310   RecordFieldPtr<uInt> ifnoRF( trec, "IFNO" ) ;
    311   RecordFieldPtr<uInt> molidRF( trec, "MOLECULE_ID" ) ;
    312   RecordFieldPtr<uInt> freqidRF( trec, "FREQ_ID" ) ;
    313   RecordFieldPtr<uInt> scanRF( trec, "SCANNO" ) ;
    314   RecordFieldPtr<Int> srctypeRF( trec, "SRCTYPE" ) ;
    315   RecordFieldPtr<String> srcnameRF( trec, "SRCNAME" ) ;
    316   RecordFieldPtr<String> fieldnameRF( trec, "FIELDNAME" ) ;
    317   RecordFieldPtr< Array<Double> > srcpmRF( trec, "SRCPROPERMOTION" ) ;
    318   RecordFieldPtr< Array<Double> > srcdirRF( trec, "SRCDIRECTION" ) ;
    319   RecordFieldPtr<Double> sysvelRF( trec, "SRCVELOCITY" ) ;
    320 
    321   // REFBEAMNO
    322   *refbRF = -1 ;
    323 
    324   // FIT_ID
    325   *fitidRF = -1 ;
    326 
    327   // OPACITY
    328   *tauRF = 0.0 ;
    329 
    330   //
    331   // ITERATION: OBSERVATION_ID
    332   //
    333   TableIterator iter0( mstable_, "OBSERVATION_ID" ) ;
    334   while( !iter0.pastEnd() ) {
    335     //time0 = mathutil::gettimeofday_sec() ;
    336     //os_ << "start 0th iteration: " << time0 << LogIO::POST ;
    337     Table t0 = iter0.table() ;
    338     Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ;
    339     if ( sdh.observer == "" ) {
    340       sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ;
    341     }
    342     if ( sdh.project == "" ) {
    343       sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ;
    344     }
    345     ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ;
    346     MEpoch me = (*tmpMeasCol)( obsId )( IPosition(1,0) ) ;
    347     delete tmpMeasCol ;
    348     if ( sdh.utc == 0.0 ) {
    349       sdh.utc = me.get( "d" ).getValue() ;
    350     }
    351     if ( telescopeName == "" ) {
    352       telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ;
    353     }
    354     Int nbeam = 0 ;
    355     //time1 = mathutil::gettimeofday_sec() ;
    356     //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    357     //
    358     // ITERATION: FEED1
    359     //
    360     TableIterator iter1( t0, "FEED1" ) ;
    361     while( !iter1.pastEnd() ) {
    362       //time0 = mathutil::gettimeofday_sec() ;
    363       //os_ << "start 1st iteration: " << time0 << LogIO::POST ;
    364       Table t1 = iter1.table() ;
    365       // assume FEED1 == FEED2
    366       Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ;
    367       nbeam++ ;
    368 
    369       // BEAMNO
    370       *beamRF = feedId ;
    371 
    372       // FOCUS_ID
    373       *focusidRF = 0 ;
    374 
    375       //time1 = mathutil::gettimeofday_sec() ;
    376       //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    377       //
    378       // ITERATION: FIELD_ID
    379       //
    380       TableIterator iter2( t1, "FIELD_ID" ) ;
    381       while( !iter2.pastEnd() ) {
    382         //time0 = mathutil::gettimeofday_sec() ;
    383         //os_ << "start 2nd iteration: " << time0 << LogIO::POST ;
    384         Table t2 = iter2.table() ;
    385         Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ;
    386         Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ;
    387         String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ;
    388         fieldName += "__" + String::toString(fieldId) ;
    389         ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ;
    390         Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ;
    391         delete delayDirCol ;         
    392 
    393         // FIELDNAME
    394         *fieldnameRF = fieldName ;
    395 
    396 
    397         //time1 = mathutil::gettimeofday_sec() ;
    398         //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    399         //
    400         // ITERATION: DATA_DESC_ID
    401         //
    402         TableIterator iter3( t2, "DATA_DESC_ID" ) ;
    403         while( !iter3.pastEnd() ) {
    404           //time0 = mathutil::gettimeofday_sec() ;
    405           //os_ << "start 3rd iteration: " << time0 << LogIO::POST ;
    406           Table t3 = iter3.table() ;
    407           Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ;
    408           Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ;
    409           Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ;
    410 
    411           // IFNO
    412           *ifnoRF = (uInt)spwId ;
    413 
    414           // polarization information
    415           Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ;
    416           ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ;
    417           Vector<Int> corrtype = (*roArrICol)( polId ) ;
    418           delete roArrICol ;
    419           String srcName( "" ) ;
    420           Vector<Double> srcPM( 2, 0.0 ) ;
    421           Vector<Double> srcDir( 2, 0.0 ) ;
    422           MDirection md ;
    423           Vector<Double> restFreqs ;
    424           Vector<String> transitionName ;
    425           Vector<Double> sysVels ;
    426 //           os_ << "npol = " << npol << LogIO::POST ;
    427 //           os_ << "corrtype = " << corrtype << LogIO::POST ;
    428 
    429           // source and molecular transition
    430           sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ;
    431 //           os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;
    432 
    433           // SRCNAME
    434           *srcnameRF = srcName ;
    435 
    436 //           os_ << "srcName = " << srcName << LogIO::POST ;
    437 
    438           // SRCPROPERMOTION
    439           *srcpmRF = srcPM ;
    440 
    441           //os_ << "srcPM = " << srcPM << LogIO::POST ;
    442 
    443           // SRCDIRECTION
    444           *srcdirRF = md.getAngle().getValue( "rad" ) ;
    445 
    446           //os_ << "srcDir = " << srcDir << LogIO::POST ;
    447 
    448           // SRCVELOCITY
    449           Double sysVel = 0.0 ;
    450           if ( !sysVels.empty() )
    451             sysVel = sysVels[0] ;
    452           *sysvelRF = sysVel ;
    453 
    454 //           os_ << "sysVel = " << sysVel << LogIO::POST ;
    455 
    456           uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ;
    457 
    458           // MOLECULE_ID
    459           *molidRF = molId ;
    460 
    461           // spectral setup
    462           uInt freqId ;
    463           Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ;
    464           Bool iswvr = False ;
    465           if ( nchan == 4 ) iswvr = True ;
    466           sdh.nchan = max( sdh.nchan, nchan ) ;
    467           map<Int,uInt>::iterator iter = ifmap.find( spwId ) ;
    468           if ( iter == ifmap.end() ) {
    469             ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ;
    470             me = (*tmpMeasCol)( 0 ) ;
    471             delete tmpMeasCol ;
    472             Double refpix ;
    473             Double refval ;
    474             Double increment ;
    475             spectralSetup( spwId,
    476                            me,
    477                            mp,
    478                            md,
    479                            refpix,
    480                            refval,
    481                            increment,
    482                            nchan,
    483                            sdh.freqref,
    484                            sdh.reffreq,
    485                            sdh.bandwidth,
    486                            tpoolr ) ;
    487             freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
    488             ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
    489             //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
    490           }
    491           else {
    492             freqId = iter->second ;
    493           }
    494 
    495           // FREQ_ID
    496           *freqidRF = freqId ;
    497 
    498           // for TSYS and TCAL
    499           Vector<MEpoch> scTime ;
    500           Vector<Double> scInterval ;
    501           ROArrayColumn<Float> scTsysCol ;
    502           MSSysCal caltabsel ;
    503           if ( isSysCal_ ) {
    504             caltabsel = caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ;
    505             ROScalarMeasColumn<MEpoch> scTimeCol( caltabsel, "TIME" ) ;
    506             scTime.resize( caltabsel.nrow() ) ;
    507             for ( uInt irow = 0 ; irow < caltabsel.nrow() ; irow++ )
    508               scTime[irow] = scTimeCol( irow ) ;
    509             ROScalarColumn<Double> scIntervalCol( caltabsel, "INTERVAL" ) ;
    510             scIntervalCol.getColumn( scInterval ) ;
    511             if ( colTsys_ != "NONE" )
    512               scTsysCol.attach( caltabsel, colTsys_ ) ;
    513           }
    514 
    515           sdh.npol = max( sdh.npol, npol ) ;
    516           if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
    517 
    518           //time1 = mathutil::gettimeofday_sec() ;
    519           //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    520           //
    521           // ITERATION: SCAN_NUMBER
    522           //
    523           TableIterator iter4( t3, "SCAN_NUMBER" ) ;
    524           while( !iter4.pastEnd() ) {
    525             //time0 = mathutil::gettimeofday_sec() ;
    526             //os_ << "start 4th iteration: " << time0 << LogIO::POST ;
    527             Table t4 = iter4.table() ;
    528             Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ;
    529 
    530             // SCANNO
    531             *scanRF = scanNum - 1 ;
    532 
    533             uInt cycle = 0 ;
    534 
    535             //time1 = mathutil::gettimeofday_sec() ;
    536             //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    537             //
    538             // ITERATION: STATE_ID
    539             //
    540             TableIterator iter5( t4, "STATE_ID" ) ;
    541             while( !iter5.pastEnd() ) {
    542               //time0 = mathutil::gettimeofday_sec() ;
    543               //os_ << "start 5th iteration: " << time0 << LogIO::POST ;
    544               Table t5 = iter5.table() ;
    545               Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ;
    546               String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ;
    547               if ( sdh.obstype == "" ) sdh.obstype = obstype ;
    548 
    549               Int nrow = t5.nrow() ;
    550               //time1 = mathutil::gettimeofday_sec() ;
    551               //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    552 
    553               Cube<Float> spArr ;
    554               Cube<Bool> flArr ;
    555               reshapeSpectraAndFlagtra( spArr,
    556                                         flArr,
    557                                         t5,
    558                                         npol,
    559                                         nchan,
    560                                         nrow,
    561                                         corrtype ) ;
    562               if ( sdh.fluxunit == "" ) {
    563                 String colName = "FLOAT_DATA" ;
    564                 if ( isData_ ) colName = "DATA" ;
    565                 ROTableColumn dataCol( t5, colName ) ;
    566                 const TableRecord &dataColKeys = dataCol.keywordSet() ;
    567                 if ( dataColKeys.isDefined( "UNIT" ) )
    568                   sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
    569                 else if ( dataColKeys.isDefined( "QuantumUnits" ) )
    570                   sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ;
    571               }
    572 
    573               ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
    574               Block<MEpoch> mTimeB( nrow ) ;
    575               for ( Int irow = 0 ; irow < nrow ; irow++ )
    576                 mTimeB[irow] = (*mTimeCol)( irow ) ;
    577               Block<Int> sysCalIdx( nrow, -1 ) ;
    578               if ( isSysCal_ ) {
    579                 getSysCalTime( scTime, scInterval, mTimeB, sysCalIdx ) ;
    580               }
    581               delete mTimeCol ;
    582               Matrix<Float> defaulttsys( npol, 1, 1.0 ) ;
    583               Int srcType = getSrcType( stateId, tpoolr ) ;
    584               uInt diridx = 0 ;
    585               uInt wid = 0 ;
    586               Int pidx = 0 ;
    587               Bool crossOK = False ;
    588               Block<uInt> polnos( npol, 99 ) ;
    589               for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
    590                 Block<uInt> p = getPolNo( corrtype[ipol] ) ;
    591                 if ( p.size() > 1 ) {
    592                   if ( crossOK ) continue ;
    593                   else {
    594                     polnos[pidx] = p[0] ;
    595                     pidx++ ;
    596                     polnos[pidx] = p[1] ;
    597                     pidx++ ;
    598                     crossOK = True ;
    599                   }
    600                 }
    601                 else {
    602                   polnos[pidx] = p[0] ;
    603                   pidx++ ;
    604                 }
    605               }
    606              
    607               // SRCTYPE
    608               *srctypeRF = srcType ;
    609 
    610               for ( Int irow = 0 ; irow < nrow ; irow++ ) {
    611                 // CYCLENO
    612                 *cycleRF = cycle ;
    613 
    614                 // FLAGROW
    615                 *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ;
    616 
    617                 // SPECTRA, FLAG
    618                 Matrix<Float> sp = spArr.xyPlane( irow ) ;
    619                 Matrix<Bool> flb = flArr.xyPlane( irow ) ;
    620                 Matrix<uChar> fl( flb.shape() ) ;
    621                 convertArray( fl, flb ) ;
    622 
    623                 // TIME
    624                 *timeRF = mTimeB[irow].get("d").getValue() ;
    625 
    626                 // INTERVAL
    627                 *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ;
    628 
    629                 // TSYS
    630                 Matrix<Float> tsys ;
    631                 if ( sysCalIdx[irow] != -1 && colTsys_ != "NONE" )
    632                   tsys = scTsysCol( sysCalIdx[irow] ) ;
    633                 else
    634                   tsys = defaulttsys ;
    635 
    636                 // TCAL_ID
    637                 Block<uInt> tcalids( npol, 0 ) ;
    638                 if ( sysCalIdx[irow] != -1 && colTcal_ != "NONE" ) {
    639                   tcalids = getTcalId( feedId, spwId, scTime[sysCalIdx[irow]] ) ;
    640                 }
    641 
    642                 // WEATHER_ID
    643                 if ( isWeather_ ) {
    644                   wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ;
    645                   *widRF = mwIndex_[wid] ;
    646                 }
    647                 else {
    648                   *widRF = wid ;
    649                 }
    650                  
    651 
    652                 // DIRECTION, AZEL, SCANRATE
    653                 Vector<Double> dir ;
    654                 Vector<Double> azel ;
    655                 Vector<Double> scanrate = defaultScanrate ;
    656                 String refString ;
    657                 if ( getPt_ )
    658                   diridx = getDirection( diridx, dir, azel, scanrate, pt, pdcol, mTimeB[irow], mp ) ;
    659                 else
    660                   getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ;
    661                 *dirRF = dir ;
    662                 *azRF = azel[0] ;
    663                 *elRF = azel[1] ;
    664                 *scrRF = scanrate ;
    665 
    666                 // Polarization dependent things
    667                 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
    668                   // POLNO
    669                   *polnoRF = polnos[ipol] ;
    670 
    671                   spRF.define( sp.row( ipol ) ) ;
    672                   ucarrRF.define( fl.row( ipol ) ) ;
    673                   tsysRF.define( tsys.row( ipol ) ) ;
    674                   *tcalidRF = tcalids[ipol] ;
    675 
    676                   // Commit row
    677                   stab.addRow() ;
    678                   row.put( stab.nrow()-1 ) ;
    679                 }
    680 
    681                 cycle++ ;
    682               }
    683              
    684               //time1 = mathutil::gettimeofday_sec() ;
    685               //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    686 
    687               iter5.next() ;
    688             }
    689             iter4.next() ;
    690           }
    691           iter3.next() ;
    692         }
    693         iter2.next() ;
    694       }
    695       iter1.next() ;
    696     }
    697     if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
    698 
    699     iter0.next() ;
    700   }
    701 
    702 
    703   delete tpoolr ;
    704 
    705 
    706   // Table Keywords
    707   sdh.nif = ifmap.size() ;
    708   if ( ( telescopeName == "" ) || ( antennaName == telescopeName ) ) {
    709     sdh.antennaname = antennaName ;
    710   }
    711   else {
    712     sdh.antennaname = telescopeName + "//" + antennaName ;
    713   }
    714   if ( stationName != "" && stationName != antennaName ) {
    715     sdh.antennaname += "@" + stationName ;
    716   }
    717   ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ;
    718   String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ;
    719   if ( dirref == "AZELGEO" || dirref == "AZEL" ) {
    720     dirref = "J2000" ;
    721   }
    722   sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
    723   sdh.epoch = "UTC" ;
    724   if (sdh.freqref == "TOPO") {
    725     sdh.freqref = "TOPOCENT";
    726   } else if (sdh.freqref == "GEO") {
    727     sdh.freqref = "GEOCENTR";
    728   } else if (sdh.freqref == "BARY") {
    729     sdh.freqref = "BARYCENT";
    730   } else if (sdh.freqref == "GALACTO") {
    731     sdh.freqref = "GALACTOC";
    732   } else if (sdh.freqref == "LGROUP") {
    733     sdh.freqref = "LOCALGRP";
    734   } else if (sdh.freqref == "CMB") {
    735     sdh.freqref = "CMBDIPOL";
    736   } else if (sdh.freqref == "REST") {
    737     sdh.freqref = "SOURCE";
    738   }
    739 
    740   if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" )
    741     sdh.fluxunit = "K" ;
    742   table_->setHeader( sdh ) ;
     1717  /***
     1718   * Start iteration using TableVisitor
     1719   ***/
     1720  Table stab = table_->table() ;
     1721  {
     1722    static const char *cols[] = {
     1723      "OBSERVATION_ID", "FEED1", "FIELD_ID", "DATA_DESC_ID", "SCAN_NUMBER",
     1724      "STATE_ID", "TIME",
     1725      NULL
     1726    };
     1727    static const TypeManagerImpl<Int> tmInt;
     1728    static const TypeManagerImpl<Double> tmDouble;
     1729    static const TypeManager *const tms[] = {
     1730      &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmInt, &tmDouble, NULL
     1731    };
     1732    //double t0 = mathutil::gettimeofday_sec() ;
     1733    MSFillerVisitor myVisitor(mstable_, *table_ );
     1734    //double t1 = mathutil::gettimeofday_sec() ;
     1735    //cout << "MSFillerVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
     1736    myVisitor.setAntenna( antenna_ ) ;
     1737    //myVisitor.setHeader( sdh ) ;
     1738    if ( getPt_ ) {
     1739      Table ptsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort( "TIME" ) ;
     1740      myVisitor.setPointingTable( ptsel ) ;
     1741    }
     1742    if ( isWeather_ )
     1743      myVisitor.setWeatherTime( mwTime_, mwInterval_ ) ;
     1744    if ( isSysCal_ )
     1745      myVisitor.setSysCalRecord( tcalrec_ ) ;
     1746   
     1747    //double t2 = mathutil::gettimeofday_sec() ;
     1748    traverseTable(mstable_, cols, tms, &myVisitor);
     1749    //double t3 = mathutil::gettimeofday_sec() ;
     1750    //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
     1751   
     1752    sdh = myVisitor.getHeader() ;
     1753  }
     1754  /***
     1755   * End iteration using TableVisitor
     1756   ***/
     1757
     1758  // set header
     1759  //sdh = myVisitor.getHeader() ;
     1760  //table_->setHeader( sdh ) ;
    7431761
    7441762  // save path to POINTING table
     
    7531771
    7541772  // for GBT
    755   if ( antennaName.contains( "GBT" ) ) {
     1773  if ( sdh.antennaname.contains( "GBT" ) ) {
    7561774    String goTabName = datapath.absoluteName() + "/GBT_GO" ;
    7571775    stab.rwKeywordSet().define( "GBT_GO", goTabName ) ;
     
    7591777
    7601778  // for MS created from ASDM
    761   //mstable_.keywordSet().print(cout) ;
    7621779  const TableRecord &msKeys = mstable_.keywordSet() ;
    7631780  uInt nfields = msKeys.nfields() ;
     
    7821799  //tablesel_.unlock() ;
    7831800  mstable_.unlock() ;
    784 }
    785 
    786 Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
    787 {
    788   //double startSec = mathutil::gettimeofday_sec() ;
    789   //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
    790 
    791   MSState statetab = mstable_.state() ;
    792   String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ;
    793   Bool sig = asBool( "SIG", stateId, statetab, tpool ) ;
    794   Bool ref = asBool( "REF", stateId, statetab, tpool ) ;
    795   Double cal = asDouble( "CAL", stateId, statetab, tpool ) ;
    796   //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
    797 
    798   // determine separator
    799   String sep = "" ;
    800   String tmpStr = obsMode.substr( 0, obsMode.find_first_of( "," ) ) ;
    801   //os_ << "tmpStr = " << tmpStr << LogIO::POST ;
    802   //if ( obsMode.find( ":" ) != String::npos ) {
    803   if ( tmpStr.find( ":" ) != String::npos ) {
    804     sep = ":" ;
    805   }
    806   //else if ( obsMode.find( "." ) != String::npos ) {
    807   else if ( tmpStr.find( "." ) != String::npos ) {
    808     sep = "." ;
    809   }
    810   else if ( tmpStr.find( "#" ) != String::npos ) {
    811     sep = "#" ;
    812   }
    813   //else if ( obsMode.find( "_" ) != String::npos ) {
    814   else if ( tmpStr.find( "_" ) != String::npos ) {
    815     sep = "_" ;
    816   }
    817   //os_ << "separator = " << sep << LogIO::POST ;
    818 
    819   // determine SRCTYPE
    820   Int srcType = SrcType::NOTYPE ;
    821   if ( sep == ":" ) {
    822     // sep == ":"
    823     //
    824     // GBT case
    825     //
    826     // obsMode1=Nod
    827     //    NOD
    828     // obsMode1=OffOn
    829     //    obsMode2=PSWITCHON:  PSON
    830     //    obsMode2=PSWITCHOFF: PSOFF
    831     // obsMode1=??
    832     //    obsMode2=FSWITCH:
    833     //       SIG=1: FSON
    834     //       REF=1: FSOFF
    835     // Calibration scan if CAL != 0
    836     Int epos = obsMode.find_first_of( sep ) ;
    837     Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
    838     String obsMode1 = obsMode.substr( 0, epos ) ;
    839     String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
    840     if ( obsMode1 == "Nod" ) {
    841       srcType = SrcType::NOD ;
    842     }
    843     else if ( obsMode1 == "OffOn" ) {
    844       if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
    845       if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
    846     }
    847     else {
    848       if ( obsMode2 == "FSWITCH" ) {
    849         if ( sig ) srcType = SrcType::FSON ;
    850         if ( ref ) srcType = SrcType::FSOFF ;
    851       }
    852     }
    853     if ( cal > 0.0 ) {
    854       if ( srcType == SrcType::NOD )
    855         srcType = SrcType::NODCAL ;
    856       else if ( srcType == SrcType::PSON )
    857         srcType = SrcType::PONCAL ;
    858       else if ( srcType == SrcType::PSOFF )
    859         srcType = SrcType::POFFCAL ;
    860       else if ( srcType == SrcType::FSON )
    861         srcType = SrcType::FONCAL ;
    862       else if ( srcType == SrcType::FSOFF )
    863         srcType = SrcType::FOFFCAL ;
    864       else
    865         srcType = SrcType::CAL ;
    866     }
    867   }
    868   else if ( sep == "." || sep == "#" ) {
    869     // sep == "." or "#"
    870     //
    871     // ALMA & EVLA case (MS via ASDM) before3.1
    872     //
    873     // obsMode1=CALIBRATE_*
    874     //    obsMode2=ON_SOURCE: PONCAL
    875     //    obsMode2=OFF_SOURCE: POFFCAL
    876     // obsMode1=OBSERVE_TARGET
    877     //    obsMode2=ON_SOURCE: PON
    878     //    obsMode2=OFF_SOURCE: POFF
    879     string substr[2] ;
    880     int numSubstr = split( obsMode, substr, 2, "," ) ;
    881     //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
    882     //for ( int i = 0 ; i < numSubstr ; i++ )
    883     //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
    884     String obsType( substr[0] ) ;
    885     //os_ << "obsType = " << obsType << LogIO::POST ;
    886     Int epos = obsType.find_first_of( sep ) ;
    887     Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
    888     String obsMode1 = obsType.substr( 0, epos ) ;
    889     String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
    890     //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
    891     //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
    892     if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
    893       if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
    894       if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
    895     }
    896     else if ( obsMode1 == "OBSERVE_TARGET" ) {
    897       if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
    898       if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
    899     }
    900   }
    901   else if ( sep == "_" ) {
    902     // sep == "_"
    903     //
    904     // ALMA & EVLA case (MS via ASDM) after 3.2
    905     //
    906     // obsMode1=CALIBRATE_*
    907     //    obsMode2=ON_SOURCE: PONCAL
    908     //    obsMode2=OFF_SOURCE: POFFCAL
    909     // obsMode1=OBSERVE_TARGET
    910     //    obsMode2=ON_SOURCE: PON
    911     //    obsMode2=OFF_SOURCE: POFF
    912     string substr[2] ;
    913     int numSubstr = split( obsMode, substr, 2, "," ) ;
    914     //os_ << "numSubstr = " << numSubstr << LogIO::POST ;
    915     //for ( int i = 0 ; i < numSubstr ; i++ )
    916     //os_ << "substr[" << i << "] = " << substr[i] << LogIO::POST ;
    917     String obsType( substr[0] ) ;
    918     //os_ << "obsType = " << obsType << LogIO::POST ;
    919     string substr2[4] ;
    920     int numSubstr2 = split( obsType, substr2, 4, sep ) ;
    921     //Int epos = obsType.find_first_of( sep ) ;
    922     //Int nextpos = obsType.find_first_of( sep, epos+1 ) ;
    923     //String obsMode1 = obsType.substr( 0, epos ) ;
    924     //String obsMode2 = obsType.substr( epos+1, nextpos-epos-1 ) ;
    925     String obsMode1( substr2[0] ) ;
    926     String obsMode2( substr2[2] ) ;
    927     //os_ << "obsMode1 = " << obsMode1 << LogIO::POST ;
    928     //os_ << "obsMode2 = " << obsMode2 << LogIO::POST ;
    929     if ( obsMode1.find( "CALIBRATE" ) == 0 ) {
    930       if ( obsMode2 == "ON" ) srcType = SrcType::PONCAL ;
    931       if ( obsMode2 == "OFF" ) srcType = SrcType::POFFCAL ;
    932     }
    933     else if ( obsMode1 == "OBSERVE" ) {
    934       if ( obsMode2 == "ON" ) srcType = SrcType::PSON ;
    935       if ( obsMode2 == "OFF" ) srcType = SrcType::PSOFF ;
    936     }
    937   }
    938   else {
    939     if ( sig ) srcType = SrcType::SIG ;
    940     if ( ref ) srcType = SrcType::REF ;
    941   }
    942    
    943   //os_ << "srcType = " << srcType << LogIO::POST ;
    944   //double endSec = mathutil::gettimeofday_sec() ;
    945   //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    946   return srcType ;
    947 }
    948 
    949 //Vector<uInt> MSFiller::getPolNo( Int corrType )
    950 Block<uInt> MSFiller::getPolNo( Int corrType )
    951 {
    952   //double startSec = mathutil::gettimeofday_sec() ;
    953   //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
    954   Block<uInt> polno( 1 ) ;
    955 
    956   if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {
    957     polno = 0 ;
    958   }
    959   else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {
    960     polno = 1 ;
    961   }
    962   else if ( corrType == Stokes::U ) {
    963     polno = 2 ;
    964   }
    965   else if ( corrType == Stokes::V ) {
    966     polno = 3 ;
    967   }
    968   else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {
    969     polno.resize( 2 ) ;
    970     polno[0] = 2 ;
    971     polno[1] = 3 ;
    972   }
    973   else if ( corrType == Stokes::Plinear ) {
    974     polno[0] = 1 ;
    975   }
    976   else if ( corrType == Stokes::Pangle ) {
    977     polno[0] = 2 ;
    978   }
    979   else {
    980     polno = 99 ;
    981   }
    982   //os_ << "polno = " << polno << LogIO::POST ;
    983   //double endSec = mathutil::gettimeofday_sec() ;
    984   //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    985  
    986   return polno ;
    987 }
    988 
    989 String MSFiller::getPolType( Int corrType )
    990 {
    991   //double startSec = mathutil::gettimeofday_sec() ;
    992   //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
    993   String poltype = "" ;
    994 
    995   if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )
    996     poltype = "stokes" ;
    997   else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )
    998     poltype = "linear" ;
    999   else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )
    1000     poltype = "circular" ;
    1001   else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )
    1002     poltype = "linpol" ;
    1003 
    1004   //double endSec = mathutil::gettimeofday_sec() ;
    1005   //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1006   return poltype ;
    10071801}
    10081802
     
    11641958}
    11651959
    1166 void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr )
     1960void MSFiller::fillTcal()
    11671961{
    11681962  //double startSec = mathutil::gettimeofday_sec() ;
     
    11891983  }
    11901984
    1191   Table sctab = mstable_.sysCal() ;
     1985  Table &sctab = mstable_.sysCal() ;
    11921986  if ( sctab.nrow() == 0 ) {
    11931987    os_ << "No SYSCAL rows" << LogIO::POST ;
    11941988    return ;
    11951989  }
    1196   Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
    1197   if ( sctabsel.nrow() == 0 ) {
     1990  ROScalarColumn<Int> antCol( sctab, "ANTENNA_ID" ) ;
     1991  Vector<Int> ant = antCol.getColumn() ;
     1992  if ( allNE( ant, antenna_ ) ) {
    11981993    os_ << "No SYSCAL rows" << LogIO::POST ;
    11991994    return ;
    12001995  }
    1201   ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, colTcal_ ) ;
    1202   // return if any rows without Tcal value exists
     1996  ROTableColumn tcalCol( sctab, colTcal_ ) ;
    12031997  Bool notDefined = False ;
    1204   for ( uInt irow = 0 ; irow < sctabsel.nrow() ; irow++ ) {
    1205     if ( !tmpTcalCol->isDefined( irow ) ) {
     1998  for ( uInt irow = 0 ; irow < sctab.nrow() ; irow++ ) {
     1999    if ( ant[irow] == antenna_ && !tcalCol.isDefined( irow ) ) {
    12062000      notDefined = True ;
    12072001      break ;
     
    12102004  if ( notDefined ) {
    12112005    os_ << "No TCAL value" << LogIO::POST ;
    1212     delete tmpTcalCol ;
    12132006    table_->tcal().table().addRow(1,True) ;
    12142007    Vector<Float> defaultTcal( 1, 1.0 ) ;
     
    12172010    return ;
    12182011  }   
    1219   uInt npol = tmpTcalCol->shape( 0 )(0) ;
    1220   delete tmpTcalCol ;
    1221   //os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
     2012 
     2013  static const char *cols[] = {
     2014    "ANTENNA_ID", "FEED_ID", "SPECTRAL_WINDOW_ID", "TIME",
     2015    NULL
     2016  };
     2017  static const TypeManagerImpl<Int> tmInt;
     2018  static const TypeManagerImpl<Double> tmDouble;
     2019  static const TypeManager *const tms[] = {
     2020    &tmInt, &tmInt, &tmInt, &tmDouble, NULL
     2021  };
    12222022  Table tab = table_->tcal().table() ;
    1223   ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
    1224   uInt oldnr = 0 ;
    1225   uInt newnr = 0 ;
    1226   TableRow row( tab ) ;
    1227   TableRecord &trec = row.record() ;
    1228   RecordFieldPtr<uInt> idRF( trec, "ID" ) ;
    1229   RecordFieldPtr<String> timeRF( trec, "TIME" ) ;
    1230   RecordFieldPtr< Array<Float> > tcalRF( trec, "TCAL" ) ;
    1231   TableIterator iter0( sctabsel, "FEED_ID" ) ;
    1232   while( !iter0.pastEnd() ) {
    1233     Table t0 = iter0.table() ;
    1234     Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ;
    1235     TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
    1236     while( !iter1.pastEnd() ) {
    1237       Table t1 = iter1.table() ;
    1238       Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ;
    1239       tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ;
    1240       ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ;
    1241       Vector<uInt> idminmax( 2, oldnr ) ;
    1242       for ( uInt irow = 0 ; irow < t1.nrow() ; irow++ ) {
    1243         String sTime = MVTime( scTimeCol(irow) ).string( MVTime::YMD ) ;
    1244         *timeRF = sTime ;
    1245         uInt idx = oldnr ;
    1246         Matrix<Float> subtcal = (*tmpTcalCol)( irow ) ;
    1247         for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
    1248           *idRF = idx++ ;
    1249           //*tcalRF = subtcal.row( ipol ) ;
    1250           tcalRF.define( subtcal.row( ipol ) ) ;
    1251 
    1252           // commit row
    1253           tab.addRow() ;
    1254           row.put( tab.nrow()-1 ) ;
    1255 
    1256           newnr++ ;
    1257         }
    1258         idminmax[0] = oldnr ;
    1259         idminmax[1] = newnr - 1 ;
    1260         oldnr = newnr ;
    1261 
    1262         String key = keyTcal( feedId, spwId, sTime ) ;
    1263         tcalrec_.define( key, idminmax ) ;
    1264       }
    1265       delete tmpTcalCol ;
    1266       iter1++ ;
    1267     }
    1268     iter0++ ;
    1269   }
     2023  TcalVisitor visitor( sctab, tab, tcalrec_, antenna_ ) ;
     2024  visitor.setTcalColumn( colTcal_ ) ;
     2025 
     2026  traverseTable(sctab, cols, tms, &visitor);
    12702027
    12712028  //tcalrec_.print( std::cout ) ;
    12722029  //double endSec = mathutil::gettimeofday_sec() ;
    12732030  //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1274 }
    1275 
    1276 uInt MSFiller::getWeatherId( uInt idx, Double wtime )
    1277 {
    1278   //double startSec = mathutil::gettimeofday_sec() ;
    1279   //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
    1280   uInt nrow = mwTime_.size() ;
    1281   if ( nrow < 2 )
    1282     return 0 ;
    1283   uInt wid = nrow ;
    1284   if ( idx == 0 ) {
    1285     wid = 0 ;
    1286     Double tStart = mwTime_[wid]-0.5*mwInterval_[wid] ;
    1287     if ( wtime < tStart )
    1288       return wid ;
    1289   }
    1290   for ( uInt i = idx ; i < nrow-1 ; i++ ) {
    1291     Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
    1292     // use of INTERVAL column is problematic
    1293     // since there are "blank" time of weather monitoring
    1294     //Double tEnd = tStart + mwInterval_[i] ;
    1295     Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
    1296     //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
    1297     if ( wtime >= tStart && wtime <= tEnd ) {
    1298       wid = i ;
    1299       break ;
    1300     }
    1301   }
    1302   if ( wid == nrow ) {
    1303     uInt i = nrow - 1 ;
    1304     Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
    1305     Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
    1306     //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
    1307     if ( wtime >= tStart && wtime <= tEnd )
    1308       wid = i-1 ;
    1309     else
    1310       wid = i ;
    1311   }
    1312 
    1313   //if ( wid == nrow )
    1314   //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
    1315 
    1316   //double endSec = mathutil::gettimeofday_sec() ;
    1317   //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1318   return wid ;
    1319 }
    1320 
    1321 void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx )
    1322 {
    1323   //double startSec = mathutil::gettimeofday_sec() ;
    1324   //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
    1325 
    1326   if ( !isSysCal_ )
    1327     return ;
    1328 
    1329   uInt nrow = tidx.nelements() ;
    1330   if ( scTime.nelements() == 0 )
    1331     return ;
    1332   else if ( scTime.nelements() == 1 ) {
    1333     tidx[0] = 0 ;
    1334     return ;
    1335   }
    1336   uInt scnrow = scTime.nelements() ;
    1337   uInt idx = 0 ;
    1338   const Double half = 0.5e0 ;
    1339   // execute  binary search
    1340   idx = binarySearch( scTime, tcol[0].get( "s" ).getValue() ) ;
    1341   if ( idx != 0 )
    1342     idx -= 1 ;
    1343   for ( uInt i = 0 ; i < nrow ; i++ ) {
    1344     Double t = tcol[i].get( "s" ).getValue() ;
    1345     Double tsc = scTime[0].get( "s" ).getValue() ;
    1346     if ( t < tsc ) {
    1347       tidx[i] = 0 ;
    1348       continue ;
    1349     }
    1350     for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
    1351       Double tsc1 = scTime[j].get( "s" ).getValue() ;
    1352       Double dt1 = scInterval[j] ;
    1353       Double tsc2 = scTime[j+1].get( "s" ).getValue() ;
    1354       Double dt2 = scInterval[j+1] ;
    1355       if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {
    1356         tidx[i] = j ;
    1357         idx = j ;
    1358         break ;
    1359       }
    1360     }
    1361     if ( tidx[i] == -1 ) {
    1362 //       Double tsc = scTime[scnrow-1].get( "s" ).getValue() ;
    1363 //       Double dt = scInterval[scnrow-1] ;
    1364 //       if ( t <= tsc+0.5*dt ) {
    1365 //         tidx[i] = scnrow-1 ;
    1366 //       }
    1367       tidx[i] = scnrow-1 ;
    1368     }
    1369   }
    1370   //double endSec = mathutil::gettimeofday_sec() ;
    1371   //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;
    1372   return ;
    1373 }
    1374 
    1375 Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t )
    1376 {
    1377   //double startSec = mathutil::gettimeofday_sec() ;
    1378   //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
    1379   //if ( table_->tcal().table().nrow() == 0 ) {
    1380   if ( !isSysCal_ ) {
    1381     os_ << "No TCAL rows" << LogIO::POST ;
    1382     Block<uInt> tcalids( 4, 0 ) ;
    1383     return  tcalids ;
    1384   }   
    1385   //String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
    1386   String sctime = MVTime( t.getValue() ).string(MVTime::YMD) ;
    1387   String key = keyTcal( fid, spwid, sctime ) ;
    1388   if ( !tcalrec_.isDefined( key ) ) {
    1389     os_ << "No TCAL rows" << LogIO::POST ;
    1390     Block<uInt> tcalids( 4, 0 ) ;
    1391     return tcalids ;
    1392   }
    1393   Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;
    1394   uInt npol = ids[1] - ids[0] + 1 ;
    1395   Block<uInt> tcalids( npol ) ;
    1396   tcalids[0] = ids[0] ;
    1397   tcalids[1] = ids[1] ;
    1398   for ( uInt ipol = 2 ; ipol < npol ; ipol++ )
    1399     tcalids[ipol] = ids[0] + ipol - 1 ;
    1400 
    1401   //double endSec = mathutil::gettimeofday_sec() ;
    1402   //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1403   return tcalids ;
    1404 }
    1405 
    1406 uInt MSFiller::getDirection( uInt idx,
    1407                              Vector<Double> &dir,
    1408                              Vector<Double> &srate,
    1409                              String &ref,
    1410                              Vector<Double> &tcol,
    1411                              ROArrayColumn<Double> &dcol,
    1412                              Double t )
    1413 {
    1414   //double startSec = mathutil::gettimeofday_sec() ;
    1415   //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ;
    1416   //double time0 = mathutil::gettimeofday_sec() ;
    1417   //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ;
    1418   // assume that cols is sorted by TIME
    1419   Bool doInterp = False ;
    1420   uInt nrow = dcol.nrow() ;
    1421   if ( nrow == 0 )
    1422     return 0 ;
    1423   if ( idx == 0 ) {
    1424     uInt nrowb = 1000 ;
    1425     if ( nrow > nrowb ) {
    1426       uInt nblock = nrow / nrowb + 1 ;
    1427       for ( uInt iblock = 0 ; iblock < nblock ; iblock++ ) {
    1428         uInt high = min( nrowb, nrow-iblock*nrowb ) ;
    1429 
    1430         if ( tcol( high-1 ) < t ) {
    1431           idx = iblock * nrowb ;
    1432           continue ;
    1433         }
    1434 
    1435         Slice slice( iblock*nrowb, nrowb ) ;
    1436         Vector<Double> tarr = tcol( slice ) ;
    1437 
    1438         uInt bidx = binarySearch( tarr, t ) ;
    1439 
    1440         idx = iblock * nrowb + bidx ;
    1441         break ;
    1442       }
    1443     }
    1444     else {
    1445       idx = binarySearch( tcol, t ) ;
    1446     }
    1447   }
    1448   //double time1 = mathutil::gettimeofday_sec() ;
    1449   //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    1450   // ensure that tcol(idx) < t
    1451   //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
    1452   //time0 = mathutil::gettimeofday_sec() ;
    1453   //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ;
    1454   //while( tcol( idx ) * factor > t && idx > 0 )
    1455   while( tcol[idx] > t && idx > 0 )
    1456     idx-- ;
    1457   //os_ << "idx = " << idx << LogIO::POST ;
    1458 
    1459   // index search
    1460   for ( uInt i = idx ; i < nrow ; i++ ) {
    1461     Double tref = tcol[i] ;
    1462     if ( tref == t ) {
    1463       idx = i ;
    1464       break ;
    1465     }
    1466     else if ( tref > t ) {
    1467       if ( i == 0 ) {
    1468         idx = i ;
    1469       }
    1470       else {
    1471         idx = i-1 ;
    1472         doInterp = True ;
    1473       }
    1474       break ;
    1475     }
    1476     else {
    1477       idx = nrow - 1 ;
    1478     }
    1479   }
    1480   //time1 = mathutil::gettimeofday_sec() ;
    1481   //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    1482   //os_ << "searched idx = " << idx << LogIO::POST ;
    1483 
    1484   //time0 = mathutil::gettimeofday_sec() ;
    1485   //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ;
    1486   //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
    1487   //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
    1488   IPosition ip( dcol(idx).shape().nelements(), 0 ) ;
    1489   //os_ << "ip = " << ip << LogIO::POST ;
    1490   //ref = dmcol(idx)(ip).getRefString() ;
    1491   TableRecord trec = dcol.keywordSet() ;
    1492   Record rec = trec.asRecord( "MEASINFO" ) ;
    1493   ref = rec.asString( "Ref" ) ;
    1494   //os_ << "ref = " << ref << LogIO::POST ;
    1495   if ( doInterp ) {
    1496     //os_ << "do interpolation" << LogIO::POST ;
    1497     //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
    1498     Double tref0 = tcol[idx] ;
    1499     Double tref1 = tcol[idx+1] ;
    1500     Matrix<Double> mdir0 = dcol( idx ) ;
    1501     Matrix<Double> mdir1 = dcol( idx+1 ) ;
    1502     Vector<Double> dir0 = mdir0.column( 0 ) ;
    1503     //os_ << "dir0 = " << dir0 << LogIO::POST ;
    1504     Vector<Double> dir1 = mdir1.column( 0 ) ;
    1505     //os_ << "dir1 = " << dir1 << LogIO::POST ;
    1506     Double dt0 = t - tref0 ;
    1507     Double dt1 = tref1 - t ;
    1508     dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;
    1509     if ( mdir0.ncolumn() > 1 ) {
    1510       if ( dt0 >= dt1 )
    1511         srate.reference( mdir0.column( 1 ) ) ;
    1512       else
    1513         srate.reference( mdir1.column( 1 ) ) ;
    1514     }
    1515     //os_ << "dir = " << dir << LogIO::POST ;
    1516   }
    1517   else {
    1518     //os_ << "no interpolation" << LogIO::POST ;
    1519     Matrix<Double> mdir0 = dcol( idx ) ;
    1520     dir.reference( mdir0.column( 0 ) ) ;
    1521     if ( mdir0.ncolumn() > 1 )
    1522       srate.reference( mdir0.column( 1 ) ) ;
    1523   }
    1524 
    1525   //time1 = mathutil::gettimeofday_sec() ;
    1526   //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;
    1527   //double endSec = mathutil::gettimeofday_sec() ;
    1528   //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1529   return idx ;
    1530 }
    1531 
    1532 String MSFiller::keyTcal( Int feedid, Int spwid, String stime )
    1533 {
    1534   String sfeed = "FEED" + String::toString( feedid ) ;
    1535   String sspw = "SPW" + String::toString( spwid ) ;
    1536   return sfeed+":"+sspw+":"+stime ;
    1537 }
    1538 
    1539 uInt MSFiller::binarySearch( Vector<MEpoch> &timeList, Double target )
    1540 {
    1541   Int low = 0 ;
    1542   Int high = timeList.nelements() ;
    1543   uInt idx = 0 ;
    1544 
    1545   while ( low <= high ) {
    1546     idx = (Int)( 0.5 * ( low + high ) ) ;
    1547     Double t = timeList[idx].get( "s" ).getValue() ;
    1548     if ( t < target )
    1549       low = idx + 1 ;
    1550     else if ( t > target )
    1551       high = idx - 1 ;
    1552     else {
    1553       return idx ;
    1554     }
    1555   }
    1556 
    1557   idx = max( 0, min( low, high ) ) ;
    1558 
    1559   return idx ;
    1560 }
    1561 
    1562 uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target )
    1563 {
    1564   Int low = 0 ;
    1565   Int high = timeList.nelements() ;
    1566   uInt idx = 0 ;
    1567 
    1568   while ( low <= high ) {
    1569     idx = (Int)( 0.5 * ( low + high ) ) ;
    1570     Double t = timeList[idx] ;
    1571     if ( t < target )
    1572       low = idx + 1 ;
    1573     else if ( t > target )
    1574       high = idx - 1 ;
    1575     else {
    1576       return idx ;
    1577     }
    1578   }
    1579 
    1580   idx = max( 0, min( low, high ) ) ;
    1581 
    1582   return idx ;
    15832031}
    15842032
     
    16032051}
    16042052
    1605 void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp,
    1606                                          Cube<Bool> &fl,
    1607                                          Table &tab,
    1608                                          Int &npol,
    1609                                          Int &nchan,
    1610                                          Int &nrow,
    1611                                          Vector<Int> &corrtype )
    1612 {
    1613   //double startSec = mathutil::gettimeofday_sec() ;
    1614   //os_ << "start MSFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ; 
    1615   if ( isFloatData_ ) {
    1616     ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
    1617     ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ;
    1618     mFloatDataCol.getColumn( sp ) ;
    1619     mFlagCol.getColumn( fl ) ;
    1620   }
    1621   else if ( isData_ ) {
    1622     sp.resize( npol, nchan, nrow ) ;
    1623     fl.resize( npol, nchan, nrow ) ;
    1624     ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ;
    1625     ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ;
    1626     if ( npol < 3 ) {
    1627       Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ;
    1628       IPosition start( 3, 0, 0, 0 ) ;
    1629       IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ;
    1630       IPosition inc( 3, 2, 1, 1 ) ;
    1631       sp = tmp( start, end, inc ) ;
    1632       fl = mFlagCol.getColumn() ;
    1633     }
    1634     else {
    1635       for ( Int irow = 0 ; irow < nrow ; irow++ ) {
    1636         Bool crossOK = False ;
    1637         Matrix<Complex> mSp = mDataCol( irow ) ;
    1638         Matrix<Bool> mFl = mFlagCol( irow ) ;
    1639         Matrix<Float> spxy = sp.xyPlane( irow ) ;
    1640         Matrix<Bool> flxy = fl.xyPlane( irow ) ;
    1641         for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
    1642           if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
    1643                || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
    1644             if ( !crossOK ) {
    1645               Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
    1646               IPosition start( 1, 0 ) ;
    1647               IPosition end( 1, 2*nchan-1 ) ;
    1648               IPosition inc( 1, 2 ) ;
    1649               spxy.row( ipol ) = tmp( start, end, inc ) ;
    1650               flxy.row( ipol ) = mFl.row( ipol ) ;
    1651               start = IPosition( 1, 1 ) ;
    1652               spxy.row( ipol+1 ) = tmp( start, end, inc ) ;
    1653               flxy.row( ipol+1 ) = mFl.row( ipol ) ;
    1654               if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) {
    1655                 spxy.row( ipol+1 ) = spxy.row( ipol+1 ) * (Float)-1.0 ;
    1656               }
    1657               crossOK = True ;
    1658             }
    1659           }
    1660           else {
    1661             Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ;
    1662             IPosition start( 1, 0 ) ;
    1663             IPosition end( 1, 2*nchan-1 ) ;
    1664             IPosition inc( 1, 2 ) ;
    1665             spxy.row( ipol ) = tmp( start, end, inc ) ;
    1666             flxy.row( ipol ) = mFl.row( ipol ) ;
    1667           }
    1668         }
    1669       }
    1670     }
    1671   }
    1672   //double endSec = mathutil::gettimeofday_sec() ;
    1673   //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1674 }
    1675 
    1676 uInt MSFiller::getDirection( uInt idx,
    1677                              Vector<Double> &dir,
    1678                              Vector<Double> &azel,
    1679                              Vector<Double> &srate,
    1680                              Vector<Double> &ptcol,
    1681                              ROArrayColumn<Double> &pdcol,
    1682                              MEpoch &t,
    1683                              MPosition &antpos )
    1684 {
    1685   //double startSec = mathutil::gettimeofday_sec() ;
    1686   //os_ << "start MSFiller::getDirection2() startSec=" << startSec << LogIO::POST ; 
    1687   String refString ;
    1688   MDirection::Types dirType ;
    1689   uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ;
    1690   MDirection::getType( dirType, refString ) ;
    1691   MeasFrame mf( t, antpos ) ;
    1692   if ( refString == "J2000" ) {
    1693     MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
    1694     azel = toazel( dir ).getAngle("rad").getValue() ;
    1695   }
    1696   else if ( refString(0,4) == "AZEL" ) {
    1697     azel = dir.copy() ;
    1698     MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
    1699     dir = toj2000( dir ).getAngle("rad").getValue() ;
    1700   }
    1701   else {
    1702     MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
    1703     azel = toazel( dir ).getAngle("rad").getValue() ;
    1704     MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
    1705     dir = toj2000( dir ).getAngle("rad").getValue() ;
    1706   }
    1707   if ( srate.size() == 0 ) {
    1708     srate.resize( 2 ) ;
    1709     srate = 0.0 ;
    1710   }
    1711   //double endSec = mathutil::gettimeofday_sec() ;
    1712   //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1713   return diridx ;
    1714 }
    1715 
    1716 void MSFiller::getSourceDirection( Vector<Double> &dir,
    1717                                    Vector<Double> &azel,
    1718                                    Vector<Double> &srate,
    1719                                    MEpoch &t,
    1720                                    MPosition &antpos,
    1721                                    Vector<MDirection> &srcdir )
    1722 {
    1723   //double startSec = mathutil::gettimeofday_sec() ;
    1724   //os_ << "start MSFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ;
    1725   Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ;
    1726   if ( srcdir.nelements() > 1 )
    1727     srate = srcdir[1].getAngle( "rad" ).getValue() ;
    1728   String ref = srcdir[0].getRefString() ;
    1729   MDirection::Types dirType ;
    1730   MDirection::getType( dirType, ref ) ;
    1731   MeasFrame mf( t, antpos ) ;
    1732   if ( ref != "J2000" ) {
    1733     MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
    1734     dir = toj2000( defaultDir ).getAngle("rad").getValue() ;
    1735   }
    1736   else
    1737     dir = defaultDir ;
    1738   MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
    1739   azel = toazel( defaultDir ).getAngle("rad").getValue() ;
    1740   //double endSec = mathutil::gettimeofday_sec() ;
    1741   //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1742 }
    1743 
    17442053void MSFiller::initHeader( STHeader &header )
    17452054{
     
    17522061  header.obstype = "" ;
    17532062  header.antennaname = "" ;
    1754   header.antennaposition.resize( 0 ) ;
     2063  header.antennaposition.resize( 3 ) ;
    17552064  header.equinox = 0.0 ;
    17562065  header.freqref = "" ;
     
    17632072}
    17642073
    1765 String MSFiller::asString( String name,
    1766                            uInt idx,
    1767                            Table tab,
    1768                            boost::object_pool<ROTableColumn> *pool )
    1769 {
    1770   ROTableColumn *col = pool->construct( tab, name ) ;
    1771   String v = col->asString( idx ) ;
    1772   pool->destroy( col ) ;
    1773   return v ;
    1774 }
    1775 
    1776 Bool MSFiller::asBool( String name,
    1777                          uInt idx,
    1778                          Table &tab,
    1779                          boost::object_pool<ROTableColumn> *pool )
    1780 {
    1781   ROTableColumn *col = pool->construct( tab, name ) ;
    1782   Bool v = col->asBool( idx ) ;
    1783   pool->destroy( col ) ;
    1784   return v ;
    1785 }
    1786 
    1787 uInt MSFiller::asuInt( String name,
    1788                          uInt idx,
    1789                          Table &tab,
    1790                          boost::object_pool<ROTableColumn> *pool )
    1791 {
    1792   ROTableColumn *col = pool->construct( tab, name ) ;
    1793   uInt v = col->asuInt( idx ) ;
    1794   pool->destroy( col ) ;
    1795   return v ;
    1796 }
    1797 
    1798 Int MSFiller::asInt( String name,
    1799                         uInt idx,
    1800                         Table &tab,
    1801                         boost::object_pool<ROTableColumn> *pool )
    1802 {
    1803   ROTableColumn *col = pool->construct( tab, name ) ;
    1804   Int v = col->asInt( idx ) ;
    1805   pool->destroy( col ) ;
    1806   return v ;
    1807 }
    1808 
    1809 Float MSFiller::asFloat( String name,
    1810                           uInt idx,
    1811                           Table &tab,
    1812                           boost::object_pool<ROTableColumn> *pool )
    1813 {
    1814   ROTableColumn *col = pool->construct( tab, name ) ;
    1815   Float v = col->asfloat( idx ) ;
    1816   pool->destroy( col ) ;
    1817   return v ;
    1818 }
    1819 
    1820 Double MSFiller::asDouble( String name,
    1821                            uInt idx,
    1822                            Table &tab,
    1823                            boost::object_pool<ROTableColumn> *pool )
    1824 {
    1825   ROTableColumn *col = pool->construct( tab, name ) ;
    1826   Double v = col->asdouble( idx ) ;
    1827   pool->destroy( col ) ;
    1828   return v ;
    1829 }
    1830 
    1831 void MSFiller::sourceInfo( Int sourceId,
    1832                            Int spwId,
    1833                            String &name,
    1834                            MDirection &direction,
    1835                            Vector<casa::Double> &properMotion,
    1836                            Vector<casa::Double> &restFreqs,
    1837                            Vector<casa::String> &transitions,
    1838                            Vector<casa::Double> &sysVels,
    1839                            boost::object_pool<ROTableColumn> *tpoolr )
    1840 {
    1841   //double startSec = mathutil::gettimeofday_sec() ;
    1842   //os_ << "start MSFiller::sourceInfo() startSec=" << startSec << LogIO::POST ;
    1843 
    1844   MSSource srctab = mstable_.source() ;
    1845   MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;
    1846   if ( srctabSel.nrow() == 0 ) {
    1847     srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;
    1848   }
    1849   Int numLines = 0 ;
    1850   if ( srctabSel.nrow() > 0 ) {
    1851     // source name
    1852     name = asString( "NAME", 0, srctabSel, tpoolr ) ;
    1853    
    1854     // source proper motion
    1855     ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ;
    1856     properMotion = roArrDCol( 0 ) ;
    1857    
    1858     // source direction as MDirection object
    1859     ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ;
    1860     direction = tmpMeasCol( 0 ) ;
    1861    
    1862     // number of lines
    1863     numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ;
    1864   }
    1865   else {
    1866     name = "" ;
    1867     properMotion = Vector<Double>( 2, 0.0 ) ;
    1868     direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ;
    1869   }
    1870 
    1871   restFreqs.resize( numLines ) ;
    1872   transitions.resize( numLines ) ;
    1873   sysVels.resize( numLines ) ;
    1874   if ( numLines > 0 ) {
    1875     if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) {
    1876       ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ;
    1877       Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ;
    1878       for ( int i = 0 ; i < numLines ; i++ ) {
    1879         restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
    1880       }
    1881     }
    1882     //os_ << "restFreqs = " << restFreqs << LogIO::POST ;
    1883     if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
    1884       ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ;
    1885       if ( transitionCol.isDefined( 0 ) )
    1886         transitions = transitionCol( 0 ) ;
    1887       //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ;
    1888     }
    1889     if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) {
    1890       ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ;
    1891       sysVels = roArrDCol( 0 ) ;
    1892     }
    1893   }
    1894  
    1895   //double endSec = mathutil::gettimeofday_sec() ;
    1896   //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1897 }
    1898 
    1899 void MSFiller::spectralSetup( Int spwId,
    1900                               MEpoch &me,
    1901                               MPosition &mp,
    1902                               MDirection &md,
    1903                               Double &refpix,
    1904                               Double &refval,
    1905                               Double &increment,
    1906                               Int &nchan,
    1907                               String &freqref,
    1908                               Double &reffreq,
    1909                               Double &bandwidth,
    1910                               boost::object_pool<ROTableColumn> *tpoolr )
    1911 {
    1912   //double startSec = mathutil::gettimeofday_sec() ;
    1913   //os_ << "start MSFiller::spectralSetup() startSec=" << startSec << LogIO::POST ;
    1914 
    1915   MSSpectralWindow spwtab = mstable_.spectralWindow() ;
    1916   MeasFrame mf( me, mp, md ) ;
    1917   MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ;
    1918   Bool even = False ;
    1919   if ( (nchan/2)*2 == nchan ) even = True ;
    1920   ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ;
    1921   Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ;
    1922   if ( nchan != 4 )
    1923     bandwidth = max( bandwidth, totbw ) ;
    1924   if ( freqref == "" && nchan != 4)
    1925     //sdh.freqref = MFrequency::showType( freqRef ) ;
    1926     freqref = "LSRK" ;
    1927   if ( reffreq == -1.0 && nchan != 4 ) {
    1928     tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ;
    1929     Quantum<Double> qreffreq = tmpQuantCol( spwId ) ;
    1930     if ( freqRef == MFrequency::LSRK ) {
    1931       reffreq = qreffreq.getValue("Hz") ;
    1932     }
    1933     else {
    1934       MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
    1935       reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
    1936     }
    1937   }
    1938   Int refchan = nchan / 2 ;
    1939   IPosition refip( 1, refchan ) ;
    1940   refpix = 0.5*(nchan-1) ;
    1941   refval = 0.0 ;
    1942   ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ;
    1943   ROTableColumn netSidebandCol( spwtab, "NET_SIDEBAND" ) ;
    1944   Int netSideband = netSidebandCol.asInt( spwId ) ;
    1945   increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ;
    1946   //           os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
    1947   sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ;
    1948   Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ;
    1949   if ( ( nchan > 1 &&
    1950          chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz")  )
    1951        || ( nchan == 1 && netSideband == 1 ) ) 
    1952     increment *= -1.0 ;
    1953   if ( freqRef == MFrequency::LSRK ) {
    1954     if ( even ) {
    1955       IPosition refip0( 1, refchan-1 ) ;
    1956       Double refval0 = chanFreqs(refip0).getValue("Hz") ;
    1957       Double refval1 = chanFreqs(refip).getValue("Hz") ;
    1958       refval = 0.5 * ( refval0 + refval1 ) ;
    1959     }
    1960     else {
    1961       refval = chanFreqs(refip).getValue("Hz") ;
    1962     }
    1963   }
    1964   else {
    1965     MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
    1966     if ( even ) {
    1967       IPosition refip0( 1, refchan-1 ) ;
    1968       Double refval0 = chanFreqs(refip0).getValue("Hz") ;
    1969       Double refval1 = chanFreqs(refip).getValue("Hz") ;
    1970       refval = 0.5 * ( refval0 + refval1 ) ;
    1971       refval = tolsr( refval ).get("Hz").getValue() ;
    1972     }
    1973     else {
    1974       refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
    1975     }
    1976   }
    1977  
    1978   //double endSec = mathutil::gettimeofday_sec() ;
    1979   //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    1980 }
    1981 
    1982 } ;
    1983 
     2074};
  • trunk/src/MSFiller.h

    r2289 r2291  
    1616// STL
    1717#include <string>
    18 
    19 // Boost
    20 #include <boost/pool/object_pool.hpp>
    2118
    2219// AIPS++
     
    2724#include <casa/Arrays/Cube.h>
    2825#include <casa/Logging/LogIO.h>
    29 
     26#include <casa/Containers/RecordField.h>
    3027#include <casa/Containers/Record.h>
    3128#include <casa/Containers/Block.h>
     29#include <casa/Quanta/MVTime.h>
    3230
    3331#include <ms/MeasurementSets/MeasurementSet.h>
    3432#include <ms/MeasurementSets/MSPointing.h>
    3533
     34#include <tables/Tables/ScalarColumn.h>
     35#include <tables/Tables/ArrayColumn.h>
     36#include <tables/Tables/TableRow.h>
     37
     38#include <measures/TableMeasures/ScalarMeasColumn.h>
     39#include <measures/TableMeasures/ArrayMeasColumn.h>
     40#include <measures/TableMeasures/ScalarQuantColumn.h>
     41#include <measures/TableMeasures/ArrayQuantColumn.h>
     42
     43#include "TableTraverse.h"
    3644#include "Scantable.h"
     45#include "MathUtils.h"
     46
     47using namespace casa;
    3748
    3849namespace asap
    3950{
     51
     52class MSFillerUtils {
     53protected:
     54  template<class T> void getScalar( const String &name,
     55                                    const uInt &idx,
     56                                    const Table &tab,
     57                                    T &val )
     58  {
     59    ROScalarColumn<T> col( tab, name ) ;
     60    val = col( idx ) ;
     61  }
     62  template<class T> void getArray( const String &name,
     63                                   const uInt &idx,
     64                                   const Table &tab,
     65                                   Array<T> &val )
     66  {
     67    ROArrayColumn<T> col( tab, name ) ;
     68    val = col( idx ) ;
     69  }
     70  template<class T> void getScalarMeas( const String &name,
     71                                        const uInt &idx,
     72                                        const Table &tab,
     73                                        T &val )
     74  {
     75    ROScalarMeasColumn<T> measCol( tab, name ) ;
     76    val = measCol( idx ) ;
     77  }
     78  template<class T> void getArrayMeas( const String &name,
     79                                       const uInt &idx,
     80                                       const Table &tab,
     81                                       Array<T> &val )
     82  {
     83    ROArrayMeasColumn<T> measCol( tab, name ) ;
     84    val = measCol( idx ) ;
     85  }
     86  template<class T> void getScalarQuant( const String &name,
     87                                         const uInt &idx,
     88                                         const Table &tab,
     89                                         Quantum<T> &val )
     90  {
     91    ROScalarQuantColumn<T> quantCol( tab, name ) ;
     92    val = quantCol( idx ) ;
     93  }
     94  template<class T> void getArrayQuant( const String &name,
     95                                        const uInt &idx,
     96                                        const Table &tab,
     97                                        Array< Quantum<T> > &val )
     98  {
     99    ROArrayQuantColumn<T> quantCol( tab, name ) ;
     100    val = quantCol( idx ) ;
     101  }
     102//   template<class T> void putField( const String &name,
     103//                                    TableRecord &rec,
     104//                                    T &val )
     105//   {
     106//     RecordFieldPtr<T> rf( rec, name ) ;
     107//     *rf = val ;
     108//   }
     109//   template<class T> void defineField( const String &name,
     110//                                       TableRecord &rec,
     111//                                       T &val )
     112//   {
     113//     RecordFieldPtr<T> rf( rec, name ) ;
     114//     rf.define( val ) ;
     115//   }
     116  template<class T> T interp( Double x0, Double x1, Double x, T y0, T y1 )
     117  {
     118    Double dx0 = x - x0 ;
     119    Double dx1 = x1 - x ;
     120    return ( y0 * dx1 + y1 * dx0 ) / ( x1 - x0 ) ;
     121  }
     122  String keyTcal( const Int &feedid, const Int &spwid, const Double &time )
     123  {
     124    String stime = MVTime( Quantity(time,Unit("s")) ).string( MVTime::YMD ) ;
     125    String sfeed = "FEED" + String::toString( feedid ) ;
     126    String sspw = "SPW" + String::toString( spwid ) ;
     127    return sfeed+":"+sspw+":"+stime ;
     128  }
     129  String keyTcal( const Int &feedid, const Int &spwid, const String &stime )
     130  {
     131    String sfeed = "FEED" + String::toString( feedid ) ;
     132    String sspw = "SPW" + String::toString( spwid ) ;
     133    return sfeed+":"+sspw+":"+stime ;
     134  }
     135};
    40136
    41137class MSFiller
    42138{
    43139public:
    44   explicit MSFiller(casa::CountedPtr<Scantable> stable) ;
     140  explicit MSFiller(CountedPtr<Scantable> stable) ;
    45141  virtual ~MSFiller() ;
    46142 
    47   virtual bool open(const std::string& filename, const casa::Record& rec) ;
     143  virtual bool open(const std::string& filename, const Record& rec) ;
    48144  virtual void fill() ;
    49145  virtual void close() ;
     
    65161  //void fillHistory() ;
    66162  //void fillFit() ;
    67   void fillTcal( boost::object_pool<casa::ROTableColumn> *poolr ) ;
    68 
    69   // get SRCTYPE from STATE_ID
    70   casa::Int getSrcType( casa::Int stateId, boost::object_pool<casa::ROTableColumn> *pool ) ;
    71 
    72   // get POLNO from CORR_TYPE
    73   casa::Block<casa::uInt> getPolNo( casa::Int corrType ) ;
    74 
    75   // get poltype from CORR_TYPE
    76   casa::String getPolType( casa::Int corrType ) ;
    77 
    78   // get WEATHER_ID
    79   casa::uInt getWeatherId( casa::uInt idx, casa::Double wtime ) ;
    80 
    81   // get time stamp in SYSCAL table
    82   // assume that tab is selected by ANTENNA_ID, FEED_ID, SPECTRAL_WINDOW_ID
    83   // and sorted by TIME
    84   void getSysCalTime( casa::Vector<casa::MEpoch> &scTimeIn, casa::Vector<casa::Double> &scInterval, casa::Block<casa::MEpoch> &tcol, casa::Block<casa::Int> &tidx ) ;
    85 
    86   // get TCAL_ID
    87   casa::Block<casa::uInt> getTcalId( casa::Int feedId, casa::Int spwId, casa::MEpoch &t ) ;
    88 
    89   // get direction for DIRECTION, AZIMUTH, and ELEVATION columns
    90   casa::uInt getDirection( casa::uInt idx,
    91                            casa::Vector<casa::Double> &dir,
    92                            casa::Vector<casa::Double> &srate,
    93                            casa::String &ref,
    94                            casa::Vector<casa::Double> &ptcol,
    95                            casa::ROArrayColumn<casa::Double> &pdcol,
    96                            casa::Double t ) ;
    97   casa::uInt getDirection( casa::uInt idx,
    98                            casa::Vector<casa::Double> &dir,
    99                            casa::Vector<casa::Double> &azel,
    100                            casa::Vector<casa::Double> &srate,
    101                            casa::Vector<casa::Double> &ptcol,
    102                            casa::ROArrayColumn<casa::Double> &pdcol,
    103                            casa::MEpoch &t, casa::MPosition &antpos ) ;
    104   void getSourceDirection( casa::Vector<casa::Double> &dir,
    105                            casa::Vector<casa::Double> &azel,
    106                            casa::Vector<casa::Double> &srate,
    107                            casa::MEpoch &t,
    108                            casa::MPosition &antpos,
    109                            casa::Vector<casa::MDirection> &srcdir ) ;
     163  void fillTcal() ;
    110164
    111165  // create key for TCAL table
    112   casa::String keyTcal( casa::Int feedid, casa::Int spwid, casa::String stime ) ;
    113 
    114   // binary search
    115   casa::uInt binarySearch( casa::Vector<casa::MEpoch> &timeList, casa::Double target ) ;
    116   casa::uInt binarySearch( casa::Vector<casa::Double> &timeList, casa::Double target ) ;
    117 
    118   // tool for HPC
    119 //   double gettimeofday_sec() ;
     166  String keyTcal( Int feedid, Int spwid, String stime ) ;
    120167
    121168  // get frequency frame
    122169  std::string getFrame() ;
    123170
    124   // reshape SPECTRA and FLAGTRA
    125   void reshapeSpectraAndFlagtra( casa::Cube<casa::Float> &sp,
    126                                  casa::Cube<casa::Bool> &fl,
    127                                  casa::Table &tab,
    128                                  casa::Int &npol,
    129                                  casa::Int &nchan,
    130                                  casa::Int &nrow,
    131                                  casa::Vector<casa::Int> &corrtype ) ;
    132 
    133171  // initialize header
    134172  void initHeader( STHeader &header ) ;
    135173
    136   // get value from table column using object pool
    137   casa::String asString( casa::String name,
    138                          casa::uInt idx,
    139                          casa::Table tab,
    140                          boost::object_pool<casa::ROTableColumn> *pool ) ;
    141   casa::Bool asBool( casa::String name,
    142                      casa::uInt idx,
    143                      casa::Table &tab,
    144                      boost::object_pool<casa::ROTableColumn> *pool ) ;
    145   casa::Int asInt( casa::String name,
    146                    casa::uInt idx,
    147                    casa::Table &tab,
    148                    boost::object_pool<casa::ROTableColumn> *pool ) ;
    149   casa::uInt asuInt( casa::String name,
    150                      casa::uInt idx,
    151                      casa::Table &tab,
    152                      boost::object_pool<casa::ROTableColumn> *pool ) ;
    153   casa::Float asFloat( casa::String name,
    154                        casa::uInt idx,
    155                        casa::Table &tab,
    156                        boost::object_pool<casa::ROTableColumn> *pool ) ;
    157   casa::Double asDouble( casa::String name,
    158                          casa::uInt idx,
    159                          casa::Table &tab,
    160                          boost::object_pool<casa::ROTableColumn> *pool ) ;
    161 
    162   void sourceInfo( casa::Int sourceId,
    163                    casa::Int spwId,
    164                    casa::String &name,
    165                    casa::MDirection &direction,
    166                    casa::Vector<casa::Double> &properMotion,
    167                    casa::Vector<casa::Double> &restFreqs,
    168                    casa::Vector<casa::String> &transitions,
    169                    casa::Vector<casa::Double> &sysVels,
    170                    boost::object_pool<casa::ROTableColumn> *pool ) ;
    171 
    172   void spectralSetup( casa::Int spwId,
    173                       casa::MEpoch &me,
    174                       casa::MPosition &mp,
    175                       casa::MDirection &md,
    176                       casa::Double &refpix,
    177                       casa::Double &refval,
    178                       casa::Double &increment,
    179                       casa::Int &nchan,
    180                       casa::String &freqref,
    181                       casa::Double &reffreq,
    182                       casa::Double &bandwidth,
    183                       boost::object_pool<casa::ROTableColumn> *pool ) ;
    184 
    185   casa::CountedPtr<Scantable> table_ ;
    186   casa::MeasurementSet mstable_ ;
    187   casa::String tablename_ ;
    188   casa::Int antenna_ ;
    189   casa::String antennaStr_ ;
    190   casa::Bool getPt_ ;
    191 
    192   casa::Bool isFloatData_ ;
    193   casa::Bool isData_ ;
    194 
    195   casa::Bool isDoppler_ ;
    196   casa::Bool isFlagCmd_ ;
    197   casa::Bool isFreqOffset_ ;
    198   casa::Bool isHistory_ ;
    199   casa::Bool isProcessor_ ;
    200   casa::Bool isSysCal_ ;
    201   casa::Bool isWeather_ ;
    202 
    203   casa::String colTsys_ ;
    204   casa::String colTcal_ ;
    205 
    206   casa::LogIO os_ ;
    207  
    208   casa::Vector<casa::Double> mwTime_ ;
    209   casa::Vector<casa::Double> mwInterval_ ;
    210   casa::Vector<casa::uInt> mwIndex_ ;
     174  CountedPtr<Scantable> table_ ;
     175  MeasurementSet mstable_ ;
     176  String tablename_ ;
     177  Int antenna_ ;
     178  String antennaStr_ ;
     179  Bool getPt_ ;
     180
     181  Bool isFloatData_ ;
     182  Bool isData_ ;
     183
     184  Bool isDoppler_ ;
     185  Bool isFlagCmd_ ;
     186  Bool isFreqOffset_ ;
     187  Bool isHistory_ ;
     188  Bool isProcessor_ ;
     189  Bool isSysCal_ ;
     190  Bool isWeather_ ;
     191
     192  String colTsys_ ;
     193  String colTcal_ ;
     194
     195  LogIO os_ ;
     196 
     197  Vector<Double> mwTime_ ;
     198  Vector<Double> mwInterval_ ;
     199  Vector<uInt> mwIndex_ ;
    211200
    212201  // Record for TCAL_ID
     
    214203  //           "SPW1": Vector<uInt>
    215204  //  ...
    216   casa::Record tcalrec_ ;
    217 
    218   //casa::ROTableColumn *scCol_ ;
     205  Record tcalrec_ ;
     206  //map< String,Vector<uInt> > tcalrec_ ;
    219207};
    220208
  • trunk/src/MSWriter.cpp

    r2259 r2291  
    1111//
    1212//
     13#include <assert.h>
    1314
    1415#include <casa/OS/File.h>
     
    1718#include <casa/OS/SymLink.h>
    1819#include <casa/BasicSL/String.h>
    19 #include <casa/Containers/RecordField.h>
    2020#include <casa/Arrays/Cube.h>
    2121
     
    3939#include "STTcal.h"
    4040#include "MathUtils.h"
    41 
    42 // #include <ctime>
    43 // #include <sys/time.h>
    44 
     41#include "TableTraverse.h"
    4542
    4643using namespace casa ;
     
    4845
    4946namespace asap {
    50 // double MSWriter::gettimeofday_sec()
    51 // {
    52 //   struct timeval tv ;
    53 //   gettimeofday( &tv, NULL ) ;
    54 //   return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
    55 // }
     47
     48class PolarizedComponentHolder {
     49public:
     50  PolarizedComponentHolder()
     51    : nchan(0),
     52      maxnpol(4)
     53  {
     54    reset() ;
     55  }
     56  PolarizedComponentHolder( uInt n )
     57    : nchan(n),
     58      maxnpol(4)
     59  {
     60    reset() ;
     61  }
     62
     63  void reset()
     64  {
     65    npol = 0 ;
     66    data.clear() ;
     67    flag.clear() ;
     68    flagrow = False ;
     69    polnos.resize() ;
     70  }
     71
     72  void accumulate( uInt id, Vector<Float> sp, Vector<Bool> fl, Bool flr )
     73  {
     74    map< uInt,Vector<Float> >::iterator itr = data.find( id ) ;
     75    if ( id < maxnpol && itr == data.end() ) {
     76      addPol( id ) ;
     77      accumulateData( id, sp ) ;
     78      accumulateFlag( id, fl ) ;
     79      accumulateFlagRow( flr ) ;
     80      npol++ ;
     81    }
     82  }
     83
     84  void setNchan( uInt n ) { nchan = n ; }
     85  uInt nPol() { return npol ; }
     86  uInt nChan() { return nchan ; }
     87  Vector<uInt> polNos() { return polnos ; }
     88  Vector<Float> getWeight() { return Vector<Float>( npol, 1.0 ) ; }
     89  Vector<Float> getSigma() { return Vector<Float>( npol, 1.0 ) ; }
     90  Bool getFlagRow() { return flagrow ; }
     91  Cube<Bool> getFlagCategory() { return Cube<Bool>( npol, nchan, 1, False ) ; }
     92  Matrix<Float> getData()
     93  {
     94    Matrix<Float> v( npol, nchan ) ;
     95    for ( map< uInt,Vector<Float> >::iterator i = data.begin() ; i != data.end() ; i++ ) {
     96      v.row( i->first ) =  i->second ;
     97    }
     98    return v ;
     99  }
     100  Matrix<Bool> getFlag()
     101  {
     102    Matrix<Bool> v( npol, nchan ) ;
     103    for ( map< uInt,Vector<Bool> >::iterator i = flag.begin() ; i != flag.end() ; i++ ) {
     104      v.row( i->first ) = i->second ;
     105    }
     106    return v ;
     107  }
     108  Matrix<Complex> getComplexData()
     109  {
     110    Matrix<Complex> v( npol, nchan ) ;
     111    Matrix<Float> dummy( 2, nchan, 0.0 ) ;
     112    map< uInt,Vector<Float> >::iterator itr0 = data.find( 0 ) ;
     113    map< uInt,Vector<Float> >::iterator itr1 = data.find( 1 ) ;
     114    if ( itr0 != data.end() ) {
     115      dummy.row( 0 ) = itr0->second ;
     116      v.row( 0 ) = RealToComplex( dummy ) ;
     117    }
     118    if ( itr1 != data.end() ) {
     119      dummy.row( 0 ) = itr1->second ;
     120      v.row( npol-1 ) = RealToComplex( dummy ) ;
     121    }
     122    itr0 = data.find( 2 ) ;
     123    itr1 = data.find( 3 ) ;
     124    if ( itr0 != data.end() && itr1 != data.end() ) {
     125      dummy.row( 0 ) = itr0->second ;
     126      dummy.row( 1 ) = itr1->second ;
     127      v.row( 1 ) = RealToComplex( dummy ) ;
     128      v.row( 2 ) = conj( v.row( 1 ) ) ;
     129    }
     130    return v ;
     131  }
     132
     133  Matrix<Bool> getComplexFlag()
     134  {
     135    Matrix<Bool> tmp = getFlag() ;
     136    Matrix<Bool> v( npol, nchan ) ;
     137    v.row( 0 ) = tmp.row( 0 ) ;
     138    if ( npol == 2 ) {
     139      v.row( npol-1 ) = tmp.row( 1 ) ;
     140    }
     141    else if ( npol > 2 ) {
     142      v.row( npol-1 ) = tmp.row( 1 ) ;
     143      v.row( 1 ) = tmp.row( 2 ) || tmp.row( 3 ) ;
     144      v.row( 2 ) = v.row( 1 ) ;
     145    }
     146    return v ;
     147  }
     148 
     149private:
     150  void accumulateData( uInt &id, Vector<Float> &v )
     151  {
     152    data.insert( pair< uInt,Vector<Float> >( id, v ) ) ;
     153  }
     154    void accumulateFlag( uInt &id, Vector<Bool> &v )
     155  {
     156    flag.insert( pair< uInt,Vector<Bool> >( id, v ) ) ;
     157  }
     158  void accumulateFlagRow( Bool &v )
     159  {
     160    flagrow |= v ;
     161  }
     162  void addPol( uInt id )
     163  {
     164    uInt i = polnos.nelements() ;
     165    polnos.resize( i+1, True ) ;
     166    polnos[i] = id ;
     167  }
     168
     169  uInt nchan;
     170  const uInt maxnpol;
     171  uInt npol;
     172  Vector<uInt> polnos;
     173
     174  map< uInt,Vector<Float> > data;
     175  map< uInt,Vector<Bool> > flag;
     176  Bool flagrow;
     177};
     178
     179class BaseMSWriterVisitor: public TableVisitor {
     180  const String *lastFieldName;
     181  uInt lastRecordNo;
     182  uInt lastBeamNo, lastScanNo, lastIfNo;
     183  Int lastSrcType;
     184  uInt lastCycleNo;
     185  Double lastTime;
     186  Int lastPolNo;
     187protected:
     188  const Table &table;
     189  uInt count;
     190public:
     191  BaseMSWriterVisitor(const Table &table)
     192    : table(table)
     193  {
     194    static const String dummy;
     195    lastFieldName = &dummy;
     196    count = 0;
     197  }
     198 
     199  virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
     200  }
     201  virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
     202  }
     203  virtual void enterBeamNo(const uInt recordNo, uInt columnValue) { }
     204  virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) { }
     205  virtual void enterScanNo(const uInt recordNo, uInt columnValue) { }
     206  virtual void leaveScanNo(const uInt recordNo, uInt columnValue) { }
     207  virtual void enterIfNo(const uInt recordNo, uInt columnValue) { }
     208  virtual void leaveIfNo(const uInt recordNo, uInt columnValue) { }
     209  virtual void enterSrcType(const uInt recordNo, Int columnValue) { }
     210  virtual void leaveSrcType(const uInt recordNo, Int columnValue) { }
     211  virtual void enterCycleNo(const uInt recordNo, uInt columnValue) { }
     212  virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) { }
     213  virtual void enterTime(const uInt recordNo, Double columnValue) { }
     214  virtual void leaveTime(const uInt recordNo, Double columnValue) { }
     215  virtual void enterPolNo(const uInt recordNo, uInt columnValue) { }
     216  virtual void leavePolNo(const uInt recordNo, uInt columnValue) { }
     217
     218  virtual Bool visitRecord(const uInt recordNo,
     219                           const String &fieldName,
     220                           const uInt beamNo,
     221                           const uInt scanNo,
     222                           const uInt ifNo,
     223                           const Int srcType,
     224                           const uInt cycleNo,
     225                           const Double time,
     226                           const uInt polNo) { return True ;}
     227
     228  virtual Bool visit(Bool isFirst, const uInt recordNo,
     229                     const uInt nCols, void const *const colValues[]) {
     230    const String *fieldName = NULL;
     231    uInt beamNo, scanNo, ifNo;
     232    Int srcType;
     233    uInt cycleNo;
     234    Double time;
     235    Int polNo;
     236    { // prologue
     237      uInt i = 0;
     238      {
     239        const String *col = (const String*)colValues[i++];
     240        fieldName = &col[recordNo];
     241      }
     242      {
     243        const uInt *col = (const uInt *)colValues[i++];
     244        beamNo = col[recordNo];
     245      }
     246      {
     247        const uInt *col = (const uInt *)colValues[i++];
     248        scanNo = col[recordNo];
     249      }
     250      {
     251        const uInt *col = (const uInt *)colValues[i++];
     252        ifNo = col[recordNo];
     253      }
     254      {
     255        const Int *col = (const Int *)colValues[i++];
     256        srcType = col[recordNo];
     257      }
     258      {
     259        const uInt *col = (const uInt *)colValues[i++];
     260        cycleNo = col[recordNo];
     261      }
     262      {
     263        const Double *col = (const Double *)colValues[i++];
     264        time = col[recordNo];
     265      }
     266      {
     267        const Int *col = (const Int *)colValues[i++];
     268        polNo = col[recordNo];
     269      }
     270      assert(nCols == i);
     271    }
     272
     273    if (isFirst) {
     274      enterFieldName(recordNo, *fieldName);
     275      enterBeamNo(recordNo, beamNo);
     276      enterScanNo(recordNo, scanNo);
     277      enterIfNo(recordNo, ifNo);
     278      enterSrcType(recordNo, srcType);
     279      enterCycleNo(recordNo, cycleNo);
     280      enterTime(recordNo, time);
     281      enterPolNo(recordNo, polNo);
     282    } else {
     283      if (lastFieldName->compare(*fieldName) != 0) {
     284        leavePolNo(lastRecordNo, lastPolNo);
     285        leaveTime(lastRecordNo, lastTime);
     286        leaveCycleNo(lastRecordNo, lastCycleNo);
     287        leaveSrcType(lastRecordNo, lastSrcType);
     288        leaveIfNo(lastRecordNo, lastIfNo);
     289        leaveScanNo(lastRecordNo, lastScanNo);
     290        leaveBeamNo(lastRecordNo, lastBeamNo);
     291        leaveFieldName(lastRecordNo, *lastFieldName);
     292
     293        enterFieldName(recordNo, *fieldName);
     294        enterBeamNo(recordNo, beamNo);
     295        enterScanNo(recordNo, scanNo);
     296        enterIfNo(recordNo, ifNo);
     297        enterSrcType(recordNo, srcType);
     298        enterCycleNo(recordNo, cycleNo);
     299        enterTime(recordNo, time);
     300        enterPolNo(recordNo, polNo);
     301      } else if (lastBeamNo != beamNo) {
     302        leavePolNo(lastRecordNo, lastPolNo);
     303        leaveTime(lastRecordNo, lastTime);
     304        leaveCycleNo(lastRecordNo, lastCycleNo);
     305        leaveSrcType(lastRecordNo, lastSrcType);
     306        leaveIfNo(lastRecordNo, lastIfNo);
     307        leaveScanNo(lastRecordNo, lastScanNo);
     308        leaveBeamNo(lastRecordNo, lastBeamNo);
     309
     310        enterBeamNo(recordNo, beamNo);
     311        enterScanNo(recordNo, scanNo);
     312        enterIfNo(recordNo, ifNo);
     313        enterSrcType(recordNo, srcType);
     314        enterCycleNo(recordNo, cycleNo);
     315        enterTime(recordNo, time);
     316        enterPolNo(recordNo, polNo);
     317      } else if (lastScanNo != scanNo) {
     318        leavePolNo(lastRecordNo, lastPolNo);
     319        leaveTime(lastRecordNo, lastTime);
     320        leaveCycleNo(lastRecordNo, lastCycleNo);
     321        leaveSrcType(lastRecordNo, lastSrcType);
     322        leaveIfNo(lastRecordNo, lastIfNo);
     323        leaveScanNo(lastRecordNo, lastScanNo);
     324
     325        enterScanNo(recordNo, scanNo);
     326        enterIfNo(recordNo, ifNo);
     327        enterSrcType(recordNo, srcType);
     328        enterCycleNo(recordNo, cycleNo);
     329        enterTime(recordNo, time);
     330        enterPolNo(recordNo, polNo);
     331      } else if (lastIfNo != ifNo) {
     332        leavePolNo(lastRecordNo, lastPolNo);
     333        leaveTime(lastRecordNo, lastTime);
     334        leaveCycleNo(lastRecordNo, lastCycleNo);
     335        leaveSrcType(lastRecordNo, lastSrcType);
     336        leaveIfNo(lastRecordNo, lastIfNo);
     337
     338        enterIfNo(recordNo, ifNo);
     339        enterSrcType(recordNo, srcType);
     340        enterCycleNo(recordNo, cycleNo);
     341        enterTime(recordNo, time);
     342        enterPolNo(recordNo, polNo);
     343      } else if (lastSrcType != srcType) {
     344        leavePolNo(lastRecordNo, lastPolNo);
     345        leaveTime(lastRecordNo, lastTime);
     346        leaveCycleNo(lastRecordNo, lastCycleNo);
     347        leaveSrcType(lastRecordNo, lastSrcType);
     348
     349        enterSrcType(recordNo, srcType);
     350        enterCycleNo(recordNo, cycleNo);
     351        enterTime(recordNo, time);
     352        enterPolNo(recordNo, polNo);
     353      } else if (lastCycleNo != cycleNo) {
     354        leavePolNo(lastRecordNo, lastPolNo);
     355        leaveTime(lastRecordNo, lastTime);
     356        leaveCycleNo(lastRecordNo, lastCycleNo);
     357
     358        enterCycleNo(recordNo, cycleNo);
     359        enterTime(recordNo, time);
     360        enterPolNo(recordNo, polNo);
     361      } else if (lastTime != time) {
     362        leavePolNo(lastRecordNo, lastPolNo);
     363        leaveTime(lastRecordNo, lastTime);
     364
     365        enterTime(recordNo, time);
     366        enterPolNo(recordNo, polNo);
     367      } else if (lastPolNo != polNo) {
     368        leavePolNo(lastRecordNo, lastPolNo);
     369        enterPolNo(recordNo, polNo);
     370      }
     371    }
     372    count++;
     373    Bool result = visitRecord(recordNo, *fieldName, beamNo, scanNo, ifNo, srcType,
     374                              cycleNo, time, polNo);
     375
     376    { // epilogue
     377      lastRecordNo = recordNo;
     378
     379      lastFieldName = fieldName;
     380      lastBeamNo = beamNo;
     381      lastScanNo = scanNo;
     382      lastIfNo = ifNo;
     383      lastSrcType = srcType;
     384      lastCycleNo = cycleNo;
     385      lastTime = time;
     386      lastPolNo = polNo;
     387    }
     388    return result ;
     389  }
     390
     391  virtual void finish() {
     392    if (count > 0) {
     393      leavePolNo(lastRecordNo, lastPolNo);
     394      leaveTime(lastRecordNo, lastTime);
     395      leaveCycleNo(lastRecordNo, lastCycleNo);
     396      leaveSrcType(lastRecordNo, lastSrcType);
     397      leaveIfNo(lastRecordNo, lastIfNo);
     398      leaveScanNo(lastRecordNo, lastScanNo);
     399      leaveBeamNo(lastRecordNo, lastBeamNo);
     400      leaveFieldName(lastRecordNo, *lastFieldName);
     401    }
     402  }
     403};
     404
     405class MSWriterVisitor: public BaseMSWriterVisitor, public MSWriterUtils {
     406public:
     407  MSWriterVisitor(const Table &table, Table &mstable)
     408    : BaseMSWriterVisitor(table),
     409      ms(mstable)
     410  {
     411    rowidx = 0 ;
     412    fieldName = "" ;
     413    defaultFieldId = 0 ;
     414    spwId = -1 ;
     415    subscan = 1 ;
     416    ptName = "" ;
     417    srcId = 0 ;
     418   
     419    row = TableRow( ms ) ;
     420
     421    holder.reset() ;
     422
     423    makePolMap() ;
     424    initFrequencies() ;
     425    initTcal() ;
     426
     427    //
     428    // add rows to MS
     429    //
     430    uInt addrow = table.nrow() ;
     431    ms.addRow( addrow ) ;
     432
     433    // attach to Scantable columns
     434    spectraCol.attach( table, "SPECTRA" ) ;
     435    flagtraCol.attach( table, "FLAGTRA" ) ;
     436    flagRowCol.attach( table, "FLAGROW" ) ;
     437    tcalIdCol.attach( table, "TCAL_ID" ) ;
     438    intervalCol.attach( table, "INTERVAL" ) ;
     439    directionCol.attach( table, "DIRECTION" ) ;
     440    scanRateCol.attach( table, "SCANRATE" ) ;
     441    timeCol.attach( table, "TIME" ) ;
     442    freqIdCol.attach( table, "FREQ_ID" ) ;
     443    sourceNameCol.attach( table, "SRCNAME" ) ;
     444    sourceDirectionCol.attach( table, "SRCDIRECTION" ) ;
     445    fieldNameCol.attach( table, "FIELDNAME" ) ;
     446
     447    // MS subtables
     448    attachSubtables() ;
     449
     450    // attach to MS columns
     451    attachMain() ;
     452    attachPointing() ;
     453  }
     454 
     455  virtual void enterFieldName(const uInt recordNo, const String &columnValue) {
     456    //printf("%u: FieldName: %s\n", recordNo, columnValue.c_str());
     457    fieldName = fieldNameCol.asString( recordNo ) ;
     458    String::size_type pos = fieldName.find( "__" ) ;
     459    if ( pos != String::npos ) {
     460      fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
     461      fieldName = fieldName.substr( 0, pos ) ;
     462    }
     463    else {
     464      fieldId = defaultFieldId ;
     465      defaultFieldId++ ;
     466    }
     467    Double tSec = timeCol.asdouble( recordNo ) * 86400.0 ;
     468    Vector<Double> srcDir = sourceDirectionCol( recordNo ) ;
     469    Vector<Double> srate = scanRateCol( recordNo ) ;
     470    String srcName = sourceNameCol.asString( recordNo ) ;
     471
     472    addField( fieldId, fieldName, srcName, srcDir, srate, tSec ) ;
     473
     474    // put value
     475    *fieldIdRF = fieldId ;
     476  }
     477  virtual void leaveFieldName(const uInt recordNo, const String &columnValue) {
     478  }
     479  virtual void enterBeamNo(const uInt recordNo, uInt columnValue) {
     480    //printf("%u: BeamNo: %u\n", recordNo, columnValue);
     481   
     482    feedId = (Int)columnValue ;
     483
     484    // put value
     485    *feed1RF = feedId ;
     486    *feed2RF = feedId ;
     487  }
     488  virtual void leaveBeamNo(const uInt recordNo, uInt columnValue) {
     489  }
     490  virtual void enterScanNo(const uInt recordNo, uInt columnValue) {
     491    //printf("%u: ScanNo: %u\n", recordNo, columnValue);
     492
     493    // put value
     494    // SCAN_NUMBER is 0-based in Scantable while 1-based in MS
     495    *scanNumberRF = (Int)columnValue + 1 ;
     496  }
     497  virtual void leaveScanNo(const uInt recordNo, uInt columnValue) {
     498    subscan = 1 ;
     499  }
     500  virtual void enterIfNo(const uInt recordNo, uInt columnValue) {
     501    //printf("%u: IfNo: %u\n", recordNo, columnValue);
     502
     503    spwId = (Int)columnValue ;
     504    uInt freqId = freqIdCol.asuInt( recordNo ) ;
     505
     506    Vector<Float> sp = spectraCol( recordNo ) ;
     507    uInt nchan = sp.nelements() ;
     508    holder.setNchan( nchan ) ;
     509
     510    addSpectralWindow( spwId, freqId ) ;
     511
     512    addFeed( feedId, spwId ) ;
     513  }
     514  virtual void leaveIfNo(const uInt recordNo, uInt columnValue) {
     515  }
     516  virtual void enterSrcType(const uInt recordNo, Int columnValue) {
     517    //printf("%u: SrcType: %d\n", recordNo, columnValue);
     518
     519    Int stateId = addState( columnValue ) ;
     520
     521    // put value
     522    *stateIdRF = stateId ;
     523  }
     524  virtual void leaveSrcType(const uInt recordNo, Int columnValue) {
     525  }
     526  virtual void enterCycleNo(const uInt recordNo, uInt columnValue) {
     527    //printf("%u: CycleNo: %u\n", recordNo, columnValue);
     528  }
     529  virtual void leaveCycleNo(const uInt recordNo, uInt columnValue) {
     530  }
     531  virtual void enterTime(const uInt recordNo, Double columnValue) {
     532    //printf("%u: Time: %f\n", recordNo, columnValue);
     533
     534    Double timeSec = columnValue * 86400.0 ;
     535    Double interval = intervalCol.asdouble( recordNo ) ;
     536
     537    if ( ptName.empty() ) {
     538      Vector<Double> dir = directionCol( recordNo ) ;
     539      Vector<Double> rate = scanRateCol( recordNo ) ;
     540      if ( anyNE( rate, 0.0 ) ) {
     541        Matrix<Double> msdir( 2, 2 ) ;
     542        msdir.column( 0 ) = dir ;
     543        msdir.column( 1 ) = rate ;
     544        addPointing( timeSec, interval, msdir ) ;
     545      }
     546      else {
     547        Matrix<Double> msdir( 2, 1 ) ;
     548        msdir.column( 0 ) = dir ;
     549        addPointing( timeSec, interval, msdir ) ;
     550      }
     551    }
     552
     553    // put value
     554    *timeRF = timeSec ;
     555    *timeCentroidRF = timeSec ;
     556    *intervalRF = interval ;
     557    *exposureRF = interval ;
     558  }
     559  virtual void leaveTime(const uInt recordNo, Double columnValue) {
     560    if ( holder.nPol() > 0 ) {
     561      Vector<Float> w = holder.getWeight() ;
     562      Cube<Bool> c = holder.getFlagCategory() ;
     563      Bool flr = holder.getFlagRow() ;
     564      Matrix<Bool> fl = holder.getFlag() ;
     565      Vector<uInt> polnos = holder.polNos() ;
     566      Int polId = addPolarization( polnos ) ;
     567      Int ddId = addDataDescription( polId, spwId ) ;
     568       
     569      // put field
     570      *dataDescIdRF = ddId ;
     571      *flagRowRF = flr ;
     572      weightRF.define( w ) ;
     573      sigmaRF.define( w ) ;
     574      flagCategoryRF.define( c ) ;
     575      flagRF.define( fl ) ;
     576      if ( useFloat ) {
     577        Matrix<Float> sp = holder.getData() ;
     578        floatDataRF.define( sp ) ;
     579      }
     580      else {
     581        Matrix<Complex> sp = holder.getComplexData() ;
     582        dataRF.define( sp ) ;
     583      }
     584     
     585      // commit row
     586      row.put( rowidx ) ;
     587      rowidx++ ;
     588
     589      // reset holder
     590      holder.reset() ;
     591    }
     592    if ( tcalKey != -1 ) {
     593      tcalNotYet[tcalKey] = False ;
     594      tcalKey = -1 ;
     595    }
     596  }
     597  virtual void enterPolNo(const uInt recordNo, uInt columnValue) {
     598    //printf("%u: PolNo: %d\n", recordNo, columnValue);
     599    uInt tcalId = tcalIdCol.asuInt( recordNo ) ;
     600    if ( tcalKey == -1 ) {
     601      tcalKey = tcalId ;
     602    }
     603    if ( tcalNotYet[tcalKey] ) {
     604      map< Int,Vector<uInt> >::iterator itr = tcalIdRec.find( tcalKey ) ;
     605      if ( itr != tcalIdRec.end() ) {
     606        Vector<uInt> ids = itr->second ;
     607        uInt nrow = ids.nelements() ;
     608        ids.resize( nrow+1, True ) ;
     609        ids[nrow] = tcalId ;
     610        tcalIdRec.erase( tcalKey ) ;
     611        tcalIdRec[tcalKey] = ids ;
     612      }
     613      else {
     614        Vector<uInt> rows( 1, tcalId ) ;
     615        tcalIdRec[tcalKey] = rows ;
     616      }
     617    }
     618    map< Int,Vector<uInt> >::iterator itr = tcalRowRec.find( tcalKey ) ;
     619    if ( itr != tcalRowRec.end() ) {
     620      Vector<uInt> rows = itr->second ;
     621      uInt nrow = rows.nelements() ;
     622      rows.resize( nrow+1, True ) ;
     623      rows[nrow] = recordNo ;
     624      tcalRowRec.erase( tcalKey ) ;
     625      tcalRowRec[tcalKey] = rows ;
     626    }
     627    else {
     628      Vector<uInt> rows( 1, recordNo ) ;
     629      tcalRowRec[tcalKey] = rows ;
     630    }
     631  }
     632  virtual void leavePolNo(const uInt recordNo, uInt columnValue) {
     633  }
     634
     635  virtual Bool visitRecord(const uInt recordNo,
     636                           const String &fieldName,
     637                           const uInt beamNo,
     638                           const uInt scanNo,
     639                           const uInt ifNo,
     640                           const Int srcType,
     641                           const uInt cycleNo,
     642                           const Double time,
     643                           const uInt polNo) {
     644    //printf("%u: %s, %u, %u, %u, %d, %u, %f, %d\n", recordNo,
     645    //       fieldName.c_str(), beamNo, scanNo, ifNo, srcType, cycleNo, time, polNo);
     646
     647    Vector<Float> sp = spectraCol( recordNo ) ;
     648    Vector<uChar> tmp = flagtraCol( recordNo ) ;
     649    Vector<Bool> fl( tmp.shape() ) ;
     650    convertArray( fl, tmp ) ;
     651    Bool flr = (Bool)flagRowCol.asuInt( recordNo ) ;
     652    holder.accumulate( polNo, sp, fl, flr ) ;
     653
     654    return True ;
     655  }
     656
     657  virtual void finish() {
     658    BaseMSWriterVisitor::finish();
     659    //printf("Total: %u\n", count);
     660
     661    // remove rows
     662    if ( ms.nrow() > rowidx ) {
     663      uInt numRemove = ms.nrow() - rowidx ;
     664      //cout << "numRemove = " << numRemove << endl ;
     665      Vector<uInt> rows( numRemove ) ;
     666      indgen( rows, rowidx ) ;
     667      ms.removeRow( rows ) ;
     668    }
     669
     670    // fill empty SPECTRAL_WINDOW rows
     671    infillSpectralWindow() ;
     672  }
     673
     674  void dataColumnName( String name )
     675  {
     676    TableRecord &r = row.record() ;
     677    if ( name == "DATA" ) {
     678      useFloat = False ;
     679      dataRF.attachToRecord( r, name ) ;
     680    }
     681    else if ( name == "FLOAT_DATA" ) {
     682      useFloat = True ;
     683      floatDataRF.attachToRecord( r, name ) ;
     684    }
     685  }
     686  void pointingTableName( String name ) {
     687    ptName = name ;
     688  }
     689  void setSourceRecord( Record &r ) {
     690    srcRec = r ;
     691  }
     692  map< Int,Vector<uInt> > &getTcalIdRecord() { return tcalIdRec ; }
     693  map< Int,Vector<uInt> > &getTcalRowRecord() { return tcalRowRec ; }
     694private:
     695  void addField( Int &fid, String &fname, String &srcName,
     696                 Vector<Double> &sdir, Vector<Double> &srate,
     697                 Double &tSec )
     698  {
     699    uInt nrow = fieldtab.nrow() ;
     700    while( (Int)nrow <= fid ) {
     701      fieldtab.addRow( 1, True ) ;
     702      nrow++ ;
     703    }
     704
     705    Matrix<Double> dir ;
     706    Int numPoly = 0 ;
     707    if ( anyNE( srate, 0.0 ) ) {
     708      dir.resize( 2, 2 ) ;
     709      dir.column( 0 ) = sdir ;
     710      dir.column( 1 ) = srate ;
     711      numPoly = 1 ;
     712    }
     713    else {
     714      dir.resize( 2, 1 ) ;
     715      dir.column( 0 ) = sdir ;
     716    }
     717    srcId = srcRec.asInt( srcName ) ;
     718
     719    TableRow tr( fieldtab ) ;
     720    TableRecord &r = tr.record() ;
     721    putField( "NAME", r, fname ) ;
     722    putField( "NUM_POLY", r, numPoly ) ;
     723    putField( "TIME", r, tSec ) ;
     724    putField( "SOURCE_ID", r, srcId ) ;
     725    defineField( "DELAY_DIR", r, dir ) ;
     726    defineField( "REFERENCE_DIR", r, dir ) ;
     727    defineField( "PHASE_DIR", r, dir ) ;
     728    tr.put( fid ) ;
     729
     730    // for POINTING table
     731    *poNameRF = fname ;
     732  }
     733  Int addState( Int &id )
     734  {
     735    String obsMode ;
     736    Bool isSignal ;
     737    Double tnoise ;
     738    Double tload ;
     739    queryType( id, obsMode, isSignal, tnoise, tload ) ;
     740
     741    String key = obsMode+"_"+String::toString( subscan ) ;
     742    Int idx = -1 ;
     743    uInt nEntry = stateEntry.nelements() ;
     744    for ( uInt i = 0 ; i < nEntry ; i++ ) {
     745      if ( stateEntry[i] == key ) {
     746        idx = i ;
     747        break ;
     748      }
     749    }
     750    if ( idx == -1 ) {
     751      uInt nrow = statetab.nrow() ;
     752      statetab.addRow( 1, True ) ;
     753      TableRow tr( statetab ) ;
     754      TableRecord &r = tr.record() ;
     755      putField( "OBS_MODE", r, obsMode ) ;
     756      putField( "SIG", r, isSignal ) ;
     757      isSignal = !isSignal ;
     758      putField( "REF", r, isSignal ) ;
     759      putField( "CAL", r, tnoise ) ;
     760      putField( "LOAD", r, tload ) ;
     761      tr.put( nrow ) ;
     762      idx = nrow ;
     763
     764      stateEntry.resize( nEntry+1, True ) ;
     765      stateEntry[nEntry] = key ;
     766    }
     767    subscan++ ;
     768
     769    return idx ;
     770  }
     771  void addPointing( Double &tSec, Double &interval, Matrix<Double> &dir )
     772  {
     773    uInt nrow = potab.nrow() ;
     774    potab.addRow( 1, True ) ;
     775
     776    *poNumPolyRF = dir.ncolumn() - 1 ;
     777    *poTimeRF = tSec ;
     778    *poTimeOriginRF = tSec ;
     779    *poIntervalRF = interval ;
     780    poDirectionRF.define( dir ) ;
     781    poTargetRF.define( dir ) ;
     782    porow.put( nrow ) ;
     783  }
     784  Int addPolarization( Vector<uInt> &nos )
     785  {
     786    Int idx = -1 ;
     787    uInt nEntry = polEntry.size() ;
     788    for ( uInt i = 0 ; i < nEntry ; i++ ) {
     789      if ( polEntry[i].conform( nos ) && allEQ( polEntry[i], nos ) ) {
     790        idx = i ;
     791        break ;
     792      }
     793    }
     794   
     795    Int numCorr ;
     796    Vector<Int> corrType ;
     797    Matrix<Int> corrProduct ;
     798    polProperty( nos, numCorr, corrType, corrProduct ) ;
     799
     800    if ( idx == -1 ) {
     801      uInt nrow = poltab.nrow() ;
     802      poltab.addRow( 1, True ) ;
     803      TableRow tr( poltab ) ;
     804      TableRecord &r = tr.record() ;
     805      putField( "NUM_CORR", r, numCorr ) ;
     806      defineField( "CORR_TYPE", r, corrType ) ;
     807      defineField( "CORR_PRODUCT", r, corrProduct ) ;
     808      tr.put( nrow ) ;
     809      idx = nrow ;
     810
     811      polEntry.resize( nEntry+1 ) ;
     812      polEntry[nEntry] = nos ;
     813    }
     814
     815    return idx ;
     816  }
     817  Int addDataDescription( Int pid, Int sid )
     818  {
     819    Int idx = -1 ;
     820    uInt nEntry = ddEntry.nrow() ;
     821    Vector<Int> key( 2 ) ;
     822    key[0] = pid ;
     823    key[1] = sid ;
     824    for ( uInt i = 0 ; i < nEntry ; i++ ) {
     825      if ( allEQ( ddEntry.row(i), key ) ) {
     826        idx = i ;
     827        break ;
     828      }
     829    }
     830
     831    if ( idx == -1 ) {
     832      uInt nrow = ddtab.nrow() ;
     833      ddtab.addRow( 1, True ) ;
     834      TableRow tr( ddtab ) ;
     835      TableRecord &r = tr.record() ;
     836      putField( "POLARIZATION_ID", r, pid ) ;
     837      putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
     838      tr.put( nrow ) ;
     839      idx = nrow ;
     840
     841      ddEntry.resize( nEntry+1, 2, True ) ;
     842      ddEntry.row(nEntry) = key ;
     843    }
     844
     845    return idx ;
     846  }
     847  void infillSpectralWindow()
     848  {
     849    ROScalarColumn<Int> nchanCol( spwtab, "NUM_CHAN" ) ;
     850    Vector<Int> nchan = nchanCol.getColumn() ;
     851    TableRow tr( spwtab ) ;
     852    TableRecord &r = tr.record() ;
     853    Int mfr = 1 ;
     854    Vector<Double> dummy( 1, 0.0 ) ;
     855    putField( "MEAS_FREQ_REF", r, mfr ) ;
     856    defineField( "CHAN_FREQ", r, dummy ) ;
     857    defineField( "CHAN_WIDTH", r, dummy ) ;
     858    defineField( "EFFECTIVE_BW", r, dummy ) ;
     859    defineField( "RESOLUTION", r, dummy ) ;
     860
     861    for ( uInt i = 0 ; i < spwtab.nrow() ; i++ ) {
     862      if ( nchan[i] == 0 )
     863        tr.put( i ) ;
     864    }
     865  }
     866  void addSpectralWindow( Int sid, uInt fid )
     867  {
     868    if ( !processedFreqId[fid] ) {
     869      uInt nrow = spwtab.nrow() ;
     870      while( (Int)nrow <= sid ) {
     871        spwtab.addRow( 1, True ) ;
     872        nrow++ ;
     873      }
     874      processedFreqId[fid] = True ;
     875    }
     876
     877    Double rp = refpix[fid] ;
     878    Double rv = refval[fid] ;
     879    Double ic = increment[fid] ;
     880
     881    Int mfrInt = (Int)freqframe ;
     882    Int nchan = holder.nChan() ;
     883    Double bw = nchan * abs( ic ) ;
     884    Double reffreq = rv - rp * ic ;
     885    Int netsb = 0 ; // USB->0, LSB->1
     886    if ( ic < 0 )
     887      netsb = 1 ;
     888    Vector<Double> res( nchan, abs(ic) ) ;
     889    Vector<Double> chanf( nchan ) ;
     890    indgen( chanf, reffreq, ic ) ;
     891
     892    TableRow tr( spwtab ) ;
     893    TableRecord &r = tr.record() ;
     894    putField( "MEAS_FREQ_REF", r, mfrInt ) ;
     895    putField( "NUM_CHAN", r, nchan ) ;
     896    putField( "TOTAL_BANDWIDTH", r, bw ) ;
     897    putField( "REF_FREQUENCY", r, reffreq ) ;
     898    putField( "NET_SIDEBAND", r, netsb ) ;
     899    defineField( "RESOLUTION", r, res ) ;
     900    defineField( "CHAN_WIDTH", r, res ) ;
     901    defineField( "EFFECTIVE_BW", r, res ) ;
     902    defineField( "CHAN_FREQ", r, chanf ) ;
     903    tr.put( sid ) ;
     904  }
     905  void addFeed( Int fid, Int sid )
     906  {
     907    Int idx = -1 ;
     908    uInt nEntry = feedEntry.nrow() ;
     909    Vector<Int> key( 2 ) ;
     910    key[0] = fid ;
     911    key[1] = sid ;
     912    for ( uInt i = 0 ; i < nEntry ; i++ ) {
     913      if ( allEQ( feedEntry.row(i), key ) ) {
     914        idx = i ;
     915        break ;
     916      }
     917    }
     918
     919    if ( idx == -1 ) {
     920      uInt nrow = feedtab.nrow() ;
     921      feedtab.addRow( 1, True ) ;
     922      Int numReceptors = 2 ;
     923      Vector<String> polType( numReceptors ) ;
     924      Matrix<Double> beamOffset( 2, numReceptors, 0.0 ) ;
     925      Vector<Double> receptorAngle( numReceptors, 0.0 ) ;
     926      if ( poltype == "linear" ) {
     927        polType[0] = "X" ;
     928        polType[1] = "Y" ;
     929      }
     930      else if ( poltype == "circular" ) {
     931        polType[0] = "R" ;
     932        polType[1] = "L" ;
     933      }
     934      else {
     935        polType[0] = "X" ;
     936        polType[1] = "Y" ;
     937      }
     938      Matrix<Complex> polResponse( numReceptors, numReceptors, 0.0 ) ;
     939     
     940      TableRow tr( feedtab ) ;
     941      TableRecord &r = tr.record() ;
     942      putField( "FEED_ID", r, fid ) ;
     943      putField( "BEAM_ID", r, fid ) ;
     944      Int tmp = 0 ;
     945      putField( "ANTENNA_ID", r, tmp ) ;
     946      putField( "SPECTRAL_WINDOW_ID", r, sid ) ;
     947      putField( "NUM_RECEPTORS", r, numReceptors ) ;
     948      defineField( "POLARIZATION_TYPE", r, polType ) ;
     949      defineField( "BEAM_OFFSET", r, beamOffset ) ;
     950      defineField( "RECEPTOR_ANGLE", r, receptorAngle ) ;
     951      defineField( "POL_RESPONSE", r, polResponse ) ;
     952      tr.put( nrow ) ;
     953
     954      feedEntry.resize( nEntry+1, 2, True ) ;
     955      feedEntry.row( nEntry ) = key ;
     956    }
     957  }
     958  void makePolMap()
     959  {
     960    const TableRecord &keys = table.keywordSet() ;
     961    poltype = keys.asString( "POLTYPE" ) ;
     962
     963    if ( poltype == "stokes" ) {
     964      polmap.resize( 4 ) ;
     965      polmap[0] = Stokes::I ;
     966      polmap[1] = Stokes::Q ;
     967      polmap[2] = Stokes::U ;
     968      polmap[3] = Stokes::V ;
     969    }
     970    else if ( poltype == "linear" ) {
     971      polmap.resize( 4 ) ;
     972      polmap[0] = Stokes::XX ;
     973      polmap[1] = Stokes::YY ;
     974      polmap[2] = Stokes::XY ;
     975      polmap[3] = Stokes::YX ;
     976    }
     977    else if ( poltype == "circular" ) {
     978      polmap.resize( 4 ) ;
     979      polmap[0] = Stokes::RR ;
     980      polmap[1] = Stokes::LL ;
     981      polmap[2] = Stokes::RL ;
     982      polmap[3] = Stokes::LR ;
     983    }
     984    else if ( poltype == "linpol" ) {
     985      polmap.resize( 2 ) ;
     986      polmap[0] = Stokes::Plinear ;
     987      polmap[1] = Stokes::Pangle ;
     988    }
     989    else {
     990      polmap.resize( 0 ) ;
     991    }
     992  }
     993  void initFrequencies()
     994  {
     995    const TableRecord &keys = table.keywordSet() ;
     996    Table tab = keys.asTable( "FREQUENCIES" ) ;
     997    ROScalarColumn<uInt> idcol( tab, "ID" ) ;
     998    ROScalarColumn<Double> rpcol( tab, "REFPIX" ) ;
     999    ROScalarColumn<Double> rvcol( tab, "REFVAL" ) ;
     1000    ROScalarColumn<Double> iccol( tab, "INCREMENT" ) ;
     1001    Vector<uInt> id = idcol.getColumn() ;
     1002    Vector<Double> rp = rpcol.getColumn() ;
     1003    Vector<Double> rv = rvcol.getColumn() ;
     1004    Vector<Double> ic = iccol.getColumn() ;
     1005    for ( uInt i = 0 ; i < id.nelements() ; i++ ) {
     1006      processedFreqId.insert( pair<uInt,Bool>( id[i], False ) ) ;
     1007      refpix.insert( pair<uInt,Double>( id[i], rp[i] ) ) ;
     1008      refval.insert( pair<uInt,Double>( id[i], rv[i] ) ) ;
     1009      increment.insert( pair<uInt,Double>( id[i], ic[i] ) ) ;
     1010    }
     1011    String frameStr = tab.keywordSet().asString( "BASEFRAME" ) ;
     1012    MFrequency::getType( freqframe, frameStr ) ;
     1013  }
     1014  void attachSubtables()
     1015  {
     1016    const TableRecord &keys = table.keywordSet() ;
     1017    TableRecord &mskeys = ms.rwKeywordSet() ;
     1018
     1019    // FIELD table
     1020    fieldtab = mskeys.asTable( "FIELD" ) ;
     1021
     1022    // SPECTRAL_WINDOW table
     1023    spwtab = mskeys.asTable( "SPECTRAL_WINDOW" ) ;
     1024
     1025    // POINTING table
     1026    potab = mskeys.asTable( "POINTING" ) ;
     1027
     1028    // POLARIZATION table
     1029    poltab = mskeys.asTable( "POLARIZATION" ) ;
     1030
     1031    // DATA_DESCRIPTION table
     1032    ddtab = mskeys.asTable( "DATA_DESCRIPTION" ) ;
     1033
     1034    // STATE table
     1035    statetab = mskeys.asTable( "STATE" ) ;
     1036
     1037    // FEED table
     1038    feedtab = mskeys.asTable( "FEED" ) ;
     1039  }
     1040  void attachMain()
     1041  {
     1042    TableRecord &r = row.record() ;
     1043    dataDescIdRF.attachToRecord( r, "DATA_DESC_ID" ) ;
     1044    flagRowRF.attachToRecord( r, "FLAG_ROW" ) ;
     1045    weightRF.attachToRecord( r, "WEIGHT" ) ;
     1046    sigmaRF.attachToRecord( r, "SIGMA" ) ;
     1047    flagCategoryRF.attachToRecord( r, "FLAG_CATEGORY" ) ;
     1048    flagRF.attachToRecord( r, "FLAG" ) ;
     1049    timeRF.attachToRecord( r, "TIME" ) ;
     1050    timeCentroidRF.attachToRecord( r, "TIME_CENTROID" ) ;
     1051    intervalRF.attachToRecord( r, "INTERVAL" ) ;
     1052    exposureRF.attachToRecord( r, "EXPOSURE" ) ;
     1053    fieldIdRF.attachToRecord( r, "FIELD_ID" ) ;
     1054    feed1RF.attachToRecord( r, "FEED1" ) ;
     1055    feed2RF.attachToRecord( r, "FEED2" ) ;
     1056    scanNumberRF.attachToRecord( r, "SCAN_NUMBER" ) ;
     1057    stateIdRF.attachToRecord( r, "STATE_ID" ) ;
     1058
     1059    // constant values
     1060    Int id = 0 ;
     1061    RecordFieldPtr<Int> intRF( r, "OBSERVATION_ID" ) ;
     1062    *intRF = 0 ;
     1063    intRF.attachToRecord( r, "ANTENNA1" ) ;
     1064    *intRF = 0 ;
     1065    intRF.attachToRecord( r, "ANTENNA2" ) ;
     1066    *intRF = 0 ;
     1067    intRF.attachToRecord( r, "ARRAY_ID" ) ;
     1068    *intRF = 0 ;
     1069    intRF.attachToRecord( r, "PROCESSOR_ID" ) ;
     1070    *intRF = 0 ;
     1071    RecordFieldPtr< Vector<Double> > arrayRF( r, "UVW" ) ;
     1072    arrayRF.define( Vector<Double>( 3, 0.0 ) ) ;
     1073  }
     1074  void attachPointing()
     1075  {
     1076    porow = TableRow( potab ) ;
     1077    TableRecord &r = porow.record() ;
     1078    poNumPolyRF.attachToRecord( r, "NUM_POLY" ) ;
     1079    poTimeRF.attachToRecord( r, "TIME" ) ;
     1080    poTimeOriginRF.attachToRecord( r, "TIME_ORIGIN" ) ;
     1081    poIntervalRF.attachToRecord( r, "INTERVAL" ) ;
     1082    poNameRF.attachToRecord( r, "NAME" ) ;
     1083    poDirectionRF.attachToRecord( r, "DIRECTION" ) ;
     1084    poTargetRF.attachToRecord( r, "TARGET" ) ;
     1085   
     1086    // constant values
     1087    RecordFieldPtr<Int> antIdRF( r, "ANTENNA_ID" ) ;
     1088    *antIdRF = 0 ;
     1089    RecordFieldPtr<Bool> trackingRF( r, "TRACKING" ) ;
     1090    *trackingRF = True ;
     1091  }
     1092  void queryType( Int type, String &stype, Bool &b, Double &t, Double &l )
     1093  {
     1094    t = 0.0 ;
     1095    l = 0.0 ;
     1096
     1097    String sep1="#" ;
     1098    String sep2="," ;
     1099    String target="OBSERVE_TARGET" ;
     1100    String atmcal="CALIBRATE_TEMPERATURE" ;
     1101    String onstr="ON_SOURCE" ;
     1102    String offstr="OFF_SOURCE" ;
     1103    String pswitch="POSITION_SWITCH" ;
     1104    String nod="NOD" ;
     1105    String fswitch="FREQUENCY_SWITCH" ;
     1106    String sigstr="SIG" ;
     1107    String refstr="REF" ;
     1108    String unspecified="UNSPECIFIED" ;
     1109    String ftlow="LOWER" ;
     1110    String fthigh="HIGHER" ;
     1111    switch ( type ) {
     1112    case SrcType::PSON:
     1113      stype = target+sep1+onstr+sep2+pswitch ;
     1114      b = True ;
     1115      break ;
     1116    case SrcType::PSOFF:
     1117      stype = target+sep1+offstr+sep2+pswitch ;
     1118      b = False ;
     1119      break ;
     1120    case SrcType::NOD:
     1121      stype = target+sep1+onstr+sep2+nod ;
     1122      b = True ;
     1123      break ;
     1124    case SrcType::FSON:
     1125      stype = target+sep1+onstr+sep2+fswitch+sep1+sigstr ;
     1126      b = True ;
     1127      break ;
     1128    case SrcType::FSOFF:
     1129      stype = target+sep1+onstr+sep2+fswitch+sep1+refstr ;
     1130      b = False ;
     1131      break ;
     1132    case SrcType::SKY:
     1133      stype = atmcal+sep1+offstr+sep2+unspecified ;
     1134      b = False ;
     1135      break ;
     1136    case SrcType::HOT:
     1137      stype = atmcal+sep1+offstr+sep2+unspecified ;
     1138      b = False ;
     1139      break ;
     1140    case SrcType::WARM:
     1141      stype = atmcal+sep1+offstr+sep2+unspecified ;
     1142      b = False ;
     1143      break ;
     1144    case SrcType::COLD:
     1145      stype = atmcal+sep1+offstr+sep2+unspecified ;
     1146      b = False ;
     1147      break ;
     1148    case SrcType::PONCAL:
     1149      stype = atmcal+sep1+onstr+sep2+pswitch ;
     1150      b = True ;
     1151      break ;
     1152    case SrcType::POFFCAL:
     1153      stype = atmcal+sep1+offstr+sep2+pswitch ;
     1154      b = False ;
     1155      break ;
     1156    case SrcType::NODCAL:
     1157      stype = atmcal+sep1+onstr+sep2+nod ;
     1158      b = True ;
     1159      break ;
     1160    case SrcType::FONCAL:
     1161      stype = atmcal+sep1+onstr+sep2+fswitch+sep1+sigstr ;
     1162      b = True ;
     1163      break ;
     1164    case SrcType::FOFFCAL:
     1165      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+refstr ;
     1166      b = False ;
     1167      break ;
     1168    case SrcType::FSLO:
     1169      stype = target+sep1+onstr+sep2+fswitch+sep1+ftlow ;
     1170      b = True ;
     1171      break ;
     1172    case SrcType::FLOOFF:
     1173      stype = target+sep1+offstr+sep2+fswitch+sep1+ftlow ;
     1174      b = False ;
     1175      break ;
     1176    case SrcType::FLOSKY:
     1177      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
     1178      b = False ;
     1179      break ;
     1180    case SrcType::FLOHOT:
     1181      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
     1182      b = False ;
     1183      break ;
     1184    case SrcType::FLOWARM:
     1185      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
     1186      b = False ;
     1187      break ;
     1188    case SrcType::FLOCOLD:
     1189      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+ftlow ;
     1190      b = False ;
     1191      break ;
     1192    case SrcType::FSHI:
     1193      stype = target+sep1+onstr+sep2+fswitch+sep1+fthigh ;
     1194      b = True ;
     1195      break ;
     1196    case SrcType::FHIOFF:
     1197      stype = target+sep1+offstr+sep2+fswitch+sep1+fthigh ;
     1198      b = False ;
     1199      break ;
     1200    case SrcType::FHISKY:
     1201      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
     1202      b = False ;
     1203      break ;
     1204    case SrcType::FHIHOT:
     1205      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
     1206      b = False ;
     1207      break ;
     1208    case SrcType::FHIWARM:
     1209      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
     1210      b = False ;
     1211      break ;
     1212    case SrcType::FHICOLD:
     1213      stype = atmcal+sep1+offstr+sep2+fswitch+sep1+fthigh ;
     1214      b = False ;
     1215      break ;
     1216    case SrcType::SIG:
     1217      stype = target+sep1+onstr+sep2+unspecified ;
     1218      b = True ;
     1219      break ;
     1220    case SrcType::REF:
     1221      stype = target+sep1+offstr+sep2+unspecified ;
     1222      b = False ;
     1223      break ;
     1224    default:
     1225      stype = unspecified ;
     1226      b = True ;
     1227      break ;
     1228    }
     1229  }
     1230  void polProperty( Vector<uInt> &nos, Int &n, Vector<Int> &c, Matrix<Int> &cp )
     1231  {
     1232    n = nos.nelements() ;
     1233    c.resize( n ) ;
     1234    for ( Int i = 0 ; i < n ; i++ )
     1235      c[i] = (Int)polmap[nos[i]] ;
     1236    cp.resize( 2, n ) ;
     1237    if ( n == 1 )
     1238      cp = 0 ;
     1239    else if ( n == 2 ) {
     1240      cp.column( 0 ) = 0 ;
     1241      cp.column( 1 ) = 1 ;
     1242    }
     1243    else {
     1244      cp.column( 0 ) = 0 ;
     1245      cp.column( 1 ) = 1 ;
     1246      cp( 0, 1 ) = 0 ;
     1247      cp( 1, 1 ) = 1 ;
     1248      cp( 0, 2 ) = 1 ;
     1249      cp( 1, 2 ) = 0 ;
     1250    }
     1251  }
     1252  void initTcal()
     1253  {
     1254    const TableRecord &rec = table.keywordSet() ;
     1255    Table tcalTable = rec.asTable( "TCAL" ) ;
     1256    ROScalarColumn<uInt> idCol( tcalTable, "ID" ) ;
     1257    Vector<uInt> id = idCol.getColumn() ;
     1258    uInt maxId = max( id ) ;
     1259    tcalNotYet.resize( maxId+1 ) ;
     1260    tcalNotYet = True ;
     1261    tcalKey = -1 ;
     1262  }
     1263
     1264  Table &ms;
     1265  TableRow row;
     1266  uInt rowidx;
     1267  String fieldName;
     1268  Int fieldId;
     1269  Int srcId;
     1270  Int defaultFieldId;
     1271  Int spwId;
     1272  Int feedId;
     1273  Int subscan;
     1274  PolarizedComponentHolder holder;
     1275  String ptName;
     1276  Bool useFloat;
     1277  String poltype;
     1278  Vector<Stokes::StokesTypes> polmap;
     1279
     1280  // MS subtables
     1281  Table spwtab;
     1282  Table statetab;
     1283  Table ddtab;
     1284  Table poltab;
     1285  Table fieldtab;
     1286  Table feedtab;
     1287  Table potab;
     1288
     1289  // Scantable MAIN columns
     1290  ROArrayColumn<Float> spectraCol;
     1291  ROArrayColumn<Double> directionCol,scanRateCol,sourceDirectionCol;
     1292  ROArrayColumn<uChar> flagtraCol;
     1293  ROTableColumn tcalIdCol,intervalCol,flagRowCol,timeCol,freqIdCol,
     1294    sourceNameCol,fieldNameCol;
     1295
     1296  // MS MAIN columns
     1297  RecordFieldPtr<Int> dataDescIdRF,fieldIdRF,feed1RF,feed2RF,
     1298    scanNumberRF,stateIdRF;
     1299  RecordFieldPtr<Bool> flagRowRF;
     1300  RecordFieldPtr<Double> timeRF,timeCentroidRF,intervalRF,exposureRF;
     1301  RecordFieldPtr< Vector<Float> > weightRF,sigmaRF;
     1302  RecordFieldPtr< Cube<Bool> > flagCategoryRF;
     1303  RecordFieldPtr< Matrix<Bool> > flagRF;
     1304  RecordFieldPtr< Matrix<Float> > floatDataRF;
     1305  RecordFieldPtr< Matrix<Complex> > dataRF;
     1306
     1307  // MS POINTING columns
     1308  TableRow porow;
     1309  RecordFieldPtr<Int> poNumPolyRF ;
     1310  RecordFieldPtr<Double> poTimeRF,
     1311    poTimeOriginRF,
     1312    poIntervalRF ;
     1313  RecordFieldPtr<String> poNameRF ;
     1314  RecordFieldPtr< Matrix<Double> > poDirectionRF,
     1315    poTargetRF ;
     1316
     1317  Vector<String> stateEntry;
     1318  Matrix<Int> ddEntry;
     1319  Matrix<Int> feedEntry;
     1320  vector< Vector<uInt> > polEntry;
     1321  map<uInt,Bool> processedFreqId;
     1322  map<uInt,Double> refpix;
     1323  map<uInt,Double> refval;
     1324  map<uInt,Double> increment;
     1325  MFrequency::Types freqframe;
     1326  Record srcRec;
     1327  map< Int,Vector<uInt> > tcalIdRec;
     1328  map< Int,Vector<uInt> > tcalRowRec;
     1329  Int tcalKey;
     1330  Vector<Bool> tcalNotYet;
     1331};
    561332
    571333MSWriter::MSWriter(CountedPtr<Scantable> stable)
     
    791355    delete mstable_ ;
    801356}
    81  
     1357
    821358bool MSWriter::write(const string& filename, const Record& rec)
    831359{
    841360  os_.origin( LogOrigin( "MSWriter", "write()", WHERE ) ) ;
    85 //   double startSec = mathutil::gettimeofday_sec() ;
    86 //   os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
     1361  //double startSec = mathutil::gettimeofday_sec() ;
     1362  //os_ << "start MSWriter::write() startSec=" << startSec << LogIO::POST ;
    871363
    881364  filename_ = filename ;
     
    1341410    fillWeather() ;
    1351411
    136   // MAIN
    137   // Iterate over several ids
    138   Vector<uInt> processedFreqId( 0 ) ;
    139   Int defaultFieldId = 0 ;
    140 
    141   // row based
    142   TableRow row( *mstable_ ) ;
    143   TableRecord &trec = row.record() ;
    144   NoticeTarget *dataRF = 0 ;
    145   if ( useFloatData_ )
    146     dataRF = new RecordFieldPtr< Array<Float> >( trec, "FLOAT_DATA" ) ;
    147   else if ( useData_ )
    148     dataRF = new RecordFieldPtr< Array<Complex> >( trec, "DATA" ) ;
    149   RecordFieldPtr< Array<Bool> > flagRF( trec, "FLAG" ) ;
    150   RecordFieldPtr<Bool> flagrowRF( trec, "FLAG_ROW" ) ;
    151   RecordFieldPtr<Double> timeRF( trec, "TIME" ) ;
    152   RecordFieldPtr<Double> timecRF( trec, "TIME_CENTROID" ) ;
    153   RecordFieldPtr<Double> intervalRF( trec, "INTERVAL" ) ;
    154   RecordFieldPtr<Double> exposureRF( trec, "EXPOSURE" ) ;
    155   RecordFieldPtr< Array<Float> > weightRF( trec, "WEIGHT" ) ;
    156   RecordFieldPtr< Array<Float> > sigmaRF( trec, "SIGMA" ) ;
    157   RecordFieldPtr<Int> ddidRF( trec, "DATA_DESC_ID" ) ;
    158   RecordFieldPtr<Int> stateidRF( trec, "STATE_ID" ) ;
    159   RecordFieldPtr< Array<Bool> > flagcatRF( trec, "FLAG_CATEGORY" ) ;
    160 
    161   // OBSERVATION_ID is always 0
    162   RecordFieldPtr<Int> intRF( trec, "OBSERVATION_ID" ) ;
    163   *intRF = 0 ;
    164  
    165   // ANTENNA1 and ANTENNA2 are always 0
    166   intRF.attachToRecord( trec, "ANTENNA1" ) ;
    167   *intRF = 0 ;
    168   intRF.attachToRecord( trec, "ANTENNA2" ) ;
    169   *intRF = 0 ;
    170  
    171   // ARRAY_ID is tentatively set to 0
    172   intRF.attachToRecord( trec, "ARRAY_ID" ) ;
    173   *intRF = 0 ;
    174 
    175   // PROCESSOR_ID is tentatively set to 0
    176   intRF.attachToRecord( trec, "PROCESSOR_ID" ) ;
    177   *intRF = 0 ;
    178 
    179   // UVW is always [0,0,0]
    180   RecordFieldPtr< Array<Double> > uvwRF( trec, "UVW" ) ;
    181   *uvwRF = Vector<Double>( 3, 0.0 ) ;
    182 
    183   //
    184   // ITERATION: FIELDNAME
    185   //
    186   TableIterator iter0( table_->table(), "FIELDNAME" ) ;
    187   while( !iter0.pastEnd() ) {
    188     //Table t0( iter0.table() ) ;
    189     Table t0 = iter0.table() ;
    190     ROTableColumn sharedCol( t0, "FIELDNAME" ) ;
    191     String fieldName = sharedCol.asString( 0 ) ;
    192     sharedCol.attach( t0, "SRCNAME" ) ;
    193     String srcName = sharedCol.asString( 0 ) ;
    194     sharedCol.attach( t0, "TIME" ) ;
    195     Double minTime = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day->sec
    196     ROArrayColumn<Double> scanRateCol( t0, "SCANRATE" ) ;
    197     Vector<Double> scanRate = scanRateCol( 0 ) ;
    198     String::size_type pos = fieldName.find( "__" ) ;
    199     Int fieldId = -1 ;
    200     if ( pos != String::npos ) {
    201 //       os_ << "fieldName.substr( pos+2 )=" << fieldName.substr( pos+2 ) << LogIO::POST ;
    202       fieldId = String::toInt( fieldName.substr( pos+2 ) ) ;
    203       fieldName = fieldName.substr( 0, pos ) ;
    204     }
    205     else {
    206 //       os_ << "use default field id" << LogIO::POST ;
    207       fieldId = defaultFieldId ;
    208       defaultFieldId++ ;
    209     }
    210 //     os_ << "fieldId" << fieldId << ": " << fieldName << LogIO::POST ;
    211 
    212     // FIELD_ID
    213     intRF.attachToRecord( trec, "FIELD_ID" ) ;
    214     *intRF = fieldId ;
    215 
    216     //
    217     // ITERATION: BEAMNO
    218     //
    219     TableIterator iter1( t0, "BEAMNO" ) ;
    220     while( !iter1.pastEnd() ) {
    221       Table t1 = iter1.table() ;
    222       sharedCol.attach( t1, "BEAMNO" ) ;
    223       uInt beamNo = sharedCol.asuInt( 0 ) ;
    224 //       os_ << "beamNo = " << beamNo << LogIO::POST ;
    225 
    226       // FEED1 and FEED2
    227       intRF.attachToRecord( trec, "FEED1" ) ;
    228       *intRF = beamNo ;
    229       intRF.attachToRecord( trec, "FEED2" ) ;
    230       *intRF = beamNo ;
    231 
    232       //
    233       // ITERATION: SCANNO
    234       //
    235       TableIterator iter2( t1, "SCANNO" ) ;
    236       while( !iter2.pastEnd() ) {
    237         Table t2 = iter2.table() ;
    238         sharedCol.attach( t2, "SCANNO" ) ;
    239         uInt scanNo = sharedCol.asuInt( 0 ) ;
    240 //         os_ << "scanNo = " << scanNo << LogIO::POST ;
    241 
    242         // SCAN_NUMBER
    243         // MS: 1-based
    244         // Scantable: 0-based
    245         intRF.attachToRecord( trec, "SCAN_NUMBER" ) ;
    246         *intRF = scanNo + 1 ;
    247 
    248         //
    249         // ITERATION: IFNO
    250         //
    251         TableIterator iter3( t2, "IFNO" ) ;
    252         while( !iter3.pastEnd() ) {
    253           Table t3 = iter3.table() ;
    254           sharedCol.attach( t3, "IFNO" ) ;
    255           uInt ifNo = sharedCol.asuInt( 0 ) ;
    256 //           os_ << "ifNo = " << ifNo << LogIO::POST ;
    257           sharedCol.attach( t3, "FREQ_ID" ) ;
    258           uInt freqId = sharedCol.asuInt( 0 ) ;
    259 //           os_ << "freqId = " << freqId << LogIO::POST ;
    260           Int subscan = 1 ; // 1-base
    261           //
    262           // ITERATION: SRCTYPE
    263           //
    264           TableIterator iter4( t3, "SRCTYPE" ) ;
    265           while( !iter4.pastEnd() ) {
    266             Table t4 = iter4.table() ;
    267             sharedCol.attach( t4, "SRCTYPE" ) ;
    268             Int srcType = sharedCol.asInt( 0 ) ;
    269             Int stateId = addState( srcType, subscan ) ;
    270             *stateidRF = stateId ;
    271             //
    272             // ITERATION: CYCLENO and TIME
    273             //
    274             Block<String> cols( 2 ) ;
    275             cols[0] = "CYCLENO" ;
    276             cols[1] = "TIME" ;
    277             TableIterator iter5( t4, cols ) ;
    278             while( !iter5.pastEnd() ) {
    279               Table t5 =  iter5.table().sort("POLNO") ;
    280               //sharedCol.attach( t5, "CYCLENO" ) ;
    281               //uInt cycleNo = sharedCol.asuInt( 0 ) ;
    282               Int nrow = t5.nrow() ;
    283 //               os_ << "nrow = " << nrow << LogIO::POST ;
    284              
    285               Vector<Int> polnos( nrow ) ;
    286               indgen( polnos, 0 ) ;
    287               Int polid = addPolarization( polnos ) ;
    288 //               os_ << "polid = " << polid << LogIO::POST ;
    289              
    290               // DATA/FLOAT_DATA
    291               ROArrayColumn<Float> specCol( t5, "SPECTRA" ) ;
    292               ROArrayColumn<uChar> flagCol( t5, "FLAGTRA" ) ;
    293               uInt nchan = specCol( 0 ).size() ;
    294               IPosition cellshape( 2, nrow, nchan ) ;
    295               if ( useFloatData_ ) {
    296                 // FLOAT_DATA
    297                 Matrix<Float> dataArr( cellshape ) ;
    298                 Matrix<Bool> flagArr( cellshape ) ;
    299                 Vector<Bool> tmpB ;
    300                 for ( Int ipol = 0 ; ipol < nrow ; ipol++ ) {
    301                   dataArr.row( ipol ) = specCol( ipol ) ;
    302                   tmpB.reference( flagArr.row( ipol ) ) ;
    303                   convertArray( tmpB, flagCol( ipol ) ) ;
    304                 }
    305                 ((RecordFieldPtr< Array<Float> > *)dataRF)->define( dataArr ) ;
    306                
    307                 // FLAG
    308                 flagRF.define( flagArr ) ;
    309               }
    310               else if ( useData_ ) {
    311                 // DATA
    312                 // assume nrow = 4
    313                 Matrix<Complex> dataArr( cellshape ) ;
    314                 Vector<Float> zeroIm( nchan, 0 ) ;
    315                 Matrix<Float> dummy( IPosition( 2, 2, nchan ) ) ;
    316                 dummy.row( 0 ) = specCol( 0 ) ;
    317                 dummy.row( 1 ) = zeroIm ;
    318                 dataArr.row( 0 ) = RealToComplex( dummy ) ;
    319                 dummy.row( 0 ) = specCol( 1 ) ;
    320                 dataArr.row( 3 ) = RealToComplex( dummy ) ;
    321                 dummy.row( 0 ) = specCol( 2 ) ;
    322                 dummy.row( 1 ) = specCol( 3 ) ;
    323                 dataArr.row( 1 ) = RealToComplex( dummy ) ;
    324                 dataArr.row( 2 ) = conj( dataArr.row( 1 ) ) ;
    325                 ((RecordFieldPtr< Array<Complex> > *)dataRF)->define( dataArr ) ;
    326                
    327                
    328                 // FLAG
    329                 Matrix<Bool> flagArr( cellshape ) ;
    330                 Vector<Bool> tmpB ;
    331                 tmpB.reference( flagArr.row( 0 ) ) ;
    332                 convertArray( tmpB, flagCol( 0 ) ) ;
    333                 tmpB.reference( flagArr.row( 3 ) ) ;
    334                 convertArray( tmpB, flagCol( 3 ) ) ;
    335                 tmpB.reference( flagArr.row( 1 ) ) ;
    336                 convertArray( tmpB, ( flagCol( 2 ) | flagCol( 3 ) ) ) ;
    337                 flagArr.row( 2 ) = flagArr.row( 1 ) ;
    338                 flagRF.define( flagArr ) ;
    339               }
    340 
    341               // FLAG_ROW
    342               sharedCol.attach( t5, "FLAGROW" ) ;
    343               Vector<uInt> flagRowArr( nrow ) ;
    344               for ( Int irow = 0 ; irow < nrow ; irow++ )
    345                 flagRowArr[irow] = sharedCol.asuInt( irow ) ;
    346               *flagrowRF = anyNE( flagRowArr, (uInt)0 ) ;
    347 
    348               // TIME and TIME_CENTROID
    349               sharedCol.attach( t5, "TIME" ) ;
    350               Double mTimeV = (Double)sharedCol.asdouble( 0 ) * 86400.0 ; // day -> sec
    351               *timeRF = mTimeV ;
    352               *timecRF = mTimeV ;
    353 
    354               // INTERVAL and EXPOSURE
    355               sharedCol.attach( t5, "INTERVAL" ) ;
    356               Double interval = (Double)sharedCol.asdouble( 0 ) ;
    357               *intervalRF = interval ;
    358               *exposureRF = interval ;
    359              
    360               // WEIGHT and SIGMA
    361               // always 1 at the moment
    362               Vector<Float> wArr( nrow, 1.0 ) ;
    363               weightRF.define( wArr ) ;
    364               sigmaRF.define( wArr ) ;
    365              
    366               // add DATA_DESCRIPTION row
    367               Int ddid = addDataDescription( polid, ifNo ) ;
    368 //               os_ << "ddid = " << ddid << LogIO::POST ;
    369               *ddidRF = ddid ;
    370              
    371               // for SYSCAL table
    372               sharedCol.attach( t5, "TCAL_ID" ) ;
    373               Vector<uInt> tcalIdArr( nrow ) ;
    374               for ( Int irow = 0 ; irow < nrow ; irow++ )
    375                 tcalIdArr[irow] = sharedCol.asuInt( irow ) ;
    376 //               os_ << "tcalIdArr = " << tcalIdArr << LogIO::POST ;
    377               String key = String::toString( tcalIdArr[0] ) ;
    378               if ( !tcalIdRec_.isDefined( key ) ) {
    379                 tcalIdRec_.define( key, tcalIdArr ) ;
    380                 tcalRowRec_.define( key, t5.rowNumbers() ) ;
    381               }
    382               else {
    383                 Vector<uInt> pastrows = tcalRowRec_.asArrayuInt( key ) ;
    384                 tcalRowRec_.define( key, concatenateArray( pastrows, t5.rowNumbers() ) ) ;
    385               }
    386                            
    387               // for POINTING table
    388               if ( ptTabName_ == "" ) {
    389                 ROArrayColumn<Double> dirCol( t5, "DIRECTION" ) ;
    390                 Vector<Double> dir = dirCol( 0 ) ;
    391                 dirCol.attach( t5, "SCANRATE" ) ;
    392                 Vector<Double> rate = dirCol( 0 ) ;
    393                 Matrix<Double> msDir( 2, 1 ) ;
    394                 msDir.column( 0 ) = dir ;
    395                 if ( anyNE( rate, 0.0 ) ) {
    396                   msDir.resize( 2, 2, True ) ;
    397                   msDir.column( 1 ) = rate ;
    398                 }
    399                 addPointing( fieldName, mTimeV, interval, msDir ) ;
    400               }
    401              
    402               // FLAG_CATEGORY is tentatively set
    403               flagcatRF.define( Cube<Bool>( nrow, nchan, 1, False ) ) ;
    404              
    405               // add row
    406               mstable_->addRow( 1, True ) ;
    407               row.put( mstable_->nrow()-1 ) ;
    408              
    409               iter5.next() ;
    410             }
    411            
    412             iter4.next() ;
    413           }
    414 
    415           // add SPECTRAL_WINDOW row
    416           if ( allNE( processedFreqId, freqId ) ) {
    417             uInt vsize = processedFreqId.size() ;
    418             processedFreqId.resize( vsize+1, True ) ;
    419             processedFreqId[vsize] = freqId ;
    420             addSpectralWindow( ifNo, freqId ) ;
    421           }
    422          
    423           iter3.next() ;
    424         }
    425        
    426         iter2.next() ;
    427       }
    428      
    429       // add FEED row
    430       addFeed( beamNo ) ;
    431      
    432       iter1.next() ;
    433     }
    434    
    435     // add FIELD row
    436     addField( fieldId, fieldName, srcName, minTime, scanRate ) ;
    437 
    438     iter0.next() ;
    439   }
    440 
    441 //   delete tpoolr ;
    442   delete dataRF ;
    443 
    444   // SYSCAL
    445   if ( isTcal_ )
    446     fillSysCal() ;
    447 
    448   // fill empty SPECTRAL_WINDOW rows
    449   infillSpectralWindow() ;
     1412  /***
     1413   * Start iteration using TableVisitor
     1414   ***/
     1415  {
     1416    static const char *cols[] = {
     1417      "FIELDNAME", "BEAMNO", "SCANNO", "IFNO", "SRCTYPE", "CYCLENO", "TIME",
     1418      "POLNO",
     1419      NULL
     1420    };
     1421    static const TypeManagerImpl<uInt> tmUInt;
     1422    static const TypeManagerImpl<Int> tmInt;
     1423    static const TypeManagerImpl<Double> tmDouble;
     1424    static const TypeManagerImpl<String> tmString;
     1425    static const TypeManager *const tms[] = {
     1426      &tmString, &tmUInt, &tmUInt, &tmUInt, &tmInt, &tmUInt, &tmDouble, &tmInt, NULL
     1427    };
     1428    //double t0 = mathutil::gettimeofday_sec() ;
     1429    MSWriterVisitor myVisitor(table_->table(),*mstable_);
     1430    //double t1 = mathutil::gettimeofday_sec() ;
     1431    //cout << "MSWriterVisitor(): elapsed time " << t1-t0 << " sec" << endl ;
     1432    String dataColName = "FLOAT_DATA" ;
     1433    if ( useData_ )
     1434      dataColName = "DATA" ;
     1435    myVisitor.dataColumnName( dataColName ) ;
     1436    myVisitor.pointingTableName( ptTabName_ ) ;
     1437    myVisitor.setSourceRecord( srcRec_ ) ;
     1438    //double t2 = mathutil::gettimeofday_sec() ;
     1439    traverseTable(table_->table(), cols, tms, &myVisitor);
     1440    //double t3 = mathutil::gettimeofday_sec() ;
     1441    //cout << "traverseTable(): elapsed time " << t3-t2 << " sec" << endl ;
     1442    map< Int,Vector<uInt> > &idRec = myVisitor.getTcalIdRecord() ;
     1443    map< Int,Vector<uInt> > &rowRec = myVisitor.getTcalRowRecord() ;
     1444
     1445    // SYSCAL
     1446    if ( isTcal_ )
     1447      fillSysCal( idRec, rowRec ) ;
     1448  }
     1449  /***
     1450   * End iteration using TableVisitor
     1451   ***/
    4501452
    4511453  // ASDM tables
     
    4751477  }
    4761478
    477 //   double endSec = mathutil::gettimeofday_sec() ;
    478 //   os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1479  //double endSec = mathutil::gettimeofday_sec() ;
     1480  //os_ << "end MSWriter::write() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    4791481
    4801482  return True ;
    4811483}
    482  
     1484
    4831485void MSWriter::init()
    4841486{
     
    6971699void MSWriter::fillObservation()
    6981700{
    699 //   double startSec = mathutil::gettimeofday_sec() ;
    700 //   os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
     1701  //double startSec = mathutil::gettimeofday_sec() ;
     1702  //os_ << "start MSWriter::fillObservation() startSec=" << startSec << LogIO::POST ;
    7011703
    7021704  // only 1 row
     
    7261728  msObsCols.timeRangeMeas().put( 0, trange ) ;
    7271729
    728 //   double endSec = mathutil::gettimeofday_sec() ;
    729 //   os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1730  //double endSec = mathutil::gettimeofday_sec() ;
     1731  //os_ << "end MSWriter::fillObservation() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    7301732}
     1733
     1734void MSWriter::antennaProperty( String &name, String &m, String &t, Double &d )
     1735{
     1736  name.upcase() ;
     1737 
     1738  m = "ALT-AZ" ;
     1739  t = "GROUND-BASED" ;
     1740  if ( name.matches( Regex( "DV[0-9]+$" ) )
     1741       || name.matches( Regex( "DA[0-9]+$" ) )
     1742       || name.matches( Regex( "PM[0-9]+$" ) ) )
     1743    d = 12.0 ;
     1744  else if ( name.matches( Regex( "CM[0-9]+$" ) ) )
     1745    d = 7.0 ;
     1746  else if ( name.contains( "GBT" ) )
     1747    d = 104.9 ;
     1748  else if ( name.contains( "MOPRA" ) )
     1749    d = 22.0 ;
     1750  else if ( name.contains( "PKS" ) || name.contains( "PARKS" ) )
     1751    d = 64.0 ;
     1752  else if ( name.contains( "TIDBINBILLA" ) )
     1753    d = 70.0 ;
     1754  else if ( name.contains( "CEDUNA" ) )
     1755    d = 30.0 ;
     1756  else if ( name.contains( "HOBART" ) )
     1757    d = 26.0 ;
     1758  else if ( name.contains( "APEX" ) )
     1759    d = 12.0 ;
     1760  else if ( name.contains( "ASTE" ) )
     1761    d = 10.0 ;
     1762  else if ( name.contains( "NRO" ) )
     1763    d = 45.0 ;
     1764  else
     1765    d = 1.0 ;
     1766}
    7311767
    7321768void MSWriter::fillAntenna()
    7331769{
    734 //   double startSec = mathutil::gettimeofday_sec() ;
    735 //   os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
     1770  //double startSec = mathutil::gettimeofday_sec() ;
     1771  //os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
    7361772
    7371773  // only 1 row
    738   mstable_->antenna().addRow( 1, True ) ;
    739   MSAntennaColumns msAntCols( mstable_->antenna() ) ;
    740 
    741   String hAntennaName = header_.antennaname ;
    742   String::size_type pos = hAntennaName.find( "//" ) ;
     1774  Table anttab = mstable_->antenna() ;
     1775  anttab.addRow( 1, True ) ;
     1776 
     1777  Table &table = table_->table() ;
     1778  const TableRecord &keys = table.keywordSet() ;
     1779  String hAntName = keys.asString( "AntennaName" ) ;
     1780  String::size_type pos = hAntName.find( "//" ) ;
    7431781  String antennaName ;
    7441782  String stationName ;
    7451783  if ( pos != String::npos ) {
    746     hAntennaName = hAntennaName.substr( pos+2 ) ;
    747   }
    748   pos = hAntennaName.find( "@" ) ;
     1784    stationName = hAntName.substr( 0, pos ) ;
     1785    hAntName = hAntName.substr( pos+2 ) ;
     1786  }
     1787  pos = hAntName.find( "@" ) ;
    7491788  if ( pos != String::npos ) {
    750     antennaName = hAntennaName.substr( 0, pos ) ;
    751     stationName = hAntennaName.substr( pos+1 ) ;
     1789    antennaName = hAntName.substr( 0, pos ) ;
     1790    stationName = hAntName.substr( pos+1 ) ;
    7521791  }
    7531792  else {
    754     antennaName = hAntennaName ;
    755     stationName = hAntennaName ;
    756   }
    757 //   os_ << "antennaName = " << antennaName << LogIO::POST ;
    758 //   os_ << "stationName = " << stationName << LogIO::POST ;
     1793    antennaName = hAntName ;
     1794  }
     1795  Vector<Double> antpos = keys.asArrayDouble( "AntennaPosition" ) ;
    7591796 
    760   msAntCols.name().put( 0, antennaName ) ;
    761   msAntCols.station().put( 0, stationName ) ;
    762 
    763 //   os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
     1797  String mount, atype ;
     1798  Double diameter ;
     1799  antennaProperty( antennaName, mount, atype, diameter ) ;
    7641800 
    765   msAntCols.position().put( 0, header_.antennaposition ) ;
    766 
    767   // MOUNT is set to "ALT-AZ"
    768   msAntCols.mount().put( 0, "ALT-AZ" ) ;
    769 
    770   Double diameter = getDishDiameter( antennaName ) ;
    771   msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ;
    772 
    773 //   double endSec = mathutil::gettimeofday_sec() ;
    774 //   os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1801  TableRow tr( anttab ) ;
     1802  TableRecord &r = tr.record() ;
     1803  RecordFieldPtr<String> nameRF( r, "NAME" ) ;
     1804  RecordFieldPtr<String> stationRF( r, "STATION" ) ;
     1805  RecordFieldPtr<String> mountRF( r, "NAME" ) ;
     1806  RecordFieldPtr<String> typeRF( r, "TYPE" ) ;
     1807  RecordFieldPtr<Double> dishDiameterRF( r, "DISH_DIAMETER" ) ;
     1808  RecordFieldPtr< Vector<Double> > positionRF( r, "POSITION" ) ;
     1809  *nameRF = antennaName ;
     1810  *mountRF = mount ;
     1811  *typeRF = atype ;
     1812  *dishDiameterRF = diameter ;
     1813  *positionRF = antpos ;
     1814 
     1815  tr.put( 0 ) ;
     1816
     1817  //double endSec = mathutil::gettimeofday_sec() ;
     1818  //os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    7751819}
     1820 
     1821// void MSWriter::fillAntenna()
     1822// {
     1823// //   double startSec = mathutil::gettimeofday_sec() ;
     1824// //   os_ << "start MSWriter::fillAntenna() startSec=" << startSec << LogIO::POST ;
     1825
     1826//   // only 1 row
     1827//   mstable_->antenna().addRow( 1, True ) ;
     1828//   MSAntennaColumns msAntCols( mstable_->antenna() ) ;
     1829
     1830//   String hAntennaName = header_.antennaname ;
     1831//   String::size_type pos = hAntennaName.find( "//" ) ;
     1832//   String antennaName ;
     1833//   String stationName ;
     1834//   if ( pos != String::npos ) {
     1835//     hAntennaName = hAntennaName.substr( pos+2 ) ;
     1836//   }
     1837//   pos = hAntennaName.find( "@" ) ;
     1838//   if ( pos != String::npos ) {
     1839//     antennaName = hAntennaName.substr( 0, pos ) ;
     1840//     stationName = hAntennaName.substr( pos+1 ) ;
     1841//   }
     1842//   else {
     1843//     antennaName = hAntennaName ;
     1844//     stationName = hAntennaName ;
     1845//   }
     1846// //   os_ << "antennaName = " << antennaName << LogIO::POST ;
     1847// //   os_ << "stationName = " << stationName << LogIO::POST ;
     1848 
     1849//   msAntCols.name().put( 0, antennaName ) ;
     1850//   msAntCols.station().put( 0, stationName ) ;
     1851
     1852// //   os_ << "antennaPosition = " << header_.antennaposition << LogIO::POST ;
     1853 
     1854//   msAntCols.position().put( 0, header_.antennaposition ) ;
     1855
     1856//   // MOUNT is set to "ALT-AZ"
     1857//   msAntCols.mount().put( 0, "ALT-AZ" ) ;
     1858
     1859//   Double diameter = getDishDiameter( antennaName ) ;
     1860//   msAntCols.dishDiameterQuant().put( 0, Quantity( diameter, "m" ) ) ;
     1861
     1862// //   double endSec = mathutil::gettimeofday_sec() ;
     1863// //   os_ << "end MSWriter::fillAntenna() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     1864// }
    7761865
    7771866void MSWriter::fillProcessor()
     
    8361925    ROScalarColumn<Double> srcVelCol( t0, "SRCVELOCITY" ) ;
    8371926    Double srcVel = srcVelCol( 0 ) ;
     1927    srcRec_.define( srcName, srcId ) ;
    8381928
    8391929    // NAME
     
    9902080}
    9912081
    992 void MSWriter::fillSysCal()
     2082void MSWriter::fillSysCal( map< Int,Vector<uInt> > &idRec,
     2083                           map< Int,Vector<uInt> > &rowRec )
    9932084{
    994 //   double startSec = mathutil::gettimeofday_sec() ;
    995 //   os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
    996 
    997   //tcalIdRec_.print( cout ) ;
     2085  //double startSec = mathutil::gettimeofday_sec() ;
     2086  //os_ << "start MSWriter::fillSysCal() startSec=" << startSec << LogIO::POST ;
     2087
     2088  //idRec.print( cout ) ;
    9982089
    9992090  // access to MS SYSCAL subtable
     
    10172108    return ;
    10182109
    1019   nrow = tcalIdRec_.nfields() ;
     2110  //nrow = idRec.nfields() ;
     2111  nrow = idRec.size() ;
    10202112
    10212113  Double midTime ;
     
    10522144  ROTableColumn beamnoCol( tab, "BEAMNO" ) ;
    10532145  ROTableColumn ifnoCol( tab, "IFNO" ) ;
     2146  map< Int,Vector<uInt> >::iterator itr0 = idRec.begin() ;
     2147  map< Int,Vector<uInt> >::iterator itr1 = rowRec.begin() ;
    10542148  for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
    10552149//     double t1 = mathutil::gettimeofday_sec() ;
    1056     Vector<uInt> ids = tcalIdRec_.asArrayuInt( irow ) ;
     2150    Vector<uInt> ids = itr0->second ;
     2151    itr0++ ;
    10572152//     os_ << "ids = " << ids << LogIO::POST ;
    10582153    uInt npol = ids.size() ;
    1059     Vector<uInt> rows = tcalRowRec_.asArrayuInt( irow ) ;
     2154    Vector<uInt> rows = itr1->second ;
     2155    itr1++ ;
    10602156//     os_ << "rows = " << rows << LogIO::POST ;
    10612157    Vector<Double> atime( rows.nelements() ) ;
     
    11532249    if ( tcalSpec_ ) {
    11542250      // put TCAL_SPECTRUM
    1155       //*tcalspRF = tcal ;
    11562251      tcalspRF.define( tcal ) ;
    11572252      // set TCAL (mean of TCAL_SPECTRUM)
     
    11612256      }
    11622257      // put TCAL
    1163       *tcalRF = tcalMean ;
     2258      tcalRF.define( tcalMean ) ;
    11642259    }
    11652260    else {
    11662261      // put TCAL
    1167       *tcalRF = tcal ;
     2262      tcalRF.define( tcal ) ;
    11682263    }
    11692264   
    11702265    if ( tsysSpec_ ) {
    11712266      // put TSYS_SPECTRUM
    1172       //*tsysspRF = tsys ;
    11732267      tsysspRF.define( tsys ) ;
    11742268      // set TSYS (mean of TSYS_SPECTRUM)
     
    11782272      }
    11792273      // put TSYS
    1180       *tsysRF = tsysMean ;
     2274      tsysRF.define( tsysMean ) ;
    11812275    }
    11822276    else {
    11832277      // put TSYS
    1184       *tsysRF = tsys ;
     2278      tsysRF.define( tsys ) ;
    11852279    }
    11862280
     
    11932287  }
    11942288 
    1195 //   double endSec = mathutil::gettimeofday_sec() ;
    1196 //   os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
     2289  //double endSec = mathutil::gettimeofday_sec() ;
     2290  //os_ << "end MSWriter::fillSysCal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
    11972291}
    11982292
  • trunk/src/MSWriter.h

    r2258 r2291  
    2323#include <casa/Logging/LogIO.h>
    2424#include <casa/Containers/Record.h>
     25#include <casa/Containers/RecordField.h>
    2526
    2627#include <tables/Tables/Table.h>
     
    3940namespace asap
    4041{
     42class MSWriterUtils
     43{
     44protected:
     45  template<class T> void putField( const String &name,
     46                                   TableRecord &r,
     47                                   T &val )
     48  {
     49    RecordFieldPtr<T> rf( r, name ) ;
     50    *rf = val ;
     51  }
     52  template<class T> void defineField( const String &name,
     53                                      TableRecord &r,
     54                                      T &val )
     55  {
     56    RecordFieldPtr<T> rf( r, name ) ;
     57    rf.define( val ) ;
     58  }
     59};
    4160
    4261class MSWriter
     
    6584  void fillSource() ;
    6685  void fillWeather() ;
    67   void fillSysCal() ;
     86  void fillSysCal( std::map< casa::Int,casa::Vector<casa::uInt> > &idrec,
     87                   std::map< casa::Int,casa::Vector<casa::uInt> > &rowrec ) ;
     88//   void fillSysCal( Record &idrec, Record &rowrec ) ;
     89//   void fillSysCal() ;
    6890
    6991  // fill empty rows
     
    86108  void queryType( casa::Int type, casa::String &stype, casa::Bool &b, casa::Double &t, Double &l ) ;
    87109  casa::Double getDishDiameter( casa::String antname ) ;
     110  void antennaProperty( casa::String &name, casa::String &mount, casa::String &type, casa::Double &diameter ) ;
    88111
    89112  // tool for HPC
     
    111134  casa::LogIO os_ ;
    112135
    113   casa::Record tcalIdRec_ ;
    114   casa::Record tcalRowRec_ ;
     136//   casa::Record tcalIdRec_ ;
     137//   casa::Record tcalRowRec_ ;
     138  casa::Record srcRec_ ;
    115139 
    116140  MSWriter();
  • trunk/src/TableTraverse.cpp

    r2289 r2291  
    251251      sizes[i] = typeManagers[i]->sizeOf();
    252252      BaseColumn *baseCol = ROTableColumnBackDoor::baseColumnPtr(cols[i]);
    253       PlainColumn *col = dynamic_cast <PlainColumn *>(baseCol);
    254       if (col) {
    255         const uInt gotRows = col->dataManagerColumn()->
    256           getBlockV(0, nRows, colValues[i]);
    257         DebugAssert(gotRows == nRows, AipsError);
    258       } else {
    259         copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol);
    260       }
     253      copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol);
     254//       PlainColumn *col = dynamic_cast <PlainColumn *>(baseCol);
     255//       if (col) {
     256//      const uInt gotRows = col->dataManagerColumn()->
     257//        getBlockV(0, nRows, colValues[i]);
     258//      DebugAssert(gotRows == nRows, AipsError);
     259//       } else {
     260//      copyColumnData(colValues[i], typeManagers[i]->sizeOf(), nRows, baseCol);
     261//       }
    261262    }
    262263
Note: See TracChangeset for help on using the changeset viewer.