source: trunk/src/MSFiller.cpp@ 1989

Last change on this file since 1989 was 1988, 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...

More speed up filler.
File size: 61.9 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/TableColumn.h>
20#include <tables/Tables/ScalarColumn.h>
21#include <tables/Tables/ArrayColumn.h>
22#include <tables/Tables/RefRows.h>
23#include <tables/Tables/TableParse.h>
24#include <tables/Tables/RefRows.h>
25
26#include <casa/Containers/Block.h>
27#include <casa/Logging/LogIO.h>
28#include <casa/Arrays/Slicer.h>
29#include <casa/Quanta/MVTime.h>
30#include <casa/OS/Path.h>
31
32#include <measures/Measures/Stokes.h>
33#include <measures/Measures/MEpoch.h>
34#include <measures/Measures/MCEpoch.h>
35#include <measures/Measures/MFrequency.h>
36#include <measures/Measures/MCFrequency.h>
37#include <measures/Measures/MPosition.h>
38#include <measures/Measures/MCPosition.h>
39#include <measures/Measures/MDirection.h>
40#include <measures/Measures/MCDirection.h>
41#include <measures/Measures/MeasConvert.h>
42#include <measures/TableMeasures/ScalarMeasColumn.h>
43#include <measures/TableMeasures/ArrayMeasColumn.h>
44#include <measures/TableMeasures/ScalarQuantColumn.h>
45#include <measures/TableMeasures/ArrayQuantColumn.h>
46
47#include <atnf/PKSIO/SrcType.h>
48
49#include "MSFiller.h"
50#include "STHeader.h"
51
52#include <ctime>
53#include <sys/time.h>
54
55double gettimeofday_sec()
56{
57 struct timeval tv ;
58 gettimeofday( &tv, NULL ) ;
59 return tv.tv_sec + (double)tv.tv_usec*1.0e-6 ;
60}
61
62using namespace casa ;
63using namespace std ;
64
65namespace asap {
66MSFiller::MSFiller( casa::CountedPtr<Scantable> stable )
67 : table_( stable ),
68 tablename_( "" ),
69 antenna_( -1 ),
70 getPt_( False ),
71 isFloatData_( False ),
72 isData_( False ),
73 isDoppler_( False ),
74 isFlagCmd_( False ),
75 isFreqOffset_( False ),
76 isHistory_( False ),
77 isProcessor_( False ),
78 isSysCal_( False ),
79 isWeather_( False )
80{
81 os_ = LogIO() ;
82 os_.origin( LogOrigin( "MSFiller", "MSFiller()", WHERE ) ) ;
83}
84
85MSFiller::~MSFiller()
86{
87 os_.origin( LogOrigin( "MSFiller", "~MSFiller()", WHERE ) ) ;
88}
89
90bool MSFiller::open( const std::string &filename, const casa::Record &rec )
91{
92 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ;
93 //double startSec = gettimeofday_sec() ;
94 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;
95 //os_ << " filename = " << filename << endl ;
96
97 // parsing MS options
98 if ( rec.isDefined( "ms" ) ) {
99 Record msrec = rec.asRecord( "ms" ) ;
100 if ( msrec.isDefined( "getpt" ) ) {
101 getPt_ = msrec.asBool( "getpt" ) ;
102 }
103 if ( msrec.isDefined( "antenna" ) ) {
104 if ( msrec.type( msrec.fieldNumber( "antenna" ) ) == TpInt ) {
105 antenna_ = msrec.asInt( "antenna" ) ;
106 }
107 else {
108 antenna_ = atoi( msrec.asString( "antenna" ).c_str() ) ;
109 }
110 }
111 else {
112 antenna_ = 0 ;
113 }
114 }
115
116 os_ << "Parsing MS options" << endl ;
117 os_ << " getPt = " << getPt_ << endl ;
118 os_ << " antenna = " << antenna_ << LogIO::POST ;
119
120 MeasurementSet *tmpMS = new MeasurementSet( filename, Table::Old ) ;
121 //mstable_ = (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
122 // && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ;
123 tablename_ = tmpMS->tableName() ;
124 mstable_ = MeasurementSet( (*tmpMS)( tmpMS->col("ANTENNA1") == antenna_
125 && tmpMS->col("ANTENNA1") == tmpMS->col("ANTENNA2") ) ) ;
126// stringstream ss ;
127// ss << "SELECT FROM $1 WHERE ANTENNA1 == ANTENNA2 && ANTENNA1 == " << antenna_ ;
128// String taql( ss.str() ) ;
129// mstable_ = MeasurementSet( tableCommand( taql, *tmpMS ) ) ;
130 delete tmpMS ;
131
132 // check which data column exists
133 isFloatData_ = mstable_.tableDesc().isColumn( "FLOAT_DATA" ) ;
134 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ;
135
136 //double endSec = gettimeofday_sec() ;
137 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
138 return true ;
139}
140
141void MSFiller::fill()
142{
143 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ;
144 //double startSec = gettimeofday_sec() ;
145 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;
146
147 // Initialize header
148 STHeader sdh ;
149 sdh.nchan = 0 ;
150 sdh.npol = 0 ;
151 sdh.nif = 0 ;
152 sdh.nbeam = 0 ;
153 sdh.observer = "" ;
154 sdh.project = "" ;
155 sdh.obstype = "" ;
156 sdh.antennaname = "" ;
157 sdh.antennaposition.resize( 0 ) ;
158 sdh.equinox = 0.0 ;
159 sdh.freqref = "" ;
160 sdh.reffreq = -1.0 ;
161 sdh.bandwidth = 0.0 ;
162 sdh.utc = 0.0 ;
163 sdh.fluxunit = "" ;
164 sdh.epoch = "" ;
165 sdh.poltype = "" ;
166
167 // check if optional table exists
168 //const TableRecord msrec = tablesel_.keywordSet() ;
169 const TableRecord msrec = mstable_.keywordSet() ;
170 isDoppler_ = msrec.isDefined( "DOPPLER" ) ;
171 isFlagCmd_ = msrec.isDefined( "FLAG_CMD" ) ;
172 isFreqOffset_ = msrec.isDefined( "FREQ_OFFSET" ) ;
173 isHistory_ = msrec.isDefined( "HISTORY" ) ;
174 isProcessor_ = msrec.isDefined( "PROCESSOR" ) ;
175 isSysCal_ = msrec.isDefined( "SYSCAL" ) ;
176 isWeather_ = msrec.isDefined( "WEATHER" ) ;
177
178 // Access to MS subtables
179 MSField fieldtab = mstable_.field() ;
180 MSPolarization poltab = mstable_.polarization() ;
181 MSDataDescription ddtab = mstable_.dataDescription() ;
182 MSObservation obstab = mstable_.observation() ;
183 MSSource srctab = mstable_.source() ;
184 MSSpectralWindow spwtab = mstable_.spectralWindow() ;
185 MSSysCal caltab = mstable_.sysCal() ;
186 if ( caltab.nrow() == 0 )
187 isSysCal_ = False ;
188 MSPointing pointtab = mstable_.pointing() ;
189 if ( mstable_.weather().nrow() == 0 )
190 isWeather_ = False ;
191 MSState stattab = mstable_.state() ;
192 MSAntenna anttab = mstable_.antenna() ;
193
194 // TEST
195 // memory allocation by boost::object_pool
196 boost::object_pool<ROTableColumn> *tpoolr = new boost::object_pool<ROTableColumn> ;
197 boost::object_pool<TableColumn> *tpoolw = new boost::object_pool<TableColumn> ;
198 //
199
200 // SUBTABLES: FREQUENCIES
201 table_->frequencies().setFrame( "LSRK" ) ;
202 table_->frequencies().setFrame( "LSRK", True ) ;
203
204 // SUBTABLES: WEATHER
205 if ( isWeather_ )
206 fillWeather() ;
207
208 // SUBTABLES: FOCUS
209 fillFocus() ;
210
211 // SUBTABLES: TCAL
212 if ( isSysCal_ )
213 fillTcal( tpoolr, tpoolw ) ;
214
215 // SUBTABLES: FIT
216 //fillFit() ;
217
218 // SUBTABLES: HISTORY
219 //fillHistory() ;
220
221 // shared pointers
222 ROTableColumn *tcolr ;
223 TableColumn *tcolw ;
224
225 // Scantable columns
226 Table stab = table_->table() ;
227 TableColumn *scannoCol = tpoolw->construct( stab, "SCANNO" ) ;
228 TableColumn *cyclenoCol = tpoolw->construct( stab, "CYCLENO" ) ;
229 TableColumn *beamnoCol = tpoolw->construct( stab, "BEAMNO" ) ;
230 TableColumn *ifnoCol = tpoolw->construct( stab, "IFNO" ) ;
231 TableColumn *polnoCol = tpoolw->construct( stab, "POLNO" ) ;
232 TableColumn *freqidCol = tpoolw->construct( stab, "FREQ_ID" ) ;
233 TableColumn *molidCol = tpoolw->construct( stab, "MOLECULE_ID" ) ;
234 TableColumn *flagrowCol = tpoolw->construct( stab, "FLAGROW" ) ;
235 ScalarMeasColumn<MEpoch> *timeCol = new ScalarMeasColumn<MEpoch>( stab, "TIME" ) ;
236 TableColumn *intervalCol = tpoolw->construct( stab, "INTERVAL" ) ;
237 TableColumn *srcnameCol = tpoolw->construct( stab, "SRCNAME" ) ;
238 TableColumn *srctypeCol = tpoolw->construct( stab, "SRCTYPE" ) ;
239 TableColumn *fieldnameCol = tpoolw->construct( stab, "SRCNAME" ) ;
240 ArrayColumn<Float> *spCol = new ArrayColumn<Float>( stab, "SPECTRA" ) ;
241 ArrayColumn<uChar> *flCol = new ArrayColumn<uChar>( stab, "FLAGTRA" ) ;
242 ArrayColumn<Float> *tsysCol = new ArrayColumn<Float>( stab, "TSYS" ) ;
243 ArrayColumn<Double> *dirCol = new ArrayColumn<Double>( stab, "DIRECTION" ) ;
244 TableColumn *azCol = tpoolw->construct( stab, "AZIMUTH" ) ;
245 TableColumn *elCol = tpoolw->construct( stab, "ELEVATION" ) ;
246 TableColumn *tcalidCol = tpoolw->construct( stab, "TCAL_ID" ) ;
247 TableColumn *focusidCol = tpoolw->construct( stab, "FOCUS_ID" ) ;
248 TableColumn *weatheridCol = tpoolw->construct( stab, "WEATHER_ID" ) ;
249 ArrayColumn<Double> *srcpmCol = new ArrayColumn<Double>( stab, "SRCPROPERMOTION" ) ;
250 ArrayColumn<Double> *srcdirCol = new ArrayColumn<Double>( stab, "SRCDIRECTION" ) ;
251 TableColumn *srcvelCol = tpoolw->construct( stab, "SRCVELOCITY" ) ;
252 ArrayColumn<Double> *scanrateCol = new ArrayColumn<Double>( stab, "SCANRATE" ) ;
253
254 // MAIN
255 // Iterate over several ids
256 Int oldnr = table_->nrow() ;
257 map<Int, uInt> ifmap ; // (IFNO, FREQ_ID) pair
258 ROArrayQuantColumn<Double> *sharedQDArrCol = new ROArrayQuantColumn<Double>( anttab, "POSITION" ) ;
259 Vector< Quantum<Double> > antpos = (*sharedQDArrCol)( antenna_ ) ;
260 delete sharedQDArrCol ;
261 MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ;
262 //pointtab = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;
263 pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ;
264 tcolr = tpoolr->construct( anttab, "STATION" ) ;
265 String stationName = tcolr->asString( antenna_ ) ;
266 tpoolr->destroy( tcolr ) ;
267 tcolr = tpoolr->construct( anttab, "NAME" ) ;
268 String antennaName = tcolr->asString( antenna_ ) ;
269 tpoolr->destroy( tcolr ) ;
270 sdh.antennaposition.resize( 3 ) ;
271 for ( int i = 0 ; i < 3 ; i++ )
272 sdh.antennaposition[i] = antpos[i].getValue( "m" ) ;
273 String telescopeName = "" ;
274
275 //
276 // ITERATION: OBSERVATION_ID
277 //
278 Int added0 = 0 ;
279 Int current0 = table_->nrow() ;
280 TableIterator iter0( mstable_, "OBSERVATION_ID" ) ;
281 while( !iter0.pastEnd() ) {
282 Table t0 = iter0.table() ;
283 tcolr = tpoolr->construct( t0, "OBSERVATION_ID" ) ;
284 Int obsId = tcolr->asInt( 0 ) ;
285 tpoolr->destroy( tcolr ) ;
286 if ( sdh.observer == "" ) {
287 tcolr = tpoolr->construct( obstab, "OBSERVER" ) ;
288 sdh.observer = tcolr->asString( obsId ) ;
289 tpoolr->destroy( tcolr ) ;
290 }
291 if ( sdh.project == "" ) {
292 tcolr = tpoolr->construct( obstab, "PROJECT" ) ;
293 sdh.observer = tcolr->asString( obsId ) ;
294 tpoolr->destroy( tcolr ) ;
295 }
296 ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ;
297 MEpoch me = (*tmpMeasCol)( obsId )( IPosition(1,0) ) ;
298 delete tmpMeasCol ;
299 if ( sdh.utc == 0.0 ) {
300 sdh.utc = me.get( "s" ).getValue() ;
301 }
302 if ( telescopeName == "" ) {
303 tcolr = tpoolr->construct( obstab, "TELESCOPE_NAME" ) ;
304 sdh.observer = tcolr->asString( obsId ) ;
305 tpoolr->destroy( tcolr ) ;
306 }
307 Int nbeam = 0 ;
308 //
309 // ITERATION: FEED1
310 //
311 Int added1 = 0 ;
312 Int current1 = table_->nrow() ;
313 TableIterator iter1( t0, "FEED1" ) ;
314 while( !iter1.pastEnd() ) {
315 Table t1 = iter1.table() ;
316 // assume FEED1 == FEED2
317 tcolr = tpoolr->construct( t1, "FEED1" ) ;
318 Int feedId = tcolr->asInt( 0 ) ;
319 tpoolr->destroy( tcolr ) ;
320 nbeam++ ;
321 //
322 // ITERATION: FIELD_ID
323 //
324 Int added2 = 0 ;
325 Int current2 = table_->nrow() ;
326 TableIterator iter2( t1, "FIELD_ID" ) ;
327 while( !iter2.pastEnd() ) {
328 Table t2 = iter2.table() ;
329 tcolr = tpoolr->construct( t2, "FIELD_ID" ) ;
330 Int fieldId = tcolr->asInt( 0 ) ;
331 tpoolr->destroy( tcolr ) ;
332 tcolr = tpoolr->construct( fieldtab, "SOURCE_ID" ) ;
333 Int srcId = tcolr->asInt( fieldId ) ;
334 tpoolr->destroy( tcolr ) ;
335 tcolr = tpoolr->construct( fieldtab, "NAME" ) ;
336 String fieldName = tcolr->asString( fieldId ) + "__" + String::toString(fieldId) ;
337 tpoolr->destroy( tcolr ) ;
338 //
339 // ITERATION: DATA_DESC_ID
340 //
341 Int added3 = 0 ;
342 Int current3 = table_->nrow() ;
343 TableIterator iter3( t2, "DATA_DESC_ID" ) ;
344 while( !iter3.pastEnd() ) {
345 Table t3 = iter3.table() ;
346 tcolr = tpoolr->construct( t3, "DATA_DESC_ID" ) ;
347 Int ddId = tcolr->asInt( 0 ) ;
348 tpoolr->destroy( tcolr ) ;
349 tcolr = tpoolr->construct( ddtab, "POLARIZATION_ID" ) ;
350 Int polId = tcolr->asInt( ddId ) ;
351 tpoolr->destroy( tcolr ) ;
352 tcolr = tpoolr->construct( ddtab, "SPECTRAL_WINDOW_ID" ) ;
353 Int spwId = tcolr->asInt( ddId ) ;
354 tpoolr->destroy( tcolr ) ;
355 // polarization information
356 tcolr = tpoolr->construct( poltab, "NUM_CORR" ) ;
357 Int npol = tcolr->asInt( polId ) ;
358 tpoolr->destroy( tcolr ) ;
359 ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ;
360 Vector<Int> corrtype = (*roArrICol)( polId ) ;
361 delete roArrICol ;
362 //os_ << "npol = " << npol << LogIO::POST ;
363 //os_ << "corrtype = " << corrtype << LogIO::POST ;
364 if ( sdh.npol < npol ) sdh.npol = npol ;
365 if ( sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ;
366 // source information
367 MSSource srctabSel( srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ) ;
368 //os_ << "srcId = " << srcId << " spwId = " << spwId << " nrow = " << srctabSel.nrow() << LogIO::POST ;
369 tcolr = tpoolr->construct( srctabSel, "NAME" ) ;
370 String srcName = tcolr->asString( 0 ) ;
371 tpoolr->destroy( tcolr ) ;
372 //os_ << "srcName = " << srcName << LogIO::POST ;
373 ROArrayColumn<Double> *roArrDCol = new ROArrayColumn<Double>( srctabSel, "PROPER_MOTION" ) ;
374 Array<Double> srcPM = (*roArrDCol)( 0 ) ;
375 delete roArrDCol ;
376 //os_ << "srcPM = " << srcPM << LogIO::POST ;
377 roArrDCol = new ROArrayColumn<Double>( srctabSel, "DIRECTION" ) ;
378 Array<Double> srcDir = (*roArrDCol)( 0 ) ;
379 delete roArrDCol ;
380 //os_ << "srcDir = " << srcDir << LogIO::POST ;
381 Array<Double> sysVels ;
382 Double sysVel = 0.0 ;
383 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) {
384 roArrDCol = new ROArrayColumn<Double>( srctabSel, "SYSVEL" ) ;
385 sysVels = (*roArrDCol)( 0 ) ;
386 delete roArrDCol ;
387 }
388 if ( !sysVels.empty() ) {
389 //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ;
390 // NB: assume all SYSVEL values are the same
391 sysVel = sysVels( IPosition(1,0) ) ;
392 }
393 //delete tmpArrCol ;
394 //os_ << "sysVel = " << sysVel << LogIO::POST ;
395 ROScalarMeasColumn<MDirection> *tmpMeasCol = new ROScalarMeasColumn<MDirection>( srctabSel, "DIRECTION" ) ;
396 MDirection md = (*tmpMeasCol)( 0 ) ;
397 delete tmpMeasCol ;
398 // for MOLECULES subtable
399 tcolr = tpoolr->construct( srctabSel, "NUM_LINES" ) ;
400 Int numLines = tcolr->asInt( 0 ) ;
401 tpoolr->destroy( tcolr ) ;
402 //os_ << "numLines = " << numLines << LogIO::POST ;
403 Vector<Double> restFreqs( numLines, 0.0 ) ;
404 Vector<String> transitionName( numLines, "" ) ;
405 if ( numLines != 0 ) {
406 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) {
407 sharedQDArrCol = new ROArrayQuantColumn<Double>( srctabSel, "REST_FREQUENCY" ) ;
408 Array< Quantum<Double> > qRestFreqs = (*sharedQDArrCol)( 0 ) ;
409 delete sharedQDArrCol ;
410 for ( int i = 0 ; i < numLines ; i++ ) {
411 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ;
412 }
413 }
414 //os_ << "restFreqs = " << restFreqs << LogIO::POST ;
415 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) {
416 tcolr = tpoolr->construct( srctabSel, "TRANSITION" ) ;
417 transitionName = tcolr->asString( 0 ) ;
418 tpoolr->destroy( tcolr ) ;
419 //os_ << "transitionNameCol.nrow() = " << transitionNameCol.nrow() << LogIO::POST ;
420 }
421 }
422 uInt molId = table_->molecules().addEntry( restFreqs, transitionName, transitionName ) ;
423 // spectral setup
424 MeasFrame mf( me, mp, md ) ;
425 tcolr = tpoolr->construct( spwtab, "MEAS_FREQ_REF" ) ;
426 MFrequency::Types freqRef = MFrequency::castType((uInt)(tcolr->asInt(spwId))) ;
427 tpoolr->destroy( tcolr ) ;
428 tcolr = tpoolr->construct( spwtab, "NUM_CHAN" ) ;
429 Int nchan = tcolr->asInt( spwId ) ;
430 tpoolr->destroy( tcolr ) ;
431 Bool even = False ;
432 if ( (nchan/2)*2 == nchan ) even = True ;
433 if ( sdh.nchan < nchan ) sdh.nchan = nchan ;
434 ROScalarQuantColumn<Double> *tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "TOTAL_BANDWIDTH" ) ;
435 Double totbw = (*tmpQuantCol)( spwId ).getValue( "Hz" ) ;
436 delete tmpQuantCol ;
437 if ( sdh.bandwidth < totbw ) sdh.bandwidth = totbw ;
438 if ( sdh.freqref == "" )
439 //sdh.freqref = MFrequency::showType( freqRef ) ;
440 sdh.freqref = "LSRK" ;
441 if ( sdh.reffreq == -1.0 ) {
442 tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "REF_FREQUENCY" ) ;
443 Quantum<Double> qreffreq = (*tmpQuantCol)( spwId ) ;
444 delete tmpQuantCol ;
445 if ( freqRef == MFrequency::LSRK ) {
446 sdh.reffreq = qreffreq.getValue("Hz") ;
447 }
448 else {
449 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
450 sdh.reffreq = tolsr( qreffreq ).get("Hz").getValue() ;
451 }
452 }
453 Int refchan = nchan / 2 ;
454 IPosition refip( 1, refchan ) ;
455 Double refpix = 0.5*(nchan-1) ;
456 Double refval = 0.0 ;
457 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_WIDTH" ) ;
458 Double increment = (*sharedQDArrCol)( spwId )( refip ).getValue( "Hz" ) ;
459 delete sharedQDArrCol ;
460 //os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ;
461 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_FREQ" ) ;
462 Vector< Quantum<Double> > chanFreqs = (*sharedQDArrCol)( spwId ) ;
463 delete sharedQDArrCol ;
464 if ( freqRef == MFrequency::LSRK ) {
465 if ( even ) {
466 IPosition refip0( 1, refchan-1 ) ;
467 Double refval0 = chanFreqs(refip0).getValue("Hz") ;
468 Double refval1 = chanFreqs(refip).getValue("Hz") ;
469 refval = 0.5 * ( refval0 + refval1 ) ;
470 }
471 else {
472 refval = chanFreqs(refip).getValue("Hz") ;
473 }
474 }
475 else {
476 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
477 if ( even ) {
478 IPosition refip0( 1, refchan-1 ) ;
479 Double refval0 = chanFreqs(refip0).getValue("Hz") ;
480 Double refval1 = chanFreqs(refip).getValue("Hz") ;
481 refval = 0.5 * ( refval0 + refval1 ) ;
482 refval = tolsr( refval ).get("Hz").getValue() ;
483 }
484 else {
485 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ;
486 }
487 }
488 uInt freqId = table_->frequencies().addEntry( refpix, refval, increment ) ;
489 if ( ifmap.find( spwId ) == ifmap.end() ) {
490 ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ;
491 //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ;
492 }
493 // for TSYS and TCAL
494 MSSysCal caltabsel( caltab( caltab.col("ANTENNA_ID") == antenna_ && caltab.col("FEED_ID") == feedId && caltab.col("SPECTRAL_WINDOW_ID") == spwId ).sort("TIME") ) ;
495 //
496 // ITERATION: SCAN_NUMBER
497 //
498 Int added4 = 0 ;
499 Int current4 = table_->nrow() ;
500 TableIterator iter4( t3, "SCAN_NUMBER" ) ;
501 while( !iter4.pastEnd() ) {
502 Table t4 = iter4.table() ;
503 tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ;
504 Int scanNum = tcolr->asInt( 0 ) ;
505 tpoolr->destroy( tcolr ) ;
506 uInt cycle = 0 ;
507 //
508 // ITERATION: STATE_ID
509 //
510 Int added5 = 0 ;
511 Int current5 = table_->nrow() ;
512 TableIterator iter5( t4, "STATE_ID" ) ;
513 while( !iter5.pastEnd() ) {
514 Table t5 = iter5.table() ;
515 tcolr = tpoolr->construct( t5, "STATE_ID" ) ;
516 Int stateId = tcolr->asInt( 0 ) ;
517 tpoolr->destroy( tcolr ) ;
518 tcolr = tpoolr->construct( stattab, "OBS_MODE" ) ;
519 String obstype = tcolr->asString( stateId ) ;
520 tpoolr->destroy( tcolr ) ;
521 if ( sdh.obstype == "" ) sdh.obstype = obstype ;
522
523 Int nrow = t5.nrow() ;
524 Int prevnr = oldnr ;
525 Int addednr = 0 ;
526 Int nloop = 0 ;
527
528 // SPECTRA, FLAG
529 ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ;
530 if ( isFloatData_ ) {
531 //os_ << "FLOAT_DATA exists" << LogIO::POST ;
532 ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ;
533 addednr = nrow*npol ;
534 oldnr += addednr ;
535 table_->table().addRow( addednr ) ;
536 nloop = npol ;
537 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
538 Matrix<Float> sp = mFloatDataCol( irow ) ;
539 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
540 spCol->put( prevnr+ipol*nrow+irow, sp.row(ipol) ) ;
541 }
542 }
543 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
544 Matrix<Bool> flb = mFlagCol( irow ) ;
545 Matrix<uChar> fl( flb.shape() ) ;
546 convertArray( fl, flb ) ;
547 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
548 flCol->put( prevnr+ipol*nrow+irow, fl.row(ipol) ) ;
549 }
550 }
551 if ( sdh.fluxunit == "" ) {
552 const TableRecord dataColKeys = mFloatDataCol.keywordSet() ;
553 if ( dataColKeys.isDefined( "UNIT" ) )
554 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
555 }
556 }
557 else if ( isData_ ) {
558 //os_ << "DATA exists" << LogIO::POST ;
559 ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ;
560 addednr = nrow*npol ;
561 oldnr += addednr ;
562 table_->table().addRow( addednr ) ;
563 nloop = npol ;
564 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
565 Bool crossOK = False ;
566 Matrix<Complex> sp = mDataCol( irow ) ;
567 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
568 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
569 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
570 if ( !crossOK ) {
571 crossOK = True ;
572 Int pidx = prevnr + ipol * nrow + irow ;
573 spCol->put( pidx, real( sp.row(ipol) ) ) ;
574 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL )
575 spCol->put( pidx+nrow, imag( sp.row(ipol) ) ) ;
576 else
577 spCol->put( pidx+nrow, imag( conj(sp.row(ipol)) ) ) ;
578 }
579 }
580 else {
581 spCol->put( prevnr+ipol*nrow+irow, real( sp.row(ipol) ) ) ;
582 }
583 }
584 }
585 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
586 Bool crossOK = False ;
587 Matrix<Bool> flb = mFlagCol( irow ) ;
588 Matrix<uChar> fl( flb.shape() ) ;
589 convertArray( fl, flb ) ;
590 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
591 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX
592 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) {
593 if ( !crossOK ) {
594 crossOK = True ;
595 Int pidx = prevnr + ipol * nrow + irow ;
596 flCol->put( pidx, fl.row(ipol) ) ;
597 flCol->put( pidx+nrow, fl.row(ipol+1) ) ;
598 }
599 }
600 else {
601 flCol->put( prevnr+ipol*nrow+irow, fl.row(ipol) ) ;
602 }
603 }
604 }
605 if ( sdh.fluxunit == "" ) {
606 const TableRecord dataColKeys = mDataCol.keywordSet() ;
607 if ( dataColKeys.isDefined( "UNIT" ) )
608 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ;
609 }
610 }
611
612 // number of rows added in this cycle
613 //os_ << "prevnr = " << prevnr << LogIO::POST ;
614 //os_ << "addednr = " << addednr << LogIO::POST ;
615 RefRows rows( prevnr, prevnr+addednr-1 ) ;
616
617 // TIME
618 ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ;
619 Int tidx = prevnr ;
620 for ( Int i = 0 ; i < nloop ; i++ ) {
621 for ( Int j = 0 ; j < nrow ; j++ ) {
622 timeCol->put( tidx++, (*mTimeCol)( j ) ) ;
623 }
624 }
625
626 // TSYS
627 Vector<Double> sysCalTime ;
628 if ( isSysCal_ ) {
629 sysCalTime = getSysCalTime( caltabsel, *mTimeCol ) ;
630 tidx = prevnr ;
631 uInt calidx = 0 ;
632 for ( Int i = 0 ; i < nrow ; i++ ) {
633 Matrix<Float> tsys ;
634 calidx = getTsys( calidx, tsys, caltabsel, sysCalTime(i) ) ;
635 //os_ << "tsys = " << tsys << LogIO::POST ;
636 uInt ncol = tsys.ncolumn() ;
637 if ( ncol == 0 ) {
638 IPosition cShape = IPosition( 2, npol, 1 ) ;
639 tsys.resize( cShape ) ;
640 tsys = 1.0 ;
641 }
642 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
643 //floatArrCol->put( prevnr+i+nrow*ipol, tsys.row( ipol ) ) ;
644 tsysCol->put( prevnr+i+nrow*ipol, tsys.row( ipol ) ) ;
645 }
646 }
647 }
648 else {
649 Vector<Float> tsys( 1, 1.0 ) ;
650 for ( Int i = prevnr ; i < prevnr+addednr ; i++ )
651 tsysCol->put( i, tsys ) ;
652 }
653
654
655 // INTERVAL
656 tcolr = tpoolr->construct( t5, "INTERVAL" ) ;
657 //Vector<Double> integ = mIntervalCol->getColumn() ;
658 for ( int i = 0 ; i < nloop ; i++ ) {
659 //Int startrow = prevnr + i ;
660 //Int endrow = startrow + nrow - 1 ;
661 //RefRows prows( startrow, endrow ) ;
662 //intervalCol->putColumnCells( prows, integ ) ;
663 for ( int j = 0 ; j < nrow ; j++ ) {
664 intervalCol->putScalar( prevnr+i*nrow+j, (Double)(tcolr->asdouble( j )) ) ;
665 }
666 }
667 tpoolr->destroy( tcolr ) ;
668
669 // SRCTYPE
670 Int srcType = getSrcType( stateId, tpoolr ) ;
671 for ( int i = 0 ; i < addednr ; i++ ) {
672 srctypeCol->putScalar( prevnr+i, srcType ) ;
673 }
674 //Vector<Int> *srcType = new Vector<Int>( addednr, getSrcType( stateId ) ) ;
675 //srcTypeCol->putColumnCells( rows, *srcType ) ;
676
677 // DIRECTION, AZIMUTH, ELEVATION, SCANRATE
678 Vector<Double> defaultScanrate( 2, 0.0 ) ;
679 uInt diridx = 0 ;
680 MDirection::Types dirType ;
681 if ( getPt_ ) {
682 for ( Int i = 0 ; i < nrow ; i++ ) {
683 Vector<Double> dir ;
684 Vector<Double> scanrate ;
685 String refString ;
686 diridx = getDirection( diridx, dir, scanrate, refString, pointtab, (*mTimeCol)(i).get("s").getValue() ) ;
687 //os_ << "diridx = " << diridx << " dmTimeCol(" << i << ") = " << mTimeCol(i).get("s").getValue()-mTimeCol(0).get("s").getValue() << LogIO::POST ;
688 //os_ << "dir = " << dir << LogIO::POST ;
689 //os_ << "scanrate = " << scanrate << LogIO::POST ;
690 //os_ << "refString = " << refString << LogIO::POST ;
691 MDirection::getType( dirType, refString ) ;
692 //os_ << "dirType = " << dirType << LogIO::POST ;
693 mf.resetEpoch( (*mTimeCol)(i) ) ;
694 mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ;
695 if ( refString == "J2000" ) {
696 //os_ << "J2000" << LogIO::POST ;
697 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
698 dirCol->put( prevnr+i+nrow*ipol, dir ) ;
699 }
700 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
701 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
702 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
703 azCol->putScalar( prevnr+i+nrow*ipol, (Float)azel(0) ) ;
704 elCol->putScalar( prevnr+i+nrow*ipol, (Float)azel(1) ) ;
705 }
706 }
707 else if ( refString(0,4) == "AZEL" ) {
708 //os_ << "AZEL" << LogIO::POST ;
709 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
710 azCol->putScalar( prevnr+i+nrow*ipol, (Float)dir(0) ) ;
711 elCol->putScalar( prevnr+i+nrow*ipol, (Float)dir(1) ) ;
712 }
713 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
714 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
715 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
716 dirCol->put( prevnr+i+nrow*ipol, newdir ) ;
717 }
718 }
719 else {
720 //os_ << "OTHER: " << refString << LogIO::POST ;
721 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
722 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ;
723 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
724 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ;
725 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
726 dirCol->put( prevnr+i+nrow*ipol, newdir ) ;
727 azCol->putScalar( prevnr+i+nrow*ipol, (Float)dir(0) ) ;
728 elCol->putScalar( prevnr+i+nrow*ipol, (Float)dir(1) ) ;
729 }
730 }
731 if ( scanrate.size() != 0 ) {
732 //os_ << "scanrate.size() = " << scanrate.size() << LogIO::POST ;
733 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
734 scanrateCol->put( prevnr+i+nrow*ipol, scanrate ) ;
735 }
736 }
737 else {
738 //os_ << "scanrate.size() = " << scanrate.size() << LogIO::POST ;
739 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
740 scanrateCol->put( prevnr+i+nrow*ipol, defaultScanrate ) ;
741 }
742 }
743 }
744 }
745 else {
746 // All directions are set to source direction
747 ROArrayMeasColumn<MDirection> dmcol( pointtab, "DIRECTION" ) ;
748 ROArrayColumn<Double> dcol( pointtab, "DIRECTION" ) ;
749 IPosition ip( dmcol(0).shape().nelements(), 0 ) ;
750 IPosition outp( 1, 2 ) ;
751 String ref = dmcol(0)(ip).getRefString() ;
752 Slice ds( 0, 2, 1 ) ;
753 Slice ds0( 0, 1, 1 ) ;
754 Slicer dslice0( ds, ds0 ) ;
755 Vector<Double> defaultDir = dcol(0)(dslice0).reform(outp) ;
756 MDirection::getType( dirType, "J2000" ) ;
757 mf.resetDirection( MDirection( MVDirection(srcDir), dirType ) ) ;
758 if ( ref != "J2000" ) {
759 ROScalarMeasColumn<MEpoch> tmCol( pointtab, "TIME" ) ;
760 mf.resetEpoch( tmCol( 0 ) ) ;
761 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
762 defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ;
763 }
764 for ( Int i = 0 ; i < nrow ; i++ ) {
765 mf.resetEpoch( (*mTimeCol)(i) ) ;
766 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
767 Int localidx = prevnr+i+nrow*ipol ;
768 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
769 Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ;
770 azCol->putScalar( localidx, (Float)azel(0) ) ;
771 elCol->putScalar( localidx, (Float)azel(1) ) ;
772 dirCol->put( localidx, defaultDir ) ;
773 scanrateCol->put( localidx, defaultScanrate ) ;
774 }
775 }
776 }
777
778 // CYCLENO
779 for ( int i = 0 ; i < nloop ; i++ ) {
780 for ( int j = 0 ; j < nrow ; j++ ) {
781 cyclenoCol->putScalar( prevnr+nrow*i+j, cycle+j ) ;
782 }
783 }
784 cycle += nrow ;
785
786 // POLNO
787 Int pidx = 0 ;
788 Bool crossOK = False ;
789 for ( int i = 0 ; i < npol ; i++ ) {
790 Vector<uInt> polnos = getPolNo( corrtype[i] ) ;
791 if ( polnos.size() > 1 ) {
792 if ( crossOK ) continue ;
793 else crossOK = True ;
794 }
795 for ( uInt j = 0 ; j < polnos.size() ; j++ ) {
796 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
797 polnoCol->putScalar( prevnr+pidx*nrow+irow, polnos[j] ) ;
798 }
799 pidx++ ;
800 }
801 }
802
803 // FLAGROW
804 tcolr = tpoolr->construct( t5, "FLAG_ROW" ) ;
805 for ( int i = 0 ; i < nloop ; i++ ) {
806 for ( int j = 0 ; j < nrow ; j++ ) {
807 flagrowCol->putScalar( prevnr+nrow*i+j, (uInt)(tcolr->asBool( j )) ) ;
808 }
809 }
810 tpoolr->destroy( tcolr ) ;
811
812 // TCAL_ID
813 if ( isSysCal_ ) {
814 for( Int irow = 0 ; irow < nrow ; irow++ ) {
815 Vector<uInt> tcalids = getTcalId( feedId, spwId, sysCalTime[irow] ) ;
816 if ( tcalids.size() == 0 ) {
817 tcalids.resize( npol ) ;
818 tcalids = 0 ;
819 }
820 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
821 tcalidCol->putScalar( prevnr+irow+nrow*ipol, tcalids[ipol] ) ;
822 }
823 }
824 }
825 else {
826 //Vector<uInt> tcalid( addednr, 0 ) ;
827 //uIntCol->putColumnCells( rows, tcalid ) ;
828 uInt tcalid = 0 ;
829 for ( Int irow = 0 ; irow < addednr ; irow++ )
830 tcalidCol->putScalar( prevnr+irow, tcalid ) ;
831 }
832
833 // WEATHER_ID
834 if ( isWeather_ ) {
835 uInt wid = 0 ;
836 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
837 wid = getWeatherId( wid, (*mTimeCol)(irow).get("s").getValue() ) ;
838 for ( Int ipol = 0 ; ipol < nloop ; ipol++ ) {
839 weatheridCol->putScalar( prevnr+ipol*nrow+irow, wid ) ;
840 }
841 }
842 }
843
844 delete mTimeCol ;
845
846 //os_ << "field: " << fieldId << " scan: " << scanNum << " obs: " << obsId << " state: " << stateId << " ddid: " << ddId << endl ;
847 os_ << "addednr = " << addednr << endl ;
848 added5 += addednr ;
849 iter5.next() ;
850 }
851
852 // SCANNO
853 // MS: 1-base
854 // Scantable: 0-base
855 Int dest5 = current5 + added5 ;
856 scanNum -= 1 ;
857 for ( Int irow = current5 ; irow < dest5 ; irow++ )
858 scannoCol->putScalar( irow, (uInt)scanNum ) ;
859
860 os_ << "added5 = " << added5 << endl ;
861 added4 += added5 ;
862 iter4.next() ;
863 }
864
865 // IFNO
866 Int dest4 = current4 + added4 ;
867 for ( Int irow = current4 ; irow < dest4 ; irow++ )
868 ifnoCol->putScalar( irow, (uInt)spwId ) ;
869
870 // FREQ_ID
871 uInt fid = ifmap[spwId] ;
872 for ( Int irow = current4 ; irow < dest4 ; irow++ )
873 freqidCol->putScalar( irow, fid ) ;
874
875 // MOLECULE_ID
876 for ( Int irow = current4 ; irow < dest4 ; irow++ )
877 molidCol->putScalar( irow, molId ) ;
878
879 // SRCNAME
880 for ( Int irow = current4 ; irow < dest4 ; irow++ )
881 srcnameCol->putScalar( irow, srcName ) ;
882
883 // SRCVELOCITY, SRCPROPERMOTION and SRCDIRECTION
884 // no reference conversion for direction at the moment (assume J2000)
885 // no reference conversion for velocity at the moment (assume LSRK)
886 for ( Int irow = current4 ; irow < dest4 ; irow++ )
887 srcpmCol->put( irow, srcPM ) ;
888 for ( Int irow = current4 ; irow < dest4 ; irow++ )
889 srcdirCol->put( irow, srcDir ) ;
890 for ( Int irow = current4 ; irow < dest4 ; irow++ )
891 srcvelCol->putScalar( irow, sysVel ) ;
892
893 os_ << "added4 = " << added4 << endl ;
894 added3 += added4 ;
895 iter3.next() ;
896 }
897
898 // FIELDNAME
899 Int dest3 = current3 + added3 ;
900 for ( Int irow = current3 ; irow < dest3 ; irow++ )
901 fieldnameCol->putScalar( irow, fieldName ) ;
902
903 os_ << "added3 = " << added3 << endl ;
904 added2 += added3 ;
905 iter2.next() ;
906 }
907
908 // BEAMNO
909 Int dest2 = current2 + added2 ;
910 for ( Int irow = current2 ; irow < dest2 ; irow++ )
911 beamnoCol->putScalar( irow, (uInt)feedId ) ;
912
913 // FOCUS_ID
914 // tentative
915 uInt focusId = 0 ;
916 for ( Int irow = current2 ; irow < dest2 ; irow++ )
917 focusidCol->putScalar( irow, focusId ) ;
918
919 os_ << "added2 = " << added2 << endl ;
920 added1 += added2 ;
921 iter1.next() ;
922 }
923 if ( sdh.nbeam < nbeam ) sdh.nbeam = nbeam ;
924
925 os_ << "added1 = " << added1 << endl ;
926 added0 += added1 ;
927 iter0.next() ;
928 }
929
930 os_ << "added0 = " << added0 << endl ;
931
932 // REFBEAMNO
933 // set 0 at the moment
934 tcolw = tpoolw->construct( table_->table(), "REFBEAMNO" ) ;
935 for ( Int irow = current0 ; irow < added0 ; irow++ )
936 tcolw->putScalar( irow, 0 ) ;
937 tpoolw->destroy( tcolw ) ;
938
939 // FIT_ID
940 // nothing to do
941 tcolw = tpoolw->construct( table_->table(), "FIT_ID" ) ;
942 for ( Int irow = current0 ; irow < added0 ; irow++ )
943 tcolw->putScalar( irow, -1 ) ;
944 tpoolw->destroy( tcolw ) ;
945
946 // OPACITY
947 // not used?
948 tcolw = tpoolw->construct( table_->table(), "OPACITY" ) ;
949 for ( Int irow = current0 ; irow < added0 ; irow++ )
950 tcolw->putScalar( irow, 0.0 ) ;
951 tpoolw->destroy( tcolw ) ;
952
953 // delete Scantable columns
954 tpoolw->destroy( scannoCol ) ;
955 tpoolw->destroy( cyclenoCol ) ;
956 tpoolw->destroy( beamnoCol ) ;
957 tpoolw->destroy( ifnoCol ) ;
958 tpoolw->destroy( polnoCol ) ;
959 tpoolw->destroy( freqidCol ) ;
960 tpoolw->destroy( molidCol ) ;
961 tpoolw->destroy( flagrowCol ) ;
962 delete timeCol ;
963 tpoolw->destroy( intervalCol ) ;
964 tpoolw->destroy( srcnameCol ) ;
965 tpoolw->destroy( srctypeCol ) ;
966 tpoolw->destroy( fieldnameCol ) ;
967 delete spCol ;
968 delete flCol ;
969 delete tsysCol ;
970 delete dirCol ;
971 tpoolw->destroy( azCol ) ;
972 tpoolw->destroy( elCol ) ;
973 tpoolw->destroy( tcalidCol ) ;
974 tpoolw->destroy( focusidCol ) ;
975 tpoolw->destroy( weatheridCol ) ;
976 delete srcpmCol ;
977 delete srcdirCol ;
978 tpoolw->destroy( srcvelCol ) ;
979 delete scanrateCol ;
980
981 delete tpoolr ;
982 delete tpoolw ;
983
984
985 // Table Keywords
986 sdh.nif = ifmap.size() ;
987 if ( ( telescopeName == "" ) || ( antennaName == telescopeName ) ) {
988 sdh.antennaname = antennaName ;
989 }
990 else {
991 sdh.antennaname = telescopeName + "//" + antennaName ;
992 }
993 if ( stationName != "" ) {
994 sdh.antennaname += "@" + stationName ;
995 }
996 ROArrayColumn<Double> pdirCol( pointtab, "DIRECTION" ) ;
997 String dirref = pdirCol.keywordSet().asRecord("MEASINFO").asString("Ref") ;
998 if ( dirref == "AZELGEO" || dirref == "AZEL" ) {
999 dirref = "J2000" ;
1000 }
1001 sscanf( dirref.chars()+1, "%f", &sdh.equinox ) ;
1002 sdh.epoch = "UTC" ;
1003 if (sdh.freqref == "TOPO") {
1004 sdh.freqref = "TOPOCENT";
1005 } else if (sdh.freqref == "GEO") {
1006 sdh.freqref = "GEOCENTR";
1007 } else if (sdh.freqref == "BARY") {
1008 sdh.freqref = "BARYCENT";
1009 } else if (sdh.freqref == "GALACTO") {
1010 sdh.freqref = "GALACTOC";
1011 } else if (sdh.freqref == "LGROUP") {
1012 sdh.freqref = "LOCALGRP";
1013 } else if (sdh.freqref == "CMB") {
1014 sdh.freqref = "CMBDIPOL";
1015 } else if (sdh.freqref == "REST") {
1016 sdh.freqref = "SOURCE";
1017 }
1018 table_->setHeader( sdh ) ;
1019
1020 // save path to POINTING table
1021 //Path datapath(mstable_.tableName()) ;
1022 Path datapath( tablename_ ) ;
1023 String pTabName = datapath.absoluteName() + "/POINTING" ;
1024 table_->table().rwKeywordSet().define( "POINTING", pTabName ) ;
1025
1026 // for GBT
1027 if ( antennaName == "GBT" ) {
1028 String goTabName = datapath.absoluteName() + "/GBT_GO" ;
1029 table_->table().rwKeywordSet().define( "GBT_GO", goTabName ) ;
1030 }
1031 //double endSec = gettimeofday_sec() ;
1032 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1033}
1034
1035void MSFiller::close()
1036{
1037 //tablesel_.closeSubTables() ;
1038 mstable_.closeSubTables() ;
1039 //tablesel_.unlock() ;
1040 mstable_.unlock() ;
1041}
1042
1043Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool )
1044{
1045 //double startSec = gettimeofday_sec() ;
1046 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;
1047
1048 MSState statetab = mstable_.state() ;
1049 ROTableColumn *sharedCol ;
1050 sharedCol = tpool->construct( statetab, "OBS_MODE" ) ;
1051 String obsMode = sharedCol->asString( stateId ) ;
1052 tpool->destroy( sharedCol ) ;
1053 sharedCol = tpool->construct( statetab, "SIG" ) ;
1054 Bool sig = sharedCol->asBool( stateId ) ;
1055 tpool->destroy( sharedCol ) ;
1056 sharedCol = tpool->construct( statetab, "REF" ) ;
1057 Bool ref = sharedCol->asBool( stateId ) ;
1058 tpool->destroy( sharedCol ) ;
1059 sharedCol = tpool->construct( statetab, "CAL" ) ;
1060 Double cal = (Double)(sharedCol->asdouble( stateId )) ;
1061 tpool->destroy( sharedCol ) ;
1062 //os_ << "OBS_MODE = " << obsMode << LogIO::POST ;
1063
1064 // determine separator
1065 String sep = "" ;
1066 if ( obsMode.find( ":" ) != String::npos ) {
1067 sep = ":" ;
1068 }
1069 else if ( obsMode.find( "." ) != String::npos ) {
1070 sep = "." ;
1071 }
1072
1073 // determine SRCTYPE
1074 Int srcType = SrcType::NOTYPE ;
1075 if ( sep == ":" ) {
1076 // sep == ":"
1077 //
1078 // GBT case
1079 //
1080 // obsMode1=Nod
1081 // NOD
1082 // obsMode1=OffOn
1083 // obsMode2=PSWITCHON: PSON
1084 // obsMode2=PSWITCHOFF: PSOFF
1085 // obsMode1=??
1086 // obsMode2=FSWITCH:
1087 // SIG=1: FSON
1088 // REF=1: FSOFF
1089 // Calibration scan if CAL != 0
1090 Int epos = obsMode.find_first_of( sep ) ;
1091 Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
1092 String obsMode1 = obsMode.substr( 0, epos ) ;
1093 String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
1094 if ( obsMode1 == "Nod" ) {
1095 srcType = SrcType::NOD ;
1096 }
1097 else if ( obsMode1 == "OffOn" ) {
1098 if ( obsMode2 == "PSWITCHON" ) srcType = SrcType::PSON ;
1099 if ( obsMode2 == "PSWITCHOFF" ) srcType = SrcType::PSOFF ;
1100 }
1101 else {
1102 if ( obsMode2 == "FSWITCH" ) {
1103 if ( sig ) srcType = SrcType::FSON ;
1104 if ( ref ) srcType = SrcType::FSOFF ;
1105 }
1106 }
1107 if ( cal > 0.0 ) {
1108 if ( srcType == SrcType::NOD )
1109 srcType = SrcType::NODCAL ;
1110 else if ( srcType == SrcType::PSON )
1111 srcType = SrcType::PONCAL ;
1112 else if ( srcType == SrcType::PSOFF )
1113 srcType = SrcType::POFFCAL ;
1114 else if ( srcType == SrcType::FSON )
1115 srcType = SrcType::FONCAL ;
1116 else if ( srcType == SrcType::FSOFF )
1117 srcType = SrcType::FOFFCAL ;
1118 else
1119 srcType = SrcType::CAL ;
1120 }
1121 }
1122 else if ( sep == "." ) {
1123 // sep == "."
1124 //
1125 // ALMA & EVLA case (MS via ASDM)
1126 //
1127 // obsMode1=CALIBRATE_*
1128 // obsMode2=ON_SOURCE: PONCAL
1129 // obsMode2=OFF_SOURCE: POFFCAL
1130 // obsMode1=OBSERVE_TARGET
1131 // obsMode2=ON_SOURCE: PON
1132 // obsMode2=OFF_SOURCE: POFF
1133 Int epos = obsMode.find_first_of( sep ) ;
1134 Int nextpos = obsMode.find_first_of( sep, epos+1 ) ;
1135 String obsMode1 = obsMode.substr( 0, epos ) ;
1136 String obsMode2 = obsMode.substr( epos+1, nextpos-epos-1 ) ;
1137 if ( obsMode1.find( "CALIBRATE_" ) == 0 ) {
1138 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PONCAL ;
1139 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::POFFCAL ;
1140 }
1141 else if ( obsMode1 == "OBSERVE_TARGET" ) {
1142 if ( obsMode2 == "ON_SOURCE" ) srcType = SrcType::PSON ;
1143 if ( obsMode2 == "OFF_SOURCE" ) srcType = SrcType::PSOFF ;
1144 }
1145 }
1146 else {
1147 if ( sig ) srcType = SrcType::SIG ;
1148 if ( ref ) srcType = SrcType::REF ;
1149 }
1150
1151 //os_ << "srcType = " << srcType << LogIO::POST ;
1152 //double endSec = gettimeofday_sec() ;
1153 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1154 return srcType ;
1155}
1156
1157Vector<uInt> MSFiller::getPolNo( Int corrType )
1158{
1159 //double startSec = gettimeofday_sec() ;
1160 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;
1161 Vector<uInt> polno( 1 ) ;
1162
1163 if ( corrType == Stokes::I || corrType == Stokes::RR || corrType == Stokes::XX ) {
1164 polno = 0 ;
1165 }
1166 else if ( corrType == Stokes::Q || corrType == Stokes::LL || corrType == Stokes::YY ) {
1167 polno = 1 ;
1168 }
1169 else if ( corrType == Stokes::U ) {
1170 polno = 2 ;
1171 }
1172 else if ( corrType == Stokes::V ) {
1173 polno = 3 ;
1174 }
1175 else if ( corrType == Stokes::RL || corrType == Stokes::XY || corrType == Stokes::LR || corrType == Stokes::RL ) {
1176 polno.resize( 2 ) ;
1177 polno[0] = 2 ;
1178 polno[1] = 3 ;
1179 }
1180 else if ( corrType == Stokes::Plinear ) {
1181 polno[0] = 1 ;
1182 }
1183 else if ( corrType == Stokes::Pangle ) {
1184 polno[0] = 2 ;
1185 }
1186 else {
1187 polno = 99 ;
1188 }
1189 //os_ << "polno = " << polno << LogIO::POST ;
1190 //double endSec = gettimeofday_sec() ;
1191 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1192
1193 return polno ;
1194}
1195
1196String MSFiller::getPolType( Int corrType )
1197{
1198 //double startSec = gettimeofday_sec() ;
1199 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;
1200 String poltype = "" ;
1201
1202 if ( corrType == Stokes::I || corrType == Stokes::Q || corrType == Stokes::U || corrType == Stokes::V )
1203 poltype = "stokes" ;
1204 else if ( corrType == Stokes::XX || corrType == Stokes::YY || corrType == Stokes::XY || corrType == Stokes::YX )
1205 poltype = "linear" ;
1206 else if ( corrType == Stokes::RR || corrType == Stokes::LL || corrType == Stokes::RL || corrType == Stokes::LR )
1207 poltype = "circular" ;
1208 else if ( corrType == Stokes::Plinear || corrType == Stokes::Pangle )
1209 poltype = "linpol" ;
1210
1211 //double endSec = gettimeofday_sec() ;
1212 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1213 return poltype ;
1214}
1215
1216void MSFiller::fillWeather()
1217{
1218 //double startSec = gettimeofday_sec() ;
1219 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;
1220 Table mWeather = mstable_.weather() ;
1221 //Table mWeatherSel = mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ;
1222 Table mWeatherSel( mWeather( mWeather.col("ANTENNA_ID") == antenna_ ).sort("TIME") ) ;
1223 //os_ << "mWeatherSel.nrow() = " << mWeatherSel.nrow() << LogIO::POST ;
1224 if ( mWeatherSel.nrow() == 0 ) {
1225 os_ << "No rows with ANTENNA_ID = " << antenna_ << ", Try -1..." << LogIO::POST ;
1226 mWeatherSel = Table( MSWeather( mWeather( mWeather.col("ANTENNA_ID") == -1 ) ) ) ;
1227 if ( mWeatherSel.nrow() == 0 ) {
1228 os_ << "No rows in WEATHER table" << LogIO::POST ;
1229 }
1230 }
1231 uInt wnrow = mWeatherSel.nrow() ;
1232 //os_ << "wnrow = " << wnrow << LogIO::POST ;
1233
1234 if ( wnrow == 0 )
1235 return ;
1236
1237 Table wtab = table_->weather().table() ;
1238 wtab.addRow( wnrow ) ;
1239
1240 ScalarColumn<Float> *fCol ;
1241 ROScalarColumn<Float> *sharedFloatCol ;
1242 if ( mWeatherSel.tableDesc().isColumn( "TEMPERATURE" ) ) {
1243 fCol = new ScalarColumn<Float>( wtab, "TEMPERATURE" ) ;
1244 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "TEMPERATURE" ) ;
1245 fCol->putColumn( *sharedFloatCol ) ;
1246 delete sharedFloatCol ;
1247 delete fCol ;
1248 }
1249 if ( mWeatherSel.tableDesc().isColumn( "PRESSURE" ) ) {
1250 fCol = new ScalarColumn<Float>( wtab, "PRESSURE" ) ;
1251 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "PRESSURE" ) ;
1252 fCol->putColumn( *sharedFloatCol ) ;
1253 delete sharedFloatCol ;
1254 delete fCol ;
1255 }
1256 if ( mWeatherSel.tableDesc().isColumn( "REL_HUMIDITY" ) ) {
1257 fCol = new ScalarColumn<Float>( wtab, "HUMIDITY" ) ;
1258 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "REL_HUMIDITY" ) ;
1259 fCol->putColumn( *sharedFloatCol ) ;
1260 delete sharedFloatCol ;
1261 delete fCol ;
1262 }
1263 if ( mWeatherSel.tableDesc().isColumn( "WIND_SPEED" ) ) {
1264 fCol = new ScalarColumn<Float>( wtab, "WINDSPEED" ) ;
1265 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_SPEED" ) ;
1266 fCol->putColumn( *sharedFloatCol ) ;
1267 delete sharedFloatCol ;
1268 delete fCol ;
1269 }
1270 if ( mWeatherSel.tableDesc().isColumn( "WIND_DIRECTION" ) ) {
1271 fCol = new ScalarColumn<Float>( wtab, "WINDAZ" ) ;
1272 sharedFloatCol = new ROScalarColumn<Float>( mWeatherSel, "WIND_DIRECTION" ) ;
1273 fCol->putColumn( *sharedFloatCol ) ;
1274 delete sharedFloatCol ;
1275 delete fCol ;
1276 }
1277 ScalarColumn<uInt> idCol( wtab, "ID" ) ;
1278 for ( uInt irow = 0 ; irow < wnrow ; irow++ )
1279 idCol.put( irow, irow ) ;
1280
1281 ROScalarQuantColumn<Double> tqCol( mWeatherSel, "TIME" ) ;
1282 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ;
1283 String tUnit = tqCol.getUnits() ;
1284 mwTime_ = tCol.getColumn() ;
1285 if ( tUnit == "d" )
1286 mwTime_ *= 86400.0 ;
1287 tqCol.attach( mWeatherSel, "INTERVAL" ) ;
1288 tCol.attach( mWeatherSel, "INTERVAL" ) ;
1289 String iUnit = tqCol.getUnits() ;
1290 mwInterval_ = tCol.getColumn() ;
1291 if ( iUnit == "d" )
1292 mwInterval_ *= 86400.0 ;
1293 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ;
1294 //double endSec = gettimeofday_sec() ;
1295 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1296}
1297
1298void MSFiller::fillFocus()
1299{
1300 //double startSec = gettimeofday_sec() ;
1301 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;
1302 // tentative
1303 Table tab = table_->focus().table() ;
1304 tab.addRow( 1 ) ;
1305 ScalarColumn<uInt> idCol( tab, "ID" ) ;
1306 idCol.put( 0, 0 ) ;
1307 //double endSec = gettimeofday_sec() ;
1308 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1309}
1310
1311void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr,
1312 boost::object_pool<TableColumn> *tpoolw )
1313{
1314 //double startSec = gettimeofday_sec() ;
1315 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;
1316
1317 //MSSysCal sctab = mstable_.sysCal() ;
1318 Table sctab = mstable_.sysCal() ;
1319 if ( sctab.nrow() == 0 ) {
1320 os_ << "No SysCal rows" << LogIO::POST ;
1321 return ;
1322 }
1323 Bool isSp = sctab.tableDesc().isColumn( "TCAL_SPECTRUM" ) ;
1324 //Table sctabsel = sctab( sctab.col("ANTENNA_ID") == antenna_ ) ;
1325 Table sctabsel( sctab( sctab.col("ANTENNA_ID") == antenna_ ) ) ;
1326 if ( sctabsel.nrow() == 0 ) {
1327 os_ << "No SysCal rows" << LogIO::POST ;
1328 return ;
1329 }
1330 ROArrayColumn<Float> *tmpTcalCol = new ROArrayColumn<Float>( sctabsel, "TCAL" ) ;
1331 uInt npol = tmpTcalCol->shape( 0 )(0) ;
1332 delete tmpTcalCol ;
1333 //os_ << "fillTcal(): npol = " << npol << LogIO::POST ;
1334 Table tab = table_->tcal().table() ;
1335 TableColumn *idCol = tpoolw->construct( tab, "ID" ) ;
1336 TableColumn *timeCol = tpoolw->construct( tab, "TIME" ) ;
1337 ArrayColumn<Float> tcalCol( tab, "TCAL" ) ;
1338 ROTableColumn *sharedCol ;
1339 ROArrayColumn<Float> scTcalCol ;
1340 uInt oldnr = 0 ;
1341 uInt newnr = 0 ;
1342 TableIterator iter0( sctabsel, "FEED_ID" ) ;
1343 while( !iter0.pastEnd() ) {
1344 Table t0 = iter0.table() ;
1345 sharedCol = tpoolr->construct( t0, "FEED_ID" ) ;
1346 Int feedId = sharedCol->asInt( 0 ) ;
1347 tpoolr->destroy( sharedCol ) ;
1348 String ffield = "FEED" + String::toString( feedId ) ;
1349 //Record rec0 ;
1350 TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ;
1351 while( !iter1.pastEnd() ) {
1352 Table t1 = iter1.table() ;
1353 sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ;
1354 Int spwId = sharedCol->asInt( 0 ) ;
1355 tpoolr->destroy( sharedCol ) ;
1356 String spwfield = "SPW" + String::toString( spwId ) ;
1357 //Record rec1 ;
1358 TableIterator iter2( t1, "TIME" ) ;
1359 while( !iter2.pastEnd() ) {
1360 Table t2 = iter2.table() ;
1361 uInt nrow = t2.nrow() ; // must be 1
1362 //os_ << "fillTcal::t2.nrow = " << nrow << LogIO::POST ;
1363 ROScalarQuantColumn<Double> scTimeCol( t2, "TIME" ) ;
1364 IPosition newShape( 2, 1, nrow ) ;
1365 if ( isSp ) {
1366 scTcalCol.attach( t2, "TCAL_SPECTRUM" ) ;
1367 newShape[0] = scTcalCol.shape(0)(1) ;
1368 }
1369 else {
1370 scTcalCol.attach( t1, "TCAL" ) ;
1371 }
1372 tab.addRow( nrow*npol ) ;
1373 newnr += nrow*npol ;
1374 String sTime = MVTime( scTimeCol(0) ).string( MVTime::YMD ) ;
1375 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1376 timeCol->putScalar( oldnr+ipol, sTime ) ;
1377 }
1378 uInt idx = oldnr ;
1379 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1380 idCol->putScalar( oldnr+ipol, idx++ ) ;
1381 }
1382 Vector<uInt> idminmax( 2, oldnr ) ;
1383 Matrix<Float> subtcal = scTcalCol( 0 ) ;
1384 for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
1385 tcalCol.put( oldnr+ipol, subtcal.row( ipol ) ) ;
1386 }
1387 idminmax[1] = newnr - 1 ;
1388 oldnr = newnr ;
1389
1390 String key = ffield+":"+spwfield+":"+sTime ;
1391 tcalrec_.define( key, idminmax ) ;
1392 //rec1.define( sTime, idminmax ) ;
1393 iter2++ ;
1394 }
1395 //rec0.defineRecord( spwfield, rec1 ) ;
1396 iter1++ ;
1397 }
1398 //tcalrec_.defineRecord( ffield, rec0 ) ;
1399 iter0++ ;
1400 }
1401
1402 tpoolw->destroy( idCol ) ;
1403 tpoolw->destroy( timeCol ) ;
1404
1405 //tcalrec_.print( std::cout ) ;
1406 //double endSec = gettimeofday_sec() ;
1407 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1408}
1409
1410uInt MSFiller::getWeatherId( uInt idx, Double wtime )
1411{
1412 //double startSec = gettimeofday_sec() ;
1413 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;
1414 uInt nrow = mwTime_.size() ;
1415 if ( nrow == 0 )
1416 return 0 ;
1417 uInt wid = nrow ;
1418 for ( uInt i = idx ; i < nrow-1 ; i++ ) {
1419 Double tStart = mwTime_[i]-0.5*mwInterval_[i] ;
1420 // use of INTERVAL column is problematic
1421 // since there are "blank" time of weather monitoring
1422 //Double tEnd = tStart + mwInterval_[i] ;
1423 Double tEnd = mwTime_[i+1]-0.5*mwInterval_[i+1] ;
1424 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1425 if ( wtime >= tStart && wtime <= tEnd ) {
1426 wid = i ;
1427 break ;
1428 }
1429 }
1430 if ( wid == nrow ) {
1431 uInt i = nrow - 1 ;
1432 Double tStart = mwTime_[i-1]+0.5*mwInterval_[i-1] ;
1433 Double tEnd = mwTime_[i]+0.5*mwInterval_[i] ;
1434 //os_ << "tStart = " << tStart << " dtEnd = " << tEnd-tStart << " dwtime = " << wtime-tStart << LogIO::POST ;
1435 if ( wtime >= tStart && wtime <= tEnd )
1436 wid = i ;
1437 }
1438
1439 //if ( wid == nrow )
1440 //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ;
1441
1442 //double endSec = gettimeofday_sec() ;
1443 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1444 return wid ;
1445}
1446
1447Vector<Double> MSFiller::getSysCalTime( MSSysCal &tab, MEpoch::ROScalarColumn &tcol )
1448{
1449 //double startSec = gettimeofday_sec() ;
1450 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;
1451 uInt nrow = tcol.table().nrow() ;
1452 Vector<Double> tstr( nrow, -1.0 ) ;
1453 if ( tab.nrow() == 0 )
1454 return tstr ;
1455 uInt scnrow = tab.nrow() ;
1456 ROScalarMeasColumn<MEpoch> scTimeCol( tab, "TIME" ) ;
1457 ROScalarQuantColumn<Double> scIntervalCol( tab, "INTERVAL" ) ;
1458 uInt idx = 0 ;
1459 const Double half = 0.5e0 ;
1460 for ( uInt i = 0 ; i < nrow ; i++ ) {
1461 Double t = tcol( i ).get( "s" ).getValue() ;
1462 for ( uInt j = idx ; j < scnrow-1 ; j++ ) {
1463 Double tsc1 = scTimeCol( j ).get( "s" ).getValue() ;
1464 Double dt1 = scIntervalCol( j ).getValue("s") ;
1465 Double tsc2 = scTimeCol( j+1 ).get( "s" ).getValue() ;
1466 Double dt2 = scIntervalCol( j+1 ).getValue("s") ;
1467 if ( t > tsc1-half*dt1 && t <= tsc2-half*dt2 ) {
1468 tstr[i] = tsc1 ;
1469 idx = j ;
1470 break ;
1471 }
1472 }
1473 if ( tstr[i] == -1.0 ) {
1474 Double tsc = scTimeCol( scnrow-1 ).get( "s" ).getValue() ;
1475 Double dt = scIntervalCol( scnrow-1 ).getValue( "s" ) ;
1476 if ( t <= tsc+0.5*dt )
1477 tstr[i] = tsc ;
1478 }
1479 }
1480 //double endSec = gettimeofday_sec() ;
1481 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1482 return tstr ;
1483}
1484
1485uInt MSFiller::getTsys( uInt idx, Matrix<Float> &tsys, MSSysCal &tab, Double t )
1486{
1487 //double startSec = gettimeofday_sec() ;
1488 //os_ << "start MSFiller::getTsys() startSec=" << startSec << LogIO::POST ;
1489 uInt nrow = tab.nrow() ;
1490 if ( nrow == 0 ) {
1491 os_ << "No SysCal rows" << LogIO::POST ;
1492 tsys.resize( IPosition(0) ) ;
1493 return 0 ;
1494 }
1495 Bool isSp = tab.tableDesc().isColumn( "TSYS_SPECTRUM" ) ;
1496 ROScalarMeasColumn<MEpoch> scTimeCol( tab, "TIME" ) ;
1497 ROArrayColumn<Float> mTsysCol ;
1498 if ( isSp ) {
1499 mTsysCol.attach( tab, "TSYS_SPECTRUM" ) ;
1500 }
1501 else {
1502 mTsysCol.attach( tab, "TSYS" ) ;
1503 }
1504 for ( uInt i = idx ; i < nrow ; i++ ) {
1505 Double tref = scTimeCol( i ).get( "s" ).getValue() ;
1506 if ( t == tref ) {
1507 tsys.reference( mTsysCol( i ) ) ;
1508 idx = i ;
1509 break ;
1510 }
1511 }
1512 //os_ << "MSFiller::getTsys() idx = " << idx << " tsys = " << tsys << LogIO::POST ;
1513 //double endSec = gettimeofday_sec() ;
1514 //os_ << "end MSFiller::getTsys() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1515 return idx ;
1516}
1517
1518Vector<uInt> MSFiller::getTcalId( Int fid, Int spwid, Double t )
1519{
1520 //double startSec = gettimeofday_sec() ;
1521 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;
1522 if ( table_->tcal().table().nrow() == 0 ) {
1523 os_ << "No TCAL rows" << LogIO::POST ;
1524 Vector<uInt> tcalids( 0 ) ;
1525 return tcalids ;
1526 }
1527 String feed = "FEED" + String::toString(fid) ;
1528 String spw = "SPW" + String::toString(spwid) ;
1529 String sctime = MVTime( Quantum<Double>(t,"s") ).string(MVTime::YMD) ;
1530 String key = feed + ":" + spw + ":" + sctime ;
1531 //Record rec = tcalrec_.asRecord(feed).asRecord(spw) ;
1532 //if ( !rec.isDefined( sctime ) ) {
1533 if ( !tcalrec_.isDefined( key ) ) {
1534 os_ << "No TCAL rows" << LogIO::POST ;
1535 Vector<uInt> tcalids( 0 ) ;
1536 return tcalids ;
1537 }
1538 //Vector<uInt> ids = rec.asArrayuInt(sctime) ;
1539 Vector<uInt> ids = tcalrec_.asArrayuInt( key ) ;
1540 uInt npol = ids[1] - ids[0] + 1 ;
1541 Vector<uInt> tcalids( npol ) ;
1542 tcalids[0] = ids[0] ;
1543 tcalids[1] = ids[1] ;
1544 for ( uInt ipol = 2 ; ipol < npol ; ipol++ )
1545 tcalids[ipol] = ids[0] + ipol - 1 ;
1546
1547 //double endSec = gettimeofday_sec() ;
1548 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1549 return tcalids ;
1550}
1551
1552uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t )
1553{
1554 //double startSec = gettimeofday_sec() ;
1555 //os_ << "start MSFiller::getDirection() startSec=" << startSec << LogIO::POST ;
1556 // assume that cols is sorted by TIME
1557 Bool doInterp = False ;
1558 //uInt nrow = cols.nrow() ;
1559 uInt nrow = tab.nrow() ;
1560 if ( nrow == 0 )
1561 return 0 ;
1562 ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ;
1563 ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ;
1564 ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ;
1565 // ensure that tcol(idx) < t
1566 //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ;
1567 while ( tcol(idx).get("s").getValue() > t && idx > 0 )
1568 idx-- ;
1569 //os_ << "idx = " << idx << LogIO::POST ;
1570
1571 // index search
1572 for ( uInt i = idx ; i < nrow ; i++ ) {
1573 Double tref = tcol( i ).get( "s" ).getValue() ;
1574 if ( tref == t ) {
1575 idx = i ;
1576 break ;
1577 }
1578 else if ( tref > t ) {
1579 if ( i == 0 ) {
1580 idx = i ;
1581 }
1582 else {
1583 idx = i-1 ;
1584 doInterp = True ;
1585 }
1586 break ;
1587 }
1588 else {
1589 idx = nrow - 1 ;
1590 }
1591 }
1592 //os_ << "searched idx = " << idx << LogIO::POST ;
1593
1594 //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ;
1595 IPosition ip( dmcol(idx).shape().nelements(), 0 ) ;
1596 //os_ << "ip = " << ip << LogIO::POST ;
1597 ref = dmcol(idx)(ip).getRefString() ;
1598 //os_ << "ref = " << ref << LogIO::POST ;
1599 if ( doInterp ) {
1600 //os_ << "do interpolation" << LogIO::POST ;
1601 //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ;
1602 Double tref0 = tcol(idx).get("s").getValue() ;
1603 Double tref1 = tcol(idx+1).get("s").getValue() ;
1604 Matrix<Double> mdir0 = dcol( idx ) ;
1605 Matrix<Double> mdir1 = dcol( idx+1 ) ;
1606 Vector<Double> dir0 = mdir0.column( 0 ) ;
1607 //os_ << "dir0 = " << dir0 << LogIO::POST ;
1608 Vector<Double> dir1 = mdir1.column( 0 ) ;
1609 //os_ << "dir1 = " << dir1 << LogIO::POST ;
1610 Double dt0 = t - tref0 ;
1611 Double dt1 = tref1 - t ;
1612 dir.reference( (dt0*dir1+dt1*dir0)/(dt0+dt1) ) ;
1613 if ( mdir0.ncolumn() > 1 ) {
1614 if ( dt0 >= dt1 )
1615 srate.reference( mdir0.column( 1 ) ) ;
1616 else
1617 srate.reference( mdir1.column( 1 ) ) ;
1618 }
1619 //os_ << "dir = " << dir << LogIO::POST ;
1620 }
1621 else {
1622 //os_ << "no interpolation" << LogIO::POST ;
1623 Matrix<Double> mdir0 = dcol( idx ) ;
1624 dir.reference( mdir0.column( 0 ) ) ;
1625 if ( mdir0.ncolumn() > 1 )
1626 srate.reference( mdir0.column( 1 ) ) ;
1627 }
1628
1629 //double endSec = gettimeofday_sec() ;
1630 //os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;
1631 return idx ;
1632}
1633
1634} ;
1635
Note: See TracBrowser for help on using the repository browser.