source: trunk/src/MSFiller.cpp@ 1979

Last change on this file since 1979 was 1978, checked in by Takeshi Nakazato, 14 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.