source: trunk/src/MSFiller.cpp @ 1975

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

New Development: Yes

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Writer improvement.
File size: 53.6 KB
Line 
1//
2// C++ Interface: MSFiller
3//
4// Description:
5//
6// This class is specific filler for MS format
7//
8// Takeshi Nakazato <takeshi.nakazato@nao.ac.jp>, (C) 2010
9//
10// Copyright: See COPYING file that comes with this distribution
11//
12//
13
14#include <iostream>
15#include <map>
16
17#include <tables/Tables/ExprNode.h>
18#include <tables/Tables/TableIter.h>
19#include <tables/Tables/ScalarColumn.h>
20#include <tables/Tables/ArrayColumn.h>
21#include <tables/Tables/RefRows.h>
22
23#include <casa/Containers/Block.h>
24#include <casa/Logging/LogIO.h>
25#include <casa/Arrays/Slicer.h>
26#include <casa/Quanta/MVTime.h>
27
28#include <measures/Measures/Stokes.h>
29#include <measures/Measures/MEpoch.h>
30#include <measures/TableMeasures/ScalarMeasColumn.h>
31
32#include <atnf/PKSIO/SrcType.h>
33
34#include "MSFiller.h"
35#include "STHeader.h"
36
37using namespace casa ;
38using namespace std ;
39
40namespace asap {
41MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
42  : table_( stable ),
43    antenna_( -1 ),
44    getPt_( False ),
45    isFloatData_( False ),
46    isData_( False ),
47    isDoppler_( False ),
48    isFlagCmd_( False ),
49    isFreqOffset_( False ),
50    isHistory_( False ),
51    isProcessor_( False ),
52    isSysCal_( False ),
53    isWeather_( False )
54{
55  os_ = LogIO() ;
56  os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
57}
58
59MSFiller::~MSFiller()
60{
61  os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
62}
63
64bool MSFiller::open( const std::string &filename, const casa::Record &rec )
65{
66  os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
67  //os_ << "   filename = " << filename << endl ;
68  //rec.print( os_.output(), 25, "      " ) ;
69  //os_ << LogIO::POST ;
70
71  // parsing MS options
72  //Int antenna = 0 ;
73  //Bool getPt = False;
74
75  if ( rec.isDefined( "ms" ) ) {
76    Record msrec = rec.asRecord( "ms" ) ;
77    if ( msrec.isDefined( "getpt" ) ) {
78      getPt_ = msrec.asBool( "getpt" ) ;
79    }
80    if ( msrec.isDefined( "antenna" ) ) {
81      if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
82        antenna_ = msrec.asInt( "antenna" ) ;
83      }
84      else {
85        antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
86      }
87    }
88    else {
89      antenna_ = 0 ;
90    }
91  }
92
93  os_ << "Parsing MS options" << endl ;
94  os_ << "   getPt = " << getPt_ << endl ;
95  os_ << "   antenna = " << antenna_ << LogIO::POST ;
96
97  mstable_ = MeasurementSet( filename, Table::Old ) ;
98
99  // check which data column exists
100  isFloatData_ = mstable_.isColumn( MSMainEnums::FLOAT_DATA ) ;
101  isData_ = mstable_.isColumn( MSMainEnums::DATA ) ;
102
103  return true ;
104}
105
106void MSFiller::fill()
107{
108  os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
109
110  tablesel_ = mstable_( mstable_.col("ANTENNA1") == mstable_.col("ANTENNA2") 
111                        && mstable_.col("ANTENNA1") == antenna_ ) ;
112
113  // Initialize header
114  STHeader sdh ; 
115  sdh.nchan = 0 ;
116  sdh.npol = 0 ;
117  sdh.nif = 0 ;
118  sdh.nbeam = 0 ;
119  sdh.observer = "" ;
120  sdh.project = "" ;
121  sdh.obstype = "" ;
122  sdh.antennaname = "" ;
123  sdh.antennaposition.resize( 0 ) ;
124  sdh.equinox = 0.0 ;
125  sdh.freqref = "" ;
126  sdh.reffreq = -1.0 ;
127  sdh.bandwidth = 0.0 ;
128  sdh.utc = 0.0 ;
129  sdh.fluxunit = "" ;
130  sdh.epoch = "" ;
131  sdh.poltype = "" ;
132 
133  // check if optional table exists
134  const TableRecord msrec = tablesel_.keywordSet() ;
135  isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
136  isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
137  isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
138  isHistory_ = msrec.isDefined( "HISTORY" ) ;
139  isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
140  isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
141  isWeather_ = msrec.isDefined( "WEATHER" ) ;
142 
143  // Access to MS subtables
144  MSField fieldtab = tablesel_.field() ;
145  MSPolarization poltab = tablesel_.polarization() ;
146  MSDataDescription ddtab = tablesel_.dataDescription() ;
147  MSObservation obstab = tablesel_.observation() ;
148  MSSource srctab = tablesel_.source() ;
149  MSSpectralWindow spwtab = tablesel_.spectralWindow() ;
150  MSSysCal caltab = tablesel_.sysCal() ;
151  if ( caltab.nrow() == 0 )
152    isSysCal_ = False ;
153  MSPointing pointtab = tablesel_.pointing() ;
154  MSWeather weathertab = tablesel_.weather() ;
155  if ( weathertab.nrow() == 0 )
156    isWeather_ = False ;
157
158//   os_ << "isDoppler_ = " << isDoppler_ << endl
159//       << "isFlagCmd_ = " << isFlagCmd_ << endl
160//       << "isFreqOffset_ = " << isFreqOffset_ << endl
161//       << "isHistory_ = " << isHistory_ << endl
162//       << "isProcessor_ = " << isProcessor_ << endl
163//       << "isSysCal_ = " << isSysCal_ << endl
164//       << "isWeather_ = " << isWeather_ << LogIO::POST ;
165
166 // SUBTABLES: FREQUENCIES
167  //fillFrequencies() ;
168  table_->frequencies().setFrame( "LSRK" ) ;
169  table_->frequencies().setFrame( "LSRK", True ) ;
170
171  // SUBTABLES: WEATHER
172  if ( isWeather_ )
173    fillWeather() ;
174
175  // SUBTABLES: FOCUS
176  fillFocus() ;
177
178  // SUBTABLES: TCAL
179  if ( isSysCal_ )
180    fillTcal() ;
181
182  // SUBTABLES: MOLECULES
183  //fillMolecules() ;
184
185  // SUBTABLES: FIT
186  //fillFit() ;
187
188  // SUBTABLES: HISTORY
189  //fillHistory() ;
190
191  // MAIN
192  // Iterate over several ids
193  //
194  // ITERATION: OBSERVATION_ID
195  //
196  TableIterator iter0( tablesel_, "OBSERVATION_ID" ) ;
197  Int totalrow = 0 ;
198  Int oldnr = table_->nrow() ;
199  map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair
200  ROMSAntennaColumns antCols( mstable_.antenna() ) ;
201  Vector< Quantum<Double> > antpos = antCols.positionQuant()(antenna_) ;
202  MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
203  MSPointing potabsel = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;
204  String stationName = antCols.station()(antenna_) ;
205  ROMSPointingColumns pointCols( potabsel ) ;
206  String telescopeName ;
207  while( !iter0.pastEnd() ) {
208    MeasurementSet t0( iter0.table() ) ;
209    ROScalarColumn<Int> mObsIdCol( t0, "OBSERVATION_ID" ) ;
210    Int obsId = mObsIdCol( 0 ) ;
211    ROMSObservationColumns obsCols( obstab ) ;
212    if ( sdh.observer == "" ) sdh.observer = obsCols.observer()(obsId) ;
213    if ( sdh.project == "" ) sdh.project = obsCols.project()(obsId) ;
214    MEpoch me = obsCols.timeRangeMeas()(obsId)(IPosition(1,0)) ;
215    if ( sdh.utc == 0.0 ) {
216      Quantum<Double> startTime = obsCols.timeRangeQuant()(obsId)(IPosition(1,0)) ;
217      sdh.utc = startTime.getValue( "s" ) ;
218    }
219    telescopeName = obsCols.telescopeName()(obsId) ;
220    //
221    // ITERATION: FEED1
222    //
223    Int nbeam = 0 ;
224    TableIterator iter1( t0, "FEED1" ) ;
225    while( !iter1.pastEnd() ) {
226      MeasurementSet t1( iter1.table() ) ;
227      ROScalarColumn<Int> mFeedCol( t1, "FEED1" ) ;
228      Int feedId = mFeedCol( 0 ) ; // assume FEED1 == FEED2
229      nbeam++ ;
230      //
231      // ITERATION: FIELD_ID
232      //
233      TableIterator iter2( t1, "FIELD_ID" ) ;
234      while( !iter2.pastEnd() ) {
235        MeasurementSet t2( iter2.table() ) ;
236        ROScalarColumn<Int> mFieldIdCol( t2, "FIELD_ID" ) ;
237        Int fieldId = mFieldIdCol( 0 ) ;
238        ROScalarColumn<String> mFieldNameCol( fieldtab, "NAME" ) ;
239        ROScalarColumn<Int> mSrcIdCol( fieldtab, "SOURCE_ID" ) ;
240        String fieldName = mFieldNameCol( fieldId ) + "__" + String::toString(fieldId) ;
241        Int srcId = mSrcIdCol( fieldId ) ;
242        //
243        // ITERATION: DATA_DESC_ID
244        //
245        TableIterator iter3( t2, "DATA_DESC_ID" ) ;
246        while( !iter3.pastEnd() ) {
247          MeasurementSet t3( iter3.table() ) ;
248          ROScalarColumn<Int> mDDIdCol( t3, "DATA_DESC_ID" ) ;
249          Int ddId = mDDIdCol( 0 ) ;
250          ROMSDataDescColumns ddCols( ddtab ) ;
251          Int polId = ddCols.polarizationId()(ddId) ;
252          Int spwId = ddCols.spectralWindowId()(ddId) ;
253          // polarization information
254          ROMSPolarizationColumns polCols( poltab ) ;
255          Int npol = polCols.numCorr()(polId) ;
256          Vector<Int> corrtype = polCols.corrType()(polId) ;
257          //os_ << "npol = " << npol << LogIO::POST ;
258          //os_ << "corrtype = " << corrtype << LogIO::POST ;
259          if ( sdh.npol < npol ) sdh.npol = npol ;
260          if ( sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
261          // source information
262          MSSource srctabSel( srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ) ;
263          //os_ << "srcId = " << srcId << " spwId = " << spwId << " nrow = " << srctabSel.nrow() << LogIO::POST ;
264          ROMSSourceColumns srcCols( srctabSel ) ;
265          String srcName = srcCols.name()(0) ;
266          //os_ << "srcName = " << srcName << LogIO::POST ;
267          Array<Double> srcPM = srcCols.properMotion()(0) ;
268          //os_ << "srcPM = " << srcPM << LogIO::POST ;
269          Array<Double> srcDir = srcCols.direction()(0) ;
270          //os_ << "srcDir = " << srcDir << LogIO::POST ;
271          Array<Double> sysVels = srcCols.sysvel()(0) ;
272          Double sysVel = 0.0 ;
273          if ( !sysVels.empty() ) {
274            //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ;
275            // NB: assume all SYSVEL values are the same
276            Double sysVel = sysVels( IPosition(1,0) ) ;
277          }
278          //os_ << "sysVel = " << sysVel << LogIO::POST ;
279          MDirection md = srcCols.directionMeas()(0) ;
280          // for MOLECULES subtable
281          Int numLines = srcCols.numLines()(0) ;
282          //os_ << "numLines = " << numLines << LogIO::POST ;
283          Vector<Double> restFreqs( numLines, 0.0 ) ;
284          Vector<String> transitionName( numLines, "" ) ;
285          if ( numLines != 0 ) {
286            Array< Quantum<Double> > qRestFreqs = srcCols.restFrequencyQuant()(0) ;
287            for ( int i = 0 ; i < numLines ; i++ ) {
288              restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
289            }
290            //os_ << "restFreqs = " << restFreqs << LogIO::POST ;
291            if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
292              ROScalarColumn<String> transitionNameCol = srcCols.transition() ;
293              //os_ << "transitionNameCol.nrow() = " << transitionNameCol.nrow() << LogIO::POST ;
294              transitionName = transitionNameCol(0) ;
295            }
296          }
297          uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ;
298          // spectral setup
299          MeasFrame mf( me, mp, md ) ;
300          ROMSSpWindowColumns spwCols( spwtab ) ;
301          MFrequency::Types freqRef = MFrequency::castType((uInt)(spwCols.measFreqRef()(spwId))) ;
302          Int nchan = spwCols.numChan()(spwId) ;
303          Bool even = False ;
304          if ( (nchan/2)*2 == nchan ) even = True ;
305          if ( sdh.nchan < nchan ) sdh.nchan = nchan ;
306          Quantum<Double> qtotbw = spwCols.totalBandwidthQuant()(spwId) ;
307          Double totbw = qtotbw.getValue( "Hz" ) ;
308          if ( sdh.bandwidth < totbw ) sdh.bandwidth = totbw ;
309          if ( sdh.freqref == "" )
310            //sdh.freqref = MFrequency::showType( freqRef ) ;
311            sdh.freqref = "LSRK" ;
312          if ( sdh.reffreq == -1.0 ) {
313            Quantum<Double> qreffreq = spwCols.refFrequencyQuant()(spwId) ;
314            if ( freqRef == MFrequency::LSRK ) {
315              sdh.reffreq = qreffreq.getValue("Hz") ;
316            }
317            else {
318              MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
319              sdh.reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
320            }
321          }
322          Int refchan = nchan / 2 ;
323          IPosition refip( 1, refchan ) ;
324          Double refpix = 0.5*(nchan-1) ;
325          Double refval = 0.0 ;
326          Double increment = spwCols.chanWidthQuant()(spwId)(refip).getValue("Hz") ;
327          //os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
328          Vector< Quantum<Double> > chanFreqs = spwCols.chanFreqQuant()(spwId) ;
329          if ( freqRef == MFrequency::LSRK ) {
330            if ( even ) {
331              IPosition refip0( 1, refchan-1 ) ;
332              Double refval0 = chanFreqs(refip0).getValue("Hz") ;
333              Double refval1 = chanFreqs(refip).getValue("Hz") ;
334              refval = 0.5 * ( refval0 + refval1 ) ;
335            }
336            else {
337              refval = chanFreqs(refip).getValue("Hz") ;
338            }
339          }
340          else {
341            MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
342            if ( even ) {
343              IPosition refip0( 1, refchan-1 ) ;
344              Double refval0 = chanFreqs(refip0).getValue("Hz") ;
345              Double refval1 = chanFreqs(refip).getValue("Hz") ;
346              refval = 0.5 * ( refval0 + refval1 ) ;
347              refval = tolsr( refval ).get("Hz").getValue() ;
348            }
349            else {
350              refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
351            }
352          }
353          uInt freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
354          if ( ifmap.find( spwId ) == ifmap.end() ) {
355            ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
356            //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
357          }
358          // for TSYS and TCAL
359          MSSysCal caltabsel( caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ) ;
360          //Vector<uInt> tcalidrange = addTcal( caltabsel ) ;
361          //
362          // ITERATION: SCAN_NUMBER
363          //
364          TableIterator iter4( t3, "SCAN_NUMBER" ) ;
365          while( !iter4.pastEnd() ) {
366            MeasurementSet t4( iter4.table() ) ;
367            ROScalarColumn<Int> mScanNumCol( t4, "SCAN_NUMBER" ) ;
368            Int scanNum = mScanNumCol( 0 ) ;
369            uInt cycle = 0 ;
370            //
371            // ITERATION: STATE_ID
372            //
373            TableIterator iter5( t4, "STATE_ID" ) ;
374            while( !iter5.pastEnd() ) {
375              MeasurementSet t5( iter5.table().sort( "TIME" ) ) ;
376              ROScalarColumn<Int> mStateIdCol( t5, "STATE_ID" ) ;
377              Int stateId = mStateIdCol( 0 ) ;
378              ROMSStateColumns stateCols( t5.state() ) ;
379              String obstype = stateCols.obsMode()(stateId) ;
380              if ( sdh.obstype == "" ) sdh.obstype = obstype ;
381
382              Int nrow = t5.nrow() ;
383              Int prevnr = oldnr ;
384              Int addednr = 0 ;
385           
386              // SPECTRA, FLAG
387              ArrayColumn<Float> spCol( table_->table(), "SPECTRA" ) ;
388              ArrayColumn<uChar> flagCol( table_->table(), "FLAGTRA" ) ;
389              ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
390              if ( isFloatData_ ) {
391                //os_ << "FLOAT_DATA exists" << LogIO::POST ;
392                ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ;
393                IPosition cShape = mFloatDataCol.shape( 0 ) ;
394                IPosition newShape( 2, cShape[1], nrow ) ;
395                for ( int ipol = 0 ; ipol < npol ; ipol++ ) {
396                  table_->table().addRow( nrow ) ;
397                  addednr += nrow ;
398                  Int newnr = oldnr + nrow ;
399                  RefRows rows( oldnr, newnr-1 ) ;
400                  Slice paxis( ipol, 1, 1 ) ;
401                  Slice caxis( 0, cShape[1], 1 ) ;
402                  Slicer slicer( paxis, caxis ) ;
403                  spCol.putColumnCells( rows, mFloatDataCol.getColumn(slicer).reform(newShape) ) ;
404                  Array<Bool> flags = mFlagCol.getColumn(slicer).reform(newShape) ;
405                  Array<uChar> flagtra( flags.shape() ) ;
406                  convertArray( flagtra, flags ) ;
407                  flagCol.putColumnCells( rows, flagtra ) ;
408                  oldnr = newnr ;
409                }
410                if ( sdh.fluxunit == "" ) {
411                  const TableRecord dataColKeys = mFloatDataCol.keywordSet() ;
412                  if ( dataColKeys.isDefined( "UNIT" ) )
413                    sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
414                }
415              }
416              else if ( isData_ ) {
417                //os_ << "DATA exists" << LogIO::POST ;
418                ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ;
419                IPosition cShape = mDataCol.shape( 0 ) ;
420                IPosition newShape( 2, cShape[1], nrow ) ;
421                Bool crossOK = False ;
422                for ( int ipol = 0 ; ipol < npol ; ipol++ ) {
423                  //os_ << "corrtype[" << ipol << "] = " << corrtype[ipol] << LogIO::POST ;
424                  if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
425                       || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
426                    if ( !crossOK ) {
427                      //os_ << "cross polarization data" << LogIO::POST ;
428                      table_->table().addRow( nrow, True ) ;
429                      addednr += nrow ;
430                      //os_ << "table_->nrow() = " << table_->nrow() << LogIO::POST ;
431                      Int newnr = oldnr + nrow ;
432                      RefRows rows( oldnr, newnr-1 ) ;
433                      Slice paxis( ipol, 1, 1 ) ;
434                      Slice caxis( 0, cShape[1], 1 ) ;
435                      Slicer slicer( paxis, caxis ) ;
436                      Array<Complex> data = mDataCol.getColumn(slicer).reform(newShape) ;
437                      spCol.putColumnCells( rows, real( data ) ) ;
438                      Array<Bool> flags = mFlagCol.getColumn(slicer).reform(newShape) ;
439                      Array<uChar> flagtra( flags.shape() ) ;
440                      convertArray( flagtra, flags ) ;
441                      flagCol.putColumnCells( rows, flagtra ) ;
442                      oldnr = newnr ;
443                      table_->table().addRow( nrow, True ) ;
444                      addednr += nrow ;
445                      newnr = oldnr + nrow ;
446                      rows = RefRows( oldnr, newnr-1 ) ;
447                      if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) {
448                        data = conj( data ) ;
449                      }
450                      spCol.putColumnCells( rows, imag( data ) ) ;
451                      flagCol.putColumnCells( rows, flagtra ) ;
452                      crossOK = True ;
453                      oldnr = newnr ;
454                    }
455                  }
456                  else {
457                    table_->table().addRow( nrow, True ) ;
458                    addednr += nrow ;
459                    //os_ << "table_->nrow() = " << table_->nrow() << LogIO::POST ;
460                    Int newnr = oldnr + nrow ;
461                    RefRows rows( oldnr, newnr-1 ) ;
462                    Slice paxis( ipol, 1, 1 ) ;
463                    Slice caxis( 0, cShape[1], 1 ) ;
464                    Slicer slicer( paxis, caxis ) ;
465                    Array<Complex> data = mDataCol.getColumn(slicer).reform(newShape) ;
466                    spCol.putColumnCells( rows, real( data ) ) ;
467                    Array<Bool> flags = mFlagCol.getColumn(slicer).reform(newShape) ;
468                    Array<uChar> flagtra( flags.shape() ) ;
469                    convertArray( flagtra, flags ) ;
470                    flagCol.putColumnCells( rows, flagtra ) ;
471                    oldnr = newnr ;
472                  }
473                }
474                if ( sdh.fluxunit == "" ) {
475                  const TableRecord dataColKeys = mDataCol.keywordSet() ;
476                  if ( dataColKeys.isDefined( "UNIT" ) )
477                    sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
478                }
479              }
480
481              // number of rows added in this cycle
482              //os_ << "prevnr = " << prevnr << LogIO::POST ;
483              //os_ << "addednr = " << addednr << LogIO::POST ;
484              RefRows rows( prevnr, prevnr+addednr-1 ) ;
485
486              // SCANNO
487              ScalarColumn<uInt> scannoCol( table_->table(), "SCANNO" ) ;
488              Vector<uInt> scanno( addednr, scanNum ) ;
489              scannoCol.putColumnCells( rows, scanno ) ;
490              //fillId( (uInt)scanNum, "SCANNO", rows ) ;
491
492              // CYCLENO
493              ScalarColumn<uInt> cyclenoCol( table_->table(), "CYCLENO" ) ;
494              Vector<uInt> cycleno( nrow ) ;
495              indgen( cycleno, cycle ) ;
496              for ( int i = 0 ; i < addednr ; i += nrow ) {
497                Int startrow = prevnr + i ;
498                Int endrow = startrow + nrow - 1 ;
499                RefRows prows( startrow, endrow ) ;
500                cyclenoCol.putColumnCells( prows, cycleno ) ;
501              }
502              cycle += nrow ;
503
504              // BEAMNO
505              ScalarColumn<uInt> beamnoCol( table_->table(), "BEAMNO" ) ;
506              Vector<uInt> beamno( addednr, feedId ) ;
507              beamnoCol.putColumnCells( rows, beamno ) ;
508              //fillId( (uInt)feedId, "BEAMNO", rows ) ;
509
510              // IFNO
511              ScalarColumn<uInt> ifnoCol( table_->table(), "IFNO" ) ;
512              Vector<uInt> ifno( addednr, spwId ) ;
513              ifnoCol.putColumnCells( rows, ifno ) ;
514              //fillId( (uInt)spwId, "IFNO", rows ) ;
515
516              // POLNO
517              ScalarColumn<uInt> polNoCol( table_->table(), "POLNO" ) ;
518              Vector<uInt> polno( nrow ) ;
519              Int pidx = 0 ;
520              Bool crossOK = False ;
521              for ( int i = 0 ; i < npol ; i++ ) {
522                Vector<uInt> polnos = getPolNo( corrtype[i] ) ;
523                if ( polnos.size() > 1 ) {
524                  if ( crossOK ) continue ;
525                  else crossOK = True ;
526                }
527                for ( uInt j = 0 ; j < polnos.size() ; j++ ) {
528                  Int startrow = prevnr + pidx * nrow ;
529                  Int endrow = startrow + nrow - 1 ;
530                  RefRows prows( startrow, endrow ) ;
531                  polno = polnos[j] ;
532                  polNoCol.putColumnCells( prows, polno ) ;
533                  pidx++ ;
534                }
535              }
536
537              // FREQ_ID
538              ScalarColumn<uInt> freqIdCol( table_->table(), "FREQ_ID" ) ;
539              Vector<uInt> freqIds( addednr, ifmap[spwId] ) ;
540              freqIdCol.putColumnCells( rows, freqIds ) ;
541              //fillId( ifmap[spwId], "FREQ_ID", rows ) ;
542
543              // MOLECULE_ID
544              ScalarColumn<uInt> moleculeIdCol( table_->table(), "MOLECULE_ID" ) ;
545              Vector<uInt> moleculeId( addednr, molId ) ;
546              moleculeIdCol.putColumnCells( rows, moleculeId ) ;
547             
548              // REFBEAMNO
549              // set 0 at the moment
550              ScalarColumn<Int> refBeamCol( table_->table(), "REFBEAMNO" ) ;
551              Vector<Int> refBeam( addednr, 0 ) ;
552              refBeamCol.putColumnCells( rows, refBeam ) ;
553              //fillId( 0, "REFBEAMNO", rows ) ;
554             
555
556              // FLAGROW
557              ScalarColumn<uInt> flagRowCol( table_->table(), "FLAGROW" ) ;
558              ROScalarColumn<Bool> mFlagRowCol( t5, "FLAG_ROW" ) ;
559              Vector<uInt> fr( nrow ) ;
560              convertArray( fr, mFlagRowCol.getColumn() ) ;
561              for ( int i = 0 ; i < addednr ; i += nrow ) {
562                Int startrow = prevnr + i ;
563                Int endrow = startrow + nrow - 1 ;
564                RefRows prows( startrow, endrow ) ;
565                flagRowCol.putColumnCells( prows, fr ) ;
566              }
567
568              // TIME
569              MEpoch::ScalarColumn timeCol( table_->table(), "TIME" ) ;
570              MEpoch::ROScalarColumn mTimeCol( t5, "TIME" ) ;
571              Int tidx = prevnr ;
572              for ( Int i = 0 ; i < addednr ; i += nrow ) {
573                for ( Int j = 0 ; j < nrow ; j++ ) {
574                  timeCol.put( tidx++, mTimeCol( j ) ) ;
575                }
576              }
577           
578              // INTERVAL
579              ScalarColumn<Double> intervalCol( table_->table(), "INTERVAL" ) ;
580              ROScalarColumn<Double> mIntervalCol( t5, "INTERVAL" ) ;
581              Vector<Double> integ = mIntervalCol.getColumn() ;
582              for ( int i = 0 ; i < addednr ; i += nrow ) {
583                Int startrow = prevnr + i ;
584                Int endrow = startrow + nrow - 1 ;
585                RefRows prows( startrow, endrow ) ;
586                intervalCol.putColumnCells( prows, integ ) ;
587              }
588
589
590              // SRCNAME
591              ScalarColumn<String> srcNameCol( table_->table(), "SRCNAME" ) ;
592              Vector<String> vSrcName( addednr, srcName ) ;
593              srcNameCol.putColumnCells( rows, vSrcName ) ;
594
595              // SRCTYPE
596              ScalarColumn<Int> srcTypeCol( table_->table(), "SRCTYPE" ) ;
597              Vector<Int> srcType( addednr, getSrcType( stateId ) ) ;
598              srcTypeCol.putColumnCells( rows, srcType ) ;
599
600              // FIELDNAME
601              ScalarColumn<String> fieldNameCol( table_->table(), "FIELDNAME" ) ;
602              Vector<String> vFieldName( addednr, fieldName ) ;
603              fieldNameCol.putColumnCells( rows, vFieldName ) ;
604             
605              // TSYS
606              ArrayColumn<Float> tsysCol( table_->table(), "TSYS" ) ;
607              Vector<Double> sysCalTime ;
608              if ( isSysCal_ ) {
609                sysCalTime = getSysCalTime( caltabsel, mTimeCol ) ;
610                tidx = prevnr ;
611                uInt calidx = 0 ;
612                for ( Int i = 0 ; i < nrow ; i++ ) {
613                  Array<Float> tsys ;
614                  calidx = getTsys( calidx, tsys, caltabsel, sysCalTime(i) ) ;
615                  //os_ << "tsys = " << tsys << LogIO::POST ;
616                  IPosition cShape = tsys.shape() ;
617                  //os_ << "cShape = " << cShape << LogIO::POST ;
618                  if ( cShape.size() == 0 ) {
619                    cShape = IPosition( 1, npol ) ;
620                    tsys.resize( cShape ) ;
621                    tsys = 1.0 ;
622                  }
623                  for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
624                    if ( cShape.nelements() == 1 ) {
625                      Array<Float> subtsys( IPosition(1,1), tsys(IPosition(1,ipol)) ) ;
626                      tsysCol.put( prevnr+i+nrow*ipol, subtsys ) ;
627                    }
628                    else {
629                      Slice paxis( ipol, 1, 1 ) ;
630                      Slice caxis( 0, cShape[1], 1 ) ;
631                      Slicer slicer( paxis, caxis ) ;
632                      Array<Float> subtsys = tsys( slicer ) ;
633                      tsysCol.put( prevnr+i+nrow*ipol, subtsys ) ;
634                    }
635                  }                 
636                }
637              }
638              else {
639                Array<Float> tsys( IPosition( 2, 1, addednr ), 1.0 ) ;
640                tsysCol.putColumnCells( rows, tsys ) ;
641              }
642
643
644              // DIRECTION, AZIMUTH, ELEVATION, SCANRATE
645              ArrayColumn<Double> dirCol( table_->table(), "DIRECTION" ) ;
646              ScalarColumn<Float> azCol( table_->table(), "AZIMUTH" ) ;
647              ScalarColumn<Float> elCol( table_->table(), "ELEVATION" ) ;
648              ArrayColumn<Double> scanRateCol( table_->table(), "SCANRATE" ) ;
649              Vector<Double> defaultScanrate( 2, 0.0 ) ;
650              uInt diridx = 0 ;
651              MDirection::Types dirType ;
652              if ( getPt_ ) {
653                for ( Int i = 0 ; i < nrow ; i++ ) {
654                  Vector<Double> dir ;
655                  Vector<Double> scanrate ;
656                  String refString ;
657                  diridx = getDirection( diridx, dir, scanrate, refString, pointCols, mTimeCol(i).get("s").getValue() ) ;
658                  //os_ << "diridx = " << diridx << " dmTimeCol(" << i << ") = " << mTimeCol(i).get("s").getValue()-mTimeCol(0).get("s").getValue() << LogIO::POST ;
659                  //os_ << "dir = " << dir << LogIO::POST ;
660                  //os_ << "scanrate = " << scanrate << LogIO::POST ;
661                  //os_ << "refString = " << refString << LogIO::POST ;
662                  MDirection::getType( dirType, refString ) ;
663                  //os_ << "dirType = " << dirType << LogIO::POST ;
664                  mf.resetEpoch( mTimeCol(i) ) ;
665                  mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ;
666                  if ( refString == "J2000" ) {
667                    //os_ << "J2000" << LogIO::POST ;
668                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
669                      dirCol.put( prevnr+i+nrow*ipol, dir ) ;
670                    }
671                    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
672                    Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
673                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
674                      azCol.put( prevnr+i+nrow*ipol, azel(0) ) ;
675                      elCol.put( prevnr+i+nrow*ipol, azel(1) ) ;
676                    }                 
677                  }
678                  else if ( refString(0,4) == "AZEL" ) {
679                    //os_ << "AZEL" << LogIO::POST ;
680                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
681                      azCol.put( prevnr+i+nrow*ipol, dir(0) ) ;
682                      elCol.put( prevnr+i+nrow*ipol, dir(1) ) ;
683                    }
684                    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
685                    Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
686                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
687                      dirCol.put( prevnr+i+nrow*ipol, newdir ) ;
688                    }                 
689                  }
690                  else {
691                    //os_ << "OTHER: " << refString << LogIO::POST ;
692                    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
693                    Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
694                    MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
695                    Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
696                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
697                      dirCol.put( prevnr+i+nrow*ipol, newdir ) ;
698                      azCol.put( prevnr+i+nrow*ipol, dir(0) ) ;
699                      elCol.put( prevnr+i+nrow*ipol, dir(1) ) ;
700                    }                 
701                  }
702                  if ( scanrate.size() != 0 ) {
703                    //os_ << "scanrate.size() = " << scanrate.size() << LogIO::POST ;
704                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
705                      scanRateCol.put( prevnr+i+nrow*ipol, scanrate ) ;
706                    }
707                  }
708                  else {
709                    //os_ << "scanrate.size() = " << scanrate.size() << LogIO::POST ;
710                    for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
711                      scanRateCol.put( prevnr+i+nrow*ipol, defaultScanrate ) ;
712                    }
713                  }
714                }
715              }
716              else {
717                // All directions are set to source direction
718                ROArrayMeasColumn<MDirection> dmcol = pointCols.directionMeasCol() ;
719                ROArrayColumn<Double> dcol = pointCols.direction() ;
720                IPosition ip( dmcol(0).shape().nelements(), 0 ) ;
721                IPosition outp( 1, 2 ) ;
722                String ref = dmcol(0)(ip).getRefString() ;
723                Slice ds( 0, 2, 1 ) ;
724                Slice ds0( 0, 1, 1 ) ;
725                Slicer dslice0( ds, ds0 ) ;
726                Vector<Double> defaultDir = dcol(0)(dslice0).reform(outp) ;
727                MDirection::getType( dirType, "J2000" ) ;
728                mf.resetDirection( MDirection( MVDirection(srcDir), dirType ) ) ;
729                if ( ref != "J2000" ) {
730                  mf.resetEpoch( pointCols.timeMeas()(0) ) ;
731                  MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
732                  defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ;
733                }
734                for ( Int i = 0 ; i < nrow ; i++ ) {
735                  mf.resetEpoch( mTimeCol(i) ) ;
736                  for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
737                    Int localidx = prevnr+i+nrow*ipol ;
738                    MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
739                    Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ;
740                    azCol.put( localidx, azel(0) ) ;
741                    elCol.put( localidx, azel(1) ) ;
742                    dirCol.put( localidx, defaultDir ) ;
743                    scanRateCol.put( localidx, defaultScanrate ) ;
744                  }
745                }
746              }
747
748              // OPACITY
749              // not used?
750              ScalarColumn<Float> opacityCol( table_->table(), "OPACITY" ) ;
751              Vector<Float> opacity( addednr, 0.0 ) ;
752              opacityCol.putColumnCells( rows, opacity ) ;
753
754
755              // TCAL_ID
756              ScalarColumn<uInt> tcalIdCol( table_->table(), "TCAL_ID" ) ;
757              if ( isSysCal_ ) {
758                for( Int irow = 0 ; irow < nrow ; irow++ ) {
759                  Vector<uInt> tcalids = getTcalId( feedId, spwId, sysCalTime[irow] ) ;
760                  if ( tcalids.size() == 0 ) {
761                    tcalids.resize( npol ) ;
762                    tcalids = 0 ;
763                  }
764                  for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
765                  tcalIdCol.put( prevnr+irow+nrow*ipol, tcalids[ipol] ) ;
766                  }
767                }
768              }
769              else {
770                Vector<uInt> tcalid( addednr, 0 ) ;
771                tcalIdCol.putColumnCells( rows, tcalid ) ;
772              }
773
774              // FIT_ID
775              // nothing to do
776              ScalarColumn<Int> fitIdCol( table_->table(), "FIT_ID" ) ;
777              Vector<Int> fitId( addednr, -1 ) ;
778              fitIdCol.putColumnCells( rows, fitId ) ;
779
780
781              // FOCUS_ID
782              // tentative
783              ScalarColumn<uInt> focusIdCol( table_->table(), "FOCUS_ID" ) ;
784              Vector<uInt> focusId( addednr, 0 ) ;
785              focusIdCol.putColumnCells( rows, focusId ) ;
786
787
788              // WEATHER_ID
789              uInt widprev = 0 ;
790              Vector<uInt> vWid( nrow, 0 ) ;
791              if ( isWeather_ ) {
792                for ( int j = 0 ; j < nrow ; j++ ) {
793                  //os_ << "TIME value = " << mTimeCol( j ).get("s").getValue() << LogIO::POST ;
794                  uInt wid = getWeatherId( widprev, mTimeCol( j ).get("s").getValue() ) ;
795                  //os_ << "wid = " << wid << LogIO::POST ;
796                  vWid[j] = wid ;
797                  widprev = wid ;
798                }
799              }
800
801              ScalarColumn<uInt> weatherIdCol( table_->table(), "WEATHER_ID" ) ;
802              for ( int i = 0 ; i < addednr ; i += nrow ) {
803                Int startrow = prevnr + i ;
804                Int endrow = startrow + nrow - 1 ;
805                RefRows prows( startrow, endrow ) ;
806                weatherIdCol.putColumnCells( prows, vWid ) ;               
807              }
808             
809              // SRCVELOCITY, SRCPROPERMOTION and SRCDIRECTION
810              // no reference conversion for direction at the moment (assume J2000)
811              // no reference conversion for velocity at the moment (assume LSRK)
812              ArrayColumn<Double> srcPMCol( table_->table(), "SRCPROPERMOTION" ) ;
813              ArrayColumn<Double> srcDirCol( table_->table(), "SRCDIRECTION" ) ;
814              ScalarColumn<Double> srcVelCol( table_->table(), "SRCVELOCITY" ) ;
815              for ( int i = 0 ; i < addednr ; i++ ) {
816                int idx = i + prevnr ;
817                srcPMCol.put( idx, srcPM ) ;   // [rad/s]
818                srcDirCol.put( idx, srcDir ) ; // [rad]
819                srcVelCol.put( idx, sysVel ) ; // [m/s]
820              }
821
822
823              //os_ << "field: " << fieldId << " scan: " << scanNum << " obs: " << obsId << " state: " << stateId << " ddid: " << ddId << endl ;
824              //os_ << "t.nrow() = " << t5.nrow() << endl ;
825              totalrow += t5.nrow() ;
826              //os_ << "totalrow = " << totalrow << LogIO::POST ;
827              iter5.next() ;
828            }
829            iter4.next() ;
830          }
831          iter3.next() ;
832        }
833        iter2.next() ;
834      }
835      iter1.next() ;
836    }
837    if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
838    iter0.next() ;
839  }
840
841  // Keywords
842  sdh.nif = ifmap.size() ;
843  String antennaName = antCols.name()(antenna_) ;
844  if ( antennaName == telescopeName ) {
845    sdh.antennaname = antennaName ;
846  }
847  else {
848    sdh.antennaname = telescopeName + "//" + antennaName ;
849  }
850  if ( stationName != "" ) {
851    sdh.antennaname += "@" + stationName ;
852  }
853  sdh.antennaposition = antCols.position()(antenna_);
854  ROMSPointingColumns pointingCols( mstable_.pointing() ) ;
855  String dirref = pointingCols.direction().keywordSet().asRecord("MEASINFO").asString("Ref") ;
856  if ( dirref == "AZELGEO" || dirref == "AZEL" ) {
857    dirref = "J2000" ;
858  }
859  sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
860  sdh.epoch = "UTC" ;
861  if (sdh.freqref == "TOPO") {
862    sdh.freqref = "TOPOCENT";
863  } else if (sdh.freqref == "GEO") {
864    sdh.freqref = "GEOCENTR";
865  } else if (sdh.freqref == "BARY") {
866    sdh.freqref = "BARYCENT";
867  } else if (sdh.freqref == "GALACTO") {
868    sdh.freqref = "GALACTOC";
869  } else if (sdh.freqref == "LGROUP") {
870    sdh.freqref = "LOCALGRP";
871  } else if (sdh.freqref == "CMB") {
872    sdh.freqref = "CMBDIPOL";
873  } else if (sdh.freqref == "REST") {
874    sdh.freqref = "SOURCE";
875  }
876  table_->setHeader( sdh ) ;
877}
878
879void MSFiller::close()
880{
881  tablesel_.closeSubTables() ;
882  mstable_.closeSubTables() ;
883  tablesel_.unlock() ;
884  mstable_.unlock() ;
885}
886
887void MSFiller::fillId( uInt idx, const char *colname, RefRows &rows )
888{
889  ScalarColumn<uInt> col( table_->table(), colname ) ;
890  Vector<uInt> ids( rows.nrow(), idx ) ;
891  col.putColumnCells( rows, ids ) ;
892}
893
894void MSFiller::fillId( Int idx, const char *colname, RefRows &rows )
895{
896  ScalarColumn<Int> col( table_->table(), colname ) ;
897  Vector<Int> ids( rows.nrow(), idx ) ;
898  col.putColumnCells( rows, ids ) ;
899}
900
901Int MSFiller::getSrcType( Int stateId )
902{
903  MSState statetab = mstable_.state() ;
904  ROScalarColumn<String> obsModeCol( statetab, "OBS_MODE" ) ;
905  String obsMode = obsModeCol( stateId ) ;
906  ROScalarColumn<Bool> sigCol( statetab, "SIG" ) ;
907  ROScalarColumn<Bool> refCol( statetab, "REF" ) ;
908  Bool sig = sigCol( stateId ) ;
909  Bool ref = refCol( stateId ) ;
910  //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
911
912  // determine separator
913  String sep = "" ;
914  if ( obsMode.find( ":" ) != String::npos ) {
915    sep = ":" ;
916  }
917  else if ( obsMode.find( "." ) != String::npos ) {
918    sep = "." ;
919  }
920 
921  // determine SRCTYPE
922  Int srcType = SrcType::NOTYPE ;
923  if ( sep == ":" ) {
924    // sep == ":"
925    //
926    // GBT case
927    //
928    // obsMode1=Nod
929    //    NOD
930    // obsMode1=OffOn
931    //    obsMode2=PSWITCHON:  PSON
932    //    obsMode2=PSWITCHOFF: PSOFF
933    // obsMode1=??
934    //    obsMode2=FSWITCH:
935    //       SIG=1: FSON
936    //       REF=1: FSOFF
937    Int epos = obsMode.find_first_of( sep ) ;
938    Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
939    String obsMode1 = obsMode.substr( 0, epos ) ;
940    String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
941    if ( obsMode1 == "Nod" ) {
942      srcType = SrcType::NOD ;
943    }
944    else if ( obsMode1 == "OffOn" ) {
945      if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
946      if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
947    }
948    else {
949      if ( obsMode2 == "FSWITCH" ) {
950        if ( sig ) srcType = SrcType::FSON ;
951        if ( ref ) srcType = SrcType::FSOFF ;
952      }
953    }
954  }
955  else if ( sep == "." ) {
956    // sep == "."
957    //
958    // ALMA & EVLA case (MS via ASDM)
959    //
960    // obsMode1=CALIBRATE_*
961    //    obsMode2=ON_SOURCE: PONCAL
962    //    obsMode2=OFF_SOURCE: POFFCAL
963    // obsMode1=OBSERVE_TARGET
964    //    obsMode2=ON_SOURCE: PON
965    //    obsMode2=OFF_SOURCE: POFF
966    Int epos = obsMode.find_first_of( sep ) ;
967    Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
968    String obsMode1 = obsMode.substr( 0, epos ) ;
969    String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
970    if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
971      if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
972      if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
973    }
974    else if ( obsMode1 == "OBSERVE_TARGET" ) {
975      if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
976      if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
977    }
978  }
979  else {
980    if ( sig ) srcType = SrcType::SIG ;
981    if ( ref ) srcType = SrcType::REF ;
982  }
983   
984  //os_ << "srcType = " << srcType << LogIO::POST ;
985
986  return srcType ;
987}
988
989Vector<uInt> MSFiller::getPolNo( Int corrType )
990{
991  Vector<uInt> polno( 1 ) ;
992
993  if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {
994    polno = 0 ;
995  }
996  else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {
997    polno = 1 ;
998  }
999  else if ( corrType == Stokes::U ) {
1000    polno = 2 ;
1001  }
1002  else if ( corrType == Stokes::V ) {
1003    polno = 3 ;
1004  }
1005  else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {
1006    polno.resize( 2 ) ;
1007    polno[0] = 2 ;
1008    polno[1] = 3 ;
1009  }
1010  else if ( corrType == Stokes::Plinear ) {
1011    polno[0] = 1 ;
1012  }
1013  else if ( corrType == Stokes::Pangle ) {
1014    polno[0] = 2 ;
1015  }
1016  else {
1017    polno = 99 ;
1018  }
1019  //os_ << "polno = " << polno << LogIO::POST ;
1020 
1021  return polno ;
1022}
1023
1024String MSFiller::getPolType( Int corrType )
1025{
1026  String poltype = "" ;
1027
1028  if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )
1029    poltype = "stokes" ;
1030  else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )
1031    poltype = "linear" ;
1032  else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )
1033    poltype = "circular" ;
1034  else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )
1035    poltype = "linpol" ;
1036
1037  return poltype ;
1038}
1039
1040void MSFiller::fillWeather()
1041{
1042  MSWeather mWeather( mstable_.weather() ) ;
1043  MSWeather mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1044  //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1045  if ( mWeatherSel.nrow() == 0 ) {
1046    os_ << "No rows with ANTENNA_ID = " << antenna_ << ", Try -1..." << LogIO::POST ;
1047    mWeatherSel = MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ;
1048    if ( mWeatherSel.nrow() == 0 ) {
1049      os_ << "No rows in WEATHER table" << LogIO::POST ;
1050    }
1051  }
1052  ROMSWeatherColumns mWeatherCols( mWeatherSel ) ;
1053  Int wnrow = mWeatherCols.nrow() ;
1054  //os_ << "wnrow = " << wnrow << LogIO::POST ;
1055
1056  if ( wnrow == 0 )
1057    return ;
1058
1059  Table wtab = table_->weather().table() ;
1060  wtab.addRow( wnrow ) ;
1061
1062  ScalarColumn<Float> tempCol( wtab, "TEMPERATURE" ) ;
1063  tempCol.putColumn( mWeatherCols.temperature() ) ;
1064  ScalarColumn<Float> pressCol( wtab, "PRESSURE" ) ;
1065  pressCol.putColumn( mWeatherCols.pressure() ) ;
1066  ScalarColumn<Float> humCol( wtab, "HUMIDITY" ) ;
1067  humCol.putColumn( mWeatherCols.relHumidity() ) ;
1068  ScalarColumn<Float> windVelCol( wtab, "WINDSPEED" ) ;
1069  windVelCol.putColumn( mWeatherCols.windSpeed() ) ;
1070  ScalarColumn<Float> windDirCol( wtab, "WINDAZ" ) ;
1071  windDirCol.putColumn( mWeatherCols.windDirection() ) ;
1072  Vector<uInt> ids( wnrow ) ;
1073  indgen( ids ) ;
1074  ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1075  idCol.putColumn( ids ) ;
1076
1077  String tUnit = mWeatherCols.timeQuant().getUnits() ;
1078  mwTime_ = mWeatherCols.time().getColumn() ;
1079  if ( tUnit == "d" )
1080    mwTime_ *= 86400.0 ;
1081  String iUnit = mWeatherCols.intervalQuant().getUnits() ;
1082  mwInterval_ = mWeatherCols.interval().getColumn() ;
1083  if ( iUnit == "d" )
1084    mwInterval_ *= 86400.0 ;
1085  //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1086}
1087
1088void MSFiller::fillFocus()
1089{
1090  // tentative
1091  Table tab = table_->focus().table() ;
1092  tab.addRow( 1 ) ;
1093  ScalarColumn<uInt> idCol( tab, "ID" ) ;
1094  idCol.put( 0, 0 ) ;
1095}
1096
1097void MSFiller::fillTcal()
1098{
1099  MSSysCal sctab = mstable_.sysCal() ;
1100  if ( sctab.nrow() == 0 ) {
1101    os_ << "No SysCal rows" << LogIO::POST ;
1102    return ;
1103  }
1104  Bool isSp = sctab.tableDesc().isColumn( "TCAL_SPECTRUM" ) ;
1105  MSSysCal sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
1106  if ( sctabsel.nrow() == 0 ) {
1107    os_ << "No SysCal rows" << LogIO::POST ;
1108    return ;
1109  }
1110  ROArrayColumn<Float> tmpTcalCol( sctabsel, "TCAL" ) ;
1111  uInt npol = tmpTcalCol.shape( 0 )(0) ;
1112  //os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
1113  Table tab = table_->tcal().table() ;
1114  ScalarColumn<uInt> idCol( tab, "ID" ) ;
1115  ScalarColumn<String> timeCol( tab, "TIME" ) ;
1116  ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
1117  uInt oldnr = 0 ;
1118  uInt newnr = 0 ;
1119  TableIterator iter0( sctabsel, "FEED_ID" ) ;
1120  // Record for TCAL_ID
1121  // "FIELD0": "SPW0": Vector<uInt>
1122  //           "SPW1": Vector<uInt>
1123  //  ...
1124  while( !iter0.pastEnd() ) {
1125    MSSysCal t0( iter0.table() ) ;
1126    ROScalarColumn<Int> feedIdCol( t0, "FEED_ID" ) ;
1127    Int feedId = feedIdCol( 0 ) ;
1128    String ffield = "FEED" + String::toString( feedId ) ;
1129    Record rec ;
1130    TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
1131    while( !iter1.pastEnd() ) {
1132      MSSysCal t1( iter1.table().sort("TIME") ) ;
1133      uInt nrow = t1.nrow() ;
1134      ROMSSysCalColumns scCols( t1 ) ;
1135      Int spwId = scCols.spectralWindowId()(0) ;
1136      String spwfield = "SPW" + String::toString( spwId ) ;
1137      ROScalarQuantColumn<Double> scTimeCol = scCols.timeQuant() ;
1138      ROArrayColumn<Float> scTcalCol ;
1139      IPosition newShape( 2, 1, nrow ) ;
1140      if ( isSp ) {
1141        scTcalCol.reference( scCols.tcalSpectrum() ) ;
1142        newShape[0] = scTcalCol.shape(0)(1) ;
1143      }
1144      else {
1145        scTcalCol.reference( scCols.tcal() ) ;
1146      }
1147      Vector<uInt> idx( nrow ) ;
1148      Vector<String> sTime( nrow ) ;
1149      for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
1150        sTime[irow] = MVTime( scTimeCol(irow) ).string(MVTime::YMD) ;
1151      }
1152      Vector<uInt> idminmax( 2, oldnr ) ;
1153      for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1154        tab.addRow( nrow ) ;
1155        newnr += nrow ;
1156        RefRows rows( oldnr, newnr-1 ) ;
1157        indgen( idx, oldnr ) ;
1158        idCol.putColumnCells( rows, idx ) ;
1159        timeCol.putColumnCells( rows, sTime ) ;
1160        Slicer slicer ;
1161        if ( isSp ) {
1162          Slice paxis( ipol, 1, 1 ) ;
1163          Slice caxis( 0, newShape[0], 1 ) ;
1164          slicer = Slicer( paxis, caxis ) ;
1165        }
1166        else {
1167          Slice paxis( ipol, 1, 1 ) ;
1168          slicer = Slicer( paxis ) ;
1169        }
1170        Array<Float> subtcal = scTcalCol.getColumn( slicer ).reform( newShape ) ;
1171        tcalCol.putColumnCells( rows, subtcal ) ;
1172        oldnr += nrow ;
1173      }
1174      idminmax[1] = newnr - 1 ;
1175      rec.define( spwfield, idminmax ) ;
1176      iter1++ ;
1177    }
1178    tcalrec_.defineRecord( ffield, rec ) ;
1179    iter0++ ;
1180  }
1181
1182  //tcalrec_.print( std::cout ) ;
1183}
1184
1185// void MSFiller::fillMolecules()
1186// {
1187//   os_ << "MSFiller::fillMolecules()" << LogIO::POST ;
1188//   // tentative
1189//   Table tab = table_->molecules().table() ;
1190//   tab.addRow( 1 ) ;
1191//   ScalarColumn<uInt> idCol( tab, "ID" ) ;
1192//   idCol.put( 0, 0 ) ;
1193// }
1194
1195// void MSFiller::fillFit()
1196// {
1197//   os_ << "MSFiller::fillFit()" << LogIO::POST ;
1198//   // tentative
1199//   Table tab = table_->fit().table() ;
1200//   tab.addRow( 1 ) ;
1201//   ScalarColumn<uInt> idCol( tab, "ID" ) ;
1202//   idCol.put( 0, 0 ) ;
1203// }
1204
1205// void MSFiller::fillFrequencies()
1206// {
1207//   os_ << "MSFiller::fillFrequencies()" << LogIO::POST ;
1208//   // tentative
1209//   Table tab = table_->frequencies().table() ;
1210//   tab.addRow( 1 ) ;
1211//   ScalarColumn<uInt> idCol( tab, "ID" ) ;
1212//   idCol.put( 0, 0 ) ;
1213// }
1214
1215// void MSFiller::fillHistory()
1216// {
1217//   os_ << "MSFiller::fillHistory()" << LogIO::POST ;
1218//   // tentative
1219//   Table tab = table_->history().table() ;
1220//   tab.addRow( 1 ) ;
1221//   ScalarColumn<uInt> idCol( tab, "ID" ) ;
1222//   idCol.put( 0, 0 ) ;
1223// }
1224
1225uInt MSFiller::getWeatherId( uInt idx, Double wtime )
1226{
1227  uInt nrow = mwTime_.size() ;
1228  if ( nrow == 0 )
1229    return 0 ;
1230  uInt wid = nrow ;
1231  for ( uInt i = idx ; i < nrow-1 ; i++ ) {
1232    Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
1233    // use of INTERVAL column is problematic
1234    // since there are "blank" time of weather monitoring
1235    //Double tEnd = tStart + mwInterval_[i] ;
1236    Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
1237    //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1238    if ( wtime >= tStart && wtime <= tEnd ) {
1239      wid = i ;
1240      break ;
1241    }
1242  }
1243  if ( wid == nrow ) {
1244    uInt i = nrow - 1 ;
1245    Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
1246    Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
1247    //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1248    if ( wtime >= tStart && wtime <= tEnd )
1249      wid = i ;
1250  }
1251
1252  //if ( wid == nrow )
1253  //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
1254
1255  return wid ;
1256}
1257
1258Vector<Double> MSFiller::getSysCalTime( MSSysCal &tab, MEpoch::ROScalarColumn &tcol )
1259{
1260  uInt nrow = tcol.table().nrow() ;
1261  Vector<Double> tstr( nrow, -1.0 ) ;
1262  if ( tab.nrow() == 0 )
1263    return tstr ;
1264  uInt scnrow = tab.nrow() ;
1265  ROMSSysCalColumns sysCalCols( tab ) ;
1266  ROScalarMeasColumn<MEpoch> scTimeCol = sysCalCols.timeMeas() ;
1267  ROScalarQuantColumn<Double> scIntervalCol = sysCalCols.intervalQuant() ;
1268  uInt idx = 0 ;
1269  const Double half = 0.5e0 ;
1270  for ( uInt i = 0 ; i < nrow ; i++ ) {
1271    Double t = tcol( i ).get( "s" ).getValue() ;
1272    for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
1273      Double tsc1 = scTimeCol( j ).get( "s" ).getValue() ;
1274      Double dt1 = scIntervalCol( j ).getValue("s") ;
1275      Double tsc2 = scTimeCol( j+1 ).get( "s" ).getValue() ;
1276      Double dt2 = scIntervalCol( j+1 ).getValue("s") ;
1277      if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {
1278        tstr[i] = tsc1 ;
1279        idx = j ;
1280        break ;
1281      }
1282    }
1283    if ( tstr[i] == -1.0 ) {
1284      Double tsc = scTimeCol( scnrow-1 ).get( "s" ).getValue() ;
1285      Double dt = scIntervalCol( scnrow-1 ).getValue( "s" ) ;
1286      if ( t <= tsc+0.5*dt )
1287        tstr[i] = tsc ;
1288    }
1289  }
1290  return tstr ;
1291}
1292
1293uInt MSFiller::getTsys( uInt idx, Array<Float> &tsys, MSSysCal &tab, Double t )
1294{
1295  uInt nrow = tab.nrow() ;
1296  if ( nrow == 0 ) {
1297    os_ << "No SysCal rows" << LogIO::POST ;
1298    tsys.resize( IPosition(0) ) ;
1299    return 0 ;
1300  }
1301  Bool isSp = tab.tableDesc().isColumn( "TSYS_SPECTRUM" ) ;
1302  ROMSSysCalColumns calCols( tab ) ;
1303  ROScalarMeasColumn<MEpoch> scTimeCol = calCols.timeMeas() ;
1304  ROArrayColumn<Float> mTsysCol ;
1305  if ( isSp ) {
1306    mTsysCol.reference( calCols.tsysSpectrum() ) ;
1307  }
1308  else {
1309    mTsysCol.reference( calCols.tsys() ) ;
1310  }
1311  for ( uInt i = idx ; i < nrow ; i++ ) {
1312    Double tref = scTimeCol( i ).get( "s" ).getValue() ;
1313    if ( t == tref ) {
1314      tsys.reference( mTsysCol( i ) ) ;
1315      idx = i ;
1316      break ;
1317    }
1318  }
1319  return idx ;
1320}
1321
1322Vector<uInt> MSFiller::getTcalId( Int fid, Int spwid, Double t )
1323{
1324  String feed = "FEED" + String::toString(fid) ;
1325  String spw = "SPW" + String::toString(spwid) ;
1326  String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
1327  Table ttab = table_->tcal().table() ;
1328  if ( ttab.nrow() == 0 ) {
1329    os_ << "No TCAL rows" << LogIO::POST ;
1330    Vector<uInt> tcalids( 0 ) ;
1331    return  tcalids ;
1332  }
1333  Vector<uInt> ids = tcalrec_.asRecord(feed).asArrayuInt(spw) ;
1334  Table ttabsel = ttab( ttab.col("TIME") == sctime && ttab.col("ID") >= ids[0] && ttab.col("ID") <= ids[1] ).sort("ID") ;
1335  uInt nrow = ttabsel.nrow() ;
1336  Vector<uInt> tcalids( nrow ) ;
1337  if ( nrow == 0 ) {
1338    os_ << "No TCAL rows" << LogIO::POST ;
1339    return tcalids ;
1340  }
1341  ROScalarColumn<uInt> idCol( ttabsel, "ID" ) ;
1342  tcalids[0] = idCol(0) ;
1343  if ( nrow == 2 ) {
1344    tcalids[1] = idCol(1) ;
1345  }
1346  else if ( nrow == 3 ) {
1347    tcalids[1] = idCol(2) ;
1348    tcalids[2] = idCol(1) ;
1349  }
1350  else if ( nrow == 4 ) {
1351    tcalids[1] = idCol(3) ;
1352    tcalids[2] = idCol(1) ;
1353    tcalids[3] = idCol(2) ;
1354  }
1355 
1356  return tcalids ;
1357}
1358
1359uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, ROMSPointingColumns &cols, Double t )
1360{
1361  // assume that cols is sorted by TIME
1362  Bool doInterp = False ;
1363  uInt nrow = cols.nrow() ;
1364  if ( nrow == 0 )
1365    return 0 ;
1366  ROScalarMeasColumn<MEpoch> tcol = cols.timeMeas() ;
1367  ROArrayMeasColumn<MDirection> dmcol = cols.directionMeasCol() ;
1368  ROArrayColumn<Double> dcol = cols.direction() ;
1369  // ensure that tcol(idx) < t
1370  //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
1371  while ( tcol(idx).get("s").getValue() > t && idx > 0 )
1372    idx-- ;
1373  //os_ << "idx = " << idx << LogIO::POST ;
1374
1375  // index search
1376  for ( uInt i = idx ; i < nrow ; i++ ) {
1377    Double tref = tcol( i ).get( "s" ).getValue() ;
1378    if ( tref == t ) {
1379      idx = i ;
1380      break ;
1381    }
1382    else if ( tref > t ) {
1383      if ( i == 0 ) {
1384        idx = i ;
1385      }
1386      else {
1387        idx = i-1 ;
1388        doInterp = True ;
1389      }
1390      break ;
1391    }
1392    else {
1393      idx = nrow - 1 ;
1394    }
1395  }
1396  //os_ << "searched idx = " << idx << LogIO::POST ;
1397
1398  Slice ds( 0, 2, 1 ) ;
1399  Slice ds0( 0, 1, 1 ) ;
1400  Slice ds1( 1, 1, 1 ) ;
1401  Slicer dslice0( ds, ds0 ) ;
1402  Slicer dslice1( ds, ds1 ) ;
1403  //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
1404  IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
1405  //os_ << "ip = " << ip << LogIO::POST ;
1406  ref = dmcol(idx)(ip).getRefString() ;
1407  //os_ << "ref = " << ref << LogIO::POST ;
1408  IPosition outp(1,2) ;
1409  if ( doInterp ) {
1410    //os_ << "do interpolation" << LogIO::POST ;
1411    //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
1412    Double tref0 = tcol(idx).get("s").getValue() ;
1413    Double tref1 = tcol(idx+1).get("s").getValue() ;
1414    Vector<Double> dir0 = dcol(idx)(dslice0).reform(outp) ;
1415    //os_ << "dir0 = " << dir0 << LogIO::POST ;
1416    Vector<Double> dir1 = dcol(idx+1)(dslice0).reform(outp) ;
1417    //os_ << "dir1 = " << dir1 << LogIO::POST ;
1418    Double dt0 = t - tref0 ;
1419    Double dt1 = tref1 - t ;
1420    dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;
1421    if ( dcol(idx).shape()(1) > 1 ) {
1422      if ( dt0 >= dt1 ) {
1423        srate.reference( dcol(idx)(dslice1).reform(outp) ) ;
1424      }
1425      else {
1426        srate.reference( dcol(idx+1)(dslice1) ) ;
1427      }
1428    }
1429    //os_ << "dir = " << dir << LogIO::POST ;
1430  }
1431  else {
1432    //os_ << "no interpolation" << LogIO::POST ;
1433    dir.reference( dcol(idx)(dslice0).reform(outp) ) ;
1434    if ( dcol(idx).shape()(1) > 1 ) {
1435      srate.reference( dcol(idx)(dslice1).reform(outp) ) ;
1436    }
1437  }
1438
1439  return idx ;
1440}
1441
1442} ;
1443
Note: See TracBrowser for help on using the repository browser.