source: trunk/src/MSFiller.cpp @ 1978

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

New Development: No

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

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