source: trunk/src/MSFiller.cpp@ 1974

Last change on this file since 1974 was 1974, checked in by Takeshi Nakazato, 15 years ago

New Development: Yes

JIRA Issue: Yes CAS-2718

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: New class msfiller and mswriter added

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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