Changeset 2240 for branches/parallel
- Timestamp:
- 07/22/11 12:26:10 (13 years ago)
- Location:
- branches/parallel
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/parallel
- Property svn:mergeinfo changed
/trunk merged: 2217-2228,2231-2239
- Property svn:mergeinfo changed
-
branches/parallel/Makefile
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/SConstruct
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/apps
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/external-alma
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/external-alma/atnf/pks/pks_maths.cc
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/getsvnrev.sh
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/python
- Property svn:mergeinfo changed (with no actual effect on merging)
-
branches/parallel/python/flagtoolbar.py
- Property svn:mergeinfo changed
/trunk/python/flagtoolbar.py merged: 2217,2219,2221-2224,2226,2228,2231-2232,2234,2237-2239
- Property svn:mergeinfo changed
-
branches/parallel/src
- Property svn:mergeinfo changed
/trunk/src merged: 2217,2219,2221-2224,2226,2228,2231-2232,2234,2237-2239
- Property svn:mergeinfo changed
-
branches/parallel/src/MSFiller.cpp
r2191 r2240 95 95 { 96 96 os_.origin( LogOrigin( "MSFiller", "open()", WHERE ) ) ; 97 //double startSec = gettimeofday_sec() ;98 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ;97 //double startSec = gettimeofday_sec() ; 98 //os_ << "start MSFiller::open() startsec=" << startSec << LogIO::POST ; 99 99 //os_ << " filename = " << filename << endl ; 100 100 … … 147 147 isData_ = mstable_.tableDesc().isColumn( "DATA" ) ; 148 148 149 //double endSec = gettimeofday_sec() ;150 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;149 //double endSec = gettimeofday_sec() ; 150 //os_ << "end MSFiller::open() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 151 151 return true ; 152 152 } … … 155 155 { 156 156 os_.origin( LogOrigin( "MSFiller", "fill()", WHERE ) ) ; 157 //double startSec = gettimeofday_sec() ;158 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ;159 160 //double time0 = gettimeofday_sec() ;161 //os_ << "start init fill: " << time0 << LogIO::POST ;157 //double startSec = gettimeofday_sec() ; 158 //os_ << "start MSFiller::fill() startSec=" << startSec << LogIO::POST ; 159 160 //double time0 = gettimeofday_sec() ; 161 //os_ << "start init fill: " << time0 << LogIO::POST ; 162 162 163 163 // Initialize header 164 164 STHeader sdh ; 165 sdh.nchan = 0 ; 166 sdh.npol = 0 ; 167 sdh.nif = 0 ; 168 sdh.nbeam = 0 ; 169 sdh.observer = "" ; 170 sdh.project = "" ; 171 sdh.obstype = "" ; 172 sdh.antennaname = "" ; 173 sdh.antennaposition.resize( 0 ) ; 174 sdh.equinox = 0.0 ; 175 sdh.freqref = "" ; 176 sdh.reffreq = -1.0 ; 177 sdh.bandwidth = 0.0 ; 178 sdh.utc = 0.0 ; 179 sdh.fluxunit = "" ; 180 sdh.epoch = "" ; 181 sdh.poltype = "" ; 165 initHeader( sdh ) ; 182 166 183 167 // check if optional table exists … … 249 233 250 234 // SUBTABLES: FREQUENCIES 251 table_->frequencies().setFrame( "LSRK" ) ; 252 table_->frequencies().setFrame( "LSRK", True ) ; 235 //string freqFrame = getFrame() ; 236 string freqFrame = "LSRK" ; 237 table_->frequencies().setFrame( freqFrame ) ; 238 table_->frequencies().setFrame( freqFrame, True ) ; 253 239 254 240 // SUBTABLES: WEATHER … … 266 252 // SUBTABLES: HISTORY 267 253 //fillHistory() ; 268 269 // shared pointers270 ROTableColumn *tcolr ;271 254 272 255 // MAIN … … 278 261 MPosition mp( MVPosition( antpos ), MPosition::ITRF ) ; 279 262 if ( getPt_ ) { 280 //pointtab = pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ;281 263 pointtab = MSPointing( pointtab( pointtab.col("ANTENNA_ID")==antenna_ ).sort("TIME") ) ; 282 264 } 283 tcolr = tpoolr->construct( anttab, "STATION" ) ; 284 String stationName = tcolr->asString( antenna_ ) ; 285 tpoolr->destroy( tcolr ) ; 286 tcolr = tpoolr->construct( anttab, "NAME" ) ; 287 String antennaName = tcolr->asString( antenna_ ) ; 288 tpoolr->destroy( tcolr ) ; 265 ROScalarColumn<Double> ptcol( pointtab, "TIME" ) ; 266 ROArrayColumn<Double> pdcol( pointtab, "DIRECTION" ) ; 267 String stationName = asString( "STATION", antenna_, anttab, tpoolr ) ; 268 String antennaName = asString( "NAME", antenna_, anttab, tpoolr ) ; 289 269 sdh.antennaposition.resize( 3 ) ; 290 270 for ( int i = 0 ; i < 3 ; i++ ) … … 292 272 String telescopeName = "" ; 293 273 294 //double time1 = gettimeofday_sec() ;295 // os_ << "end fill init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;274 //double time1 = gettimeofday_sec() ; 275 //os_ << "end init fill: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 296 276 297 277 // row based … … 313 293 RecordFieldPtr<uInt> widRF( trec, "WEATHER_ID" ) ; 314 294 RecordFieldPtr<uInt> polnoRF( trec, "POLNO" ) ; 315 295 RecordFieldPtr<Int> refbRF( trec, "REFBEAMNO" ) ; 296 RecordFieldPtr<Int> fitidRF( trec, "FIT_ID" ) ; 297 RecordFieldPtr<Float> tauRF( trec, "OPACITY" ) ; 298 RecordFieldPtr<uInt> beamRF( trec, "BEAMNO" ) ; 299 RecordFieldPtr<uInt> focusidRF( trec, "FOCUS_ID" ) ; 300 RecordFieldPtr<uInt> ifnoRF( trec, "IFNO" ) ; 301 RecordFieldPtr<uInt> molidRF( trec, "MOLECULE_ID" ) ; 302 RecordFieldPtr<uInt> freqidRF( trec, "FREQ_ID" ) ; 303 RecordFieldPtr<uInt> scanRF( trec, "SCANNO" ) ; 304 RecordFieldPtr<Int> srctypeRF( trec, "SRCTYPE" ) ; 305 RecordFieldPtr<String> srcnameRF( trec, "SRCNAME" ) ; 306 RecordFieldPtr<String> fieldnameRF( trec, "FIELDNAME" ) ; 307 RecordFieldPtr< Array<Double> > srcpmRF( trec, "SRCPROPERMOTION" ) ; 308 RecordFieldPtr< Array<Double> > srcdirRF( trec, "SRCDIRECTION" ) ; 309 RecordFieldPtr<Double> sysvelRF( trec, "SRCVELOCITY" ) ; 316 310 317 311 // REFBEAMNO 318 RecordFieldPtr<Int> intRF( trec, "REFBEAMNO" ) ; 319 *intRF = 0 ; 312 *refbRF = -1 ; 320 313 321 314 // FIT_ID 322 intRF.attachToRecord( trec, "FIT_ID" ) ; 323 *intRF = -1 ; 315 *fitidRF = -1 ; 324 316 325 317 // OPACITY 326 RecordFieldPtr<Float> floatRF( trec, "OPACITY" ) ; 327 *floatRF = 0.0 ; 318 *tauRF = 0.0 ; 328 319 329 320 // … … 332 323 TableIterator iter0( mstable_, "OBSERVATION_ID" ) ; 333 324 while( !iter0.pastEnd() ) { 334 //time0 = gettimeofday_sec() ;335 //os_ << "start 0th iteration: " << time0 << LogIO::POST ;325 //time0 = gettimeofday_sec() ; 326 //os_ << "start 0th iteration: " << time0 << LogIO::POST ; 336 327 Table t0 = iter0.table() ; 337 tcolr = tpoolr->construct( t0, "OBSERVATION_ID" ) ; 338 Int obsId = tcolr->asInt( 0 ) ; 339 tpoolr->destroy( tcolr ) ; 328 Int obsId = asInt( "OBSERVATION_ID", 0, t0, tpoolr ) ; 340 329 if ( sdh.observer == "" ) { 341 tcolr = tpoolr->construct( obstab, "OBSERVER" ) ; 342 sdh.observer = tcolr->asString( obsId ) ; 343 tpoolr->destroy( tcolr ) ; 330 sdh.observer = asString( "OBSERVER", obsId, obstab, tpoolr ) ; 344 331 } 345 332 if ( sdh.project == "" ) { 346 tcolr = tpoolr->construct( obstab, "PROJECT" ) ; 347 sdh.observer = tcolr->asString( obsId ) ; 348 tpoolr->destroy( tcolr ) ; 333 sdh.project = asString( "PROJECT", obsId, obstab, tpoolr ) ; 349 334 } 350 335 ROArrayMeasColumn<MEpoch> *tmpMeasCol = new ROArrayMeasColumn<MEpoch>( obstab, "TIME_RANGE" ) ; … … 355 340 } 356 341 if ( telescopeName == "" ) { 357 tcolr = tpoolr->construct( obstab, "TELESCOPE_NAME" ) ; 358 telescopeName = tcolr->asString( obsId ) ; 359 tpoolr->destroy( tcolr ) ; 342 telescopeName = asString( "TELESCOPE_NAME", obsId, obstab, tpoolr ) ; 360 343 } 361 344 Int nbeam = 0 ; 362 //time1 = gettimeofday_sec() ;363 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;345 //time1 = gettimeofday_sec() ; 346 //os_ << "end 0th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 364 347 // 365 348 // ITERATION: FEED1 … … 367 350 TableIterator iter1( t0, "FEED1" ) ; 368 351 while( !iter1.pastEnd() ) { 369 //time0 = gettimeofday_sec() ;370 //os_ << "start 1st iteration: " << time0 << LogIO::POST ;352 //time0 = gettimeofday_sec() ; 353 //os_ << "start 1st iteration: " << time0 << LogIO::POST ; 371 354 Table t1 = iter1.table() ; 372 355 // assume FEED1 == FEED2 373 tcolr = tpoolr->construct( t1, "FEED1" ) ; 374 Int feedId = tcolr->asInt( 0 ) ; 375 tpoolr->destroy( tcolr ) ; 356 Int feedId = asInt( "FEED1", 0, t1, tpoolr ) ; 376 357 nbeam++ ; 377 358 378 359 // BEAMNO 379 RecordFieldPtr<uInt> uintRF( trec, "BEAMNO" ) ; 380 *uintRF = feedId ; 360 *beamRF = feedId ; 381 361 382 362 // FOCUS_ID 383 uintRF.attachToRecord( trec, "FOCUS_ID" ) ; 384 *uintRF = 0 ; 385 386 // time1 = gettimeofday_sec() ; 387 // os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 363 *focusidRF = 0 ; 364 365 //time1 = gettimeofday_sec() ; 366 //os_ << "end 1st iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 388 367 // 389 368 // ITERATION: FIELD_ID … … 391 370 TableIterator iter2( t1, "FIELD_ID" ) ; 392 371 while( !iter2.pastEnd() ) { 393 //time0 = gettimeofday_sec() ;394 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ;372 //time0 = gettimeofday_sec() ; 373 //os_ << "start 2nd iteration: " << time0 << LogIO::POST ; 395 374 Table t2 = iter2.table() ; 396 tcolr = tpoolr->construct( t2, "FIELD_ID" ) ; 397 Int fieldId = tcolr->asInt( 0 ) ; 398 tpoolr->destroy( tcolr ) ; 399 tcolr = tpoolr->construct( fieldtab, "SOURCE_ID" ) ; 400 Int srcId = tcolr->asInt( fieldId ) ; 401 tpoolr->destroy( tcolr ) ; 402 tcolr = tpoolr->construct( fieldtab, "NAME" ) ; 403 String fieldName = tcolr->asString( fieldId ) + "__" + String::toString(fieldId) ; 404 tpoolr->destroy( tcolr ) ; 375 Int fieldId = asInt( "FIELD_ID", 0, t2, tpoolr ) ; 376 Int srcId = asInt( "SOURCE_ID", fieldId, fieldtab, tpoolr ) ; 377 String fieldName = asString( "NAME", fieldId, fieldtab, tpoolr ) ; 378 fieldName += "__" + String::toString(fieldId) ; 405 379 ROArrayMeasColumn<MDirection> *delayDirCol = new ROArrayMeasColumn<MDirection>( fieldtab, "DELAY_DIR" ) ; 406 380 Vector<MDirection> delayDir = (*delayDirCol)( fieldId ) ; 407 delete delayDirCol ; 408 Vector<Double> defaultScanrate( 2, 0.0 ) ; 409 Vector<Double> defaultDir = delayDir[0].getAngle( "rad" ).getValue() ; 410 if ( delayDir.nelements() > 1 ) 411 defaultScanrate = delayDir[1].getAngle( "rad" ).getValue() ; 412 381 delete delayDirCol ; 413 382 414 383 // FIELDNAME 415 RecordFieldPtr<String> strRF( trec, "FIELDNAME" ) ; 416 *strRF = fieldName ; 417 418 419 // time1 = gettimeofday_sec() ; 420 // os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 384 *fieldnameRF = fieldName ; 385 386 387 //time1 = gettimeofday_sec() ; 388 //os_ << "end 2nd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 421 389 // 422 390 // ITERATION: DATA_DESC_ID … … 424 392 TableIterator iter3( t2, "DATA_DESC_ID" ) ; 425 393 while( !iter3.pastEnd() ) { 426 //time0 = gettimeofday_sec() ;427 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ;394 //time0 = gettimeofday_sec() ; 395 //os_ << "start 3rd iteration: " << time0 << LogIO::POST ; 428 396 Table t3 = iter3.table() ; 429 tcolr = tpoolr->construct( t3, "DATA_DESC_ID" ) ; 430 Int ddId = tcolr->asInt( 0 ) ; 431 tpoolr->destroy( tcolr ) ; 432 tcolr = tpoolr->construct( ddtab, "POLARIZATION_ID" ) ; 433 Int polId = tcolr->asInt( ddId ) ; 434 tpoolr->destroy( tcolr ) ; 435 tcolr = tpoolr->construct( ddtab, "SPECTRAL_WINDOW_ID" ) ; 436 Int spwId = tcolr->asInt( ddId ) ; 437 tpoolr->destroy( tcolr ) ; 397 Int ddId = asInt( "DATA_DESC_ID", 0, t3, tpoolr ) ; 398 Int polId = asInt( "POLARIZATION_ID", ddId, ddtab, tpoolr ) ; 399 Int spwId = asInt( "SPECTRAL_WINDOW_ID", ddId, ddtab, tpoolr ) ; 438 400 439 401 // IFNO 440 uintRF.attachToRecord( trec, "IFNO" ) ; 441 *uintRF = (uInt)spwId ; 402 *ifnoRF = (uInt)spwId ; 442 403 443 404 // polarization information 444 tcolr = tpoolr->construct( poltab, "NUM_CORR" ) ; 445 Int npol = tcolr->asInt( polId ) ; 446 tpoolr->destroy( tcolr ) ; 405 Int npol = asInt( "NUM_CORR", polId, poltab, tpoolr ) ; 447 406 ROArrayColumn<Int> *roArrICol = new ROArrayColumn<Int>( poltab, "CORR_TYPE" ) ; 448 407 Vector<Int> corrtype = (*roArrICol)( polId ) ; 449 408 delete roArrICol ; 450 // os_ << "npol = " << npol << LogIO::POST ;451 // os_ << "corrtype = " << corrtype << LogIO::POST ;452 // source information453 // os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ;454 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ;455 if ( srctabSel.nrow() == 0 ) {456 srctabSel = srctab( srctab.col("SOURCE_ID") == srcId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ;457 }458 409 String srcName( "" ) ; 459 410 Vector<Double> srcPM( 2, 0.0 ) ; 460 411 Vector<Double> srcDir( 2, 0.0 ) ; 461 412 MDirection md ; 462 Int numLines = 0 ; 463 ROArrayColumn<Double> *roArrDCol = 0 ; 464 if ( srctabSel.nrow() > 0 ) { 465 // source name 466 tcolr = tpoolr->construct( srctabSel, "NAME" ) ; 467 srcName = tcolr->asString( 0 ) ; 468 tpoolr->destroy( tcolr ) ; 469 470 // source proper motion 471 roArrDCol = new ROArrayColumn<Double>( srctabSel, "PROPER_MOTION" ) ; 472 srcPM = (*roArrDCol)( 0 ) ; 473 delete roArrDCol ; 474 475 // source direction 476 roArrDCol = new ROArrayColumn<Double>( srctabSel, "DIRECTION" ) ; 477 srcDir = (*roArrDCol)( 0 ) ; 478 delete roArrDCol ; 479 480 // source direction as MDirection object 481 ROScalarMeasColumn<MDirection> *tmpMeasCol = new ROScalarMeasColumn<MDirection>( srctabSel, "DIRECTION" ) ; 482 md = (*tmpMeasCol)( 0 ) ; 483 delete tmpMeasCol ; 484 485 // number of lines 486 tcolr = tpoolr->construct( srctabSel, "NUM_LINES" ) ; 487 numLines = tcolr->asInt( 0 ) ; 488 tpoolr->destroy( tcolr ) ; 489 490 } 491 else { 492 md = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 493 } 413 Vector<Double> restFreqs ; 414 Vector<String> transitionName ; 415 Vector<Double> sysVels ; 416 // os_ << "npol = " << npol << LogIO::POST ; 417 // os_ << "corrtype = " << corrtype << LogIO::POST ; 418 419 // source and molecular transition 420 sourceInfo( srcId, spwId, srcName, md, srcPM, restFreqs, transitionName, sysVels, tpoolr ) ; 421 // os_ << "srcId = " << srcId << ", spwId = " << spwId << LogIO::POST ; 494 422 495 423 // SRCNAME 496 strRF.attachToRecord( trec, "SRCNAME" ) ; 497 *strRF = srcName ; 424 *srcnameRF = srcName ; 498 425 499 426 // os_ << "srcName = " << srcName << LogIO::POST ; 500 427 501 428 // SRCPROPERMOTION 502 RecordFieldPtr< Array<Double> > darrRF( trec, "SRCPROPERMOTION" ) ; 503 *darrRF = srcPM ; 429 *srcpmRF = srcPM ; 504 430 505 431 //os_ << "srcPM = " << srcPM << LogIO::POST ; 506 432 507 433 // SRCDIRECTION 508 darrRF.attachToRecord( trec, "SRCDIRECTION" ) ; 509 *darrRF = srcDir ; 434 *srcdirRF = md.getAngle().getValue( "rad" ) ; 510 435 511 436 //os_ << "srcDir = " << srcDir << LogIO::POST ; 512 437 513 // for MOLECULES subtable 514 // os_ << "numLines = " << numLines << LogIO::POST ; 515 516 Vector<Double> restFreqs( numLines, 0.0 ) ; 517 Vector<String> transitionName( numLines, "" ) ; 518 Vector<Double> sysVels ; 438 // SRCVELOCITY 519 439 Double sysVel = 0.0 ; 520 if ( numLines != 0 ) { 521 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) { 522 sharedQDArrCol = new ROArrayQuantColumn<Double>( srctabSel, "REST_FREQUENCY" ) ; 523 Array< Quantum<Double> > qRestFreqs = (*sharedQDArrCol)( 0 ) ; 524 delete sharedQDArrCol ; 525 for ( int i = 0 ; i < numLines ; i++ ) { 526 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ; 527 } 528 } 529 // os_ << "restFreqs = " << restFreqs << LogIO::POST ; 530 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) { 531 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ; 532 if ( transitionCol.isDefined( 0 ) ) 533 transitionName = transitionCol( 0 ) ; 534 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ; 535 } 536 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) { 537 roArrDCol = new ROArrayColumn<Double>( srctabSel, "SYSVEL" ) ; 538 sysVels = (*roArrDCol)( 0 ) ; 539 delete roArrDCol ; 540 } 541 if ( !sysVels.empty() ) { 542 //os_ << "sysVels.shape() = " << sysVels.shape() << LogIO::POST ; 543 // NB: assume all SYSVEL values are the same 544 sysVel = sysVels( IPosition(1,0) ) ; 545 } 546 } 547 548 // SRCVELOCITY 549 RecordFieldPtr<Double> doubleRF( trec, "SRCVELOCITY" ) ; 550 *doubleRF = sysVel ; 440 if ( !sysVels.empty() ) 441 sysVel = sysVels[0] ; 442 *sysvelRF = sysVel ; 551 443 552 444 // os_ << "sysVel = " << sysVel << LogIO::POST ; … … 555 447 556 448 // MOLECULE_ID 557 uintRF.attachToRecord( trec, "MOLECULE_ID" ) ; 558 *uintRF = molId ; 449 *molidRF = molId ; 559 450 560 451 // spectral setup 561 MeasFrame mf( me, mp, md ) ; 562 tcolr = tpoolr->construct( spwtab, "MEAS_FREQ_REF" ) ; 563 MFrequency::Types freqRef = MFrequency::castType((uInt)(tcolr->asInt(spwId))) ; 564 tpoolr->destroy( tcolr ) ; 565 tcolr = tpoolr->construct( spwtab, "NUM_CHAN" ) ; 566 Int nchan = tcolr->asInt( spwId ) ; 452 uInt freqId ; 453 Int nchan = asInt( "NUM_CHAN", spwId, spwtab, tpoolr ) ; 567 454 Bool iswvr = False ; 568 455 if ( nchan == 4 ) iswvr = True ; 569 tpoolr->destroy( tcolr ) ;570 Bool even = False ;571 if ( (nchan/2)*2 == nchan ) even = True ;572 456 sdh.nchan = max( sdh.nchan, nchan ) ; 573 ROScalarQuantColumn<Double> *tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "TOTAL_BANDWIDTH" ) ; 574 Double totbw = (*tmpQuantCol)( spwId ).getValue( "Hz" ) ; 575 delete tmpQuantCol ; 576 sdh.bandwidth = max( sdh.bandwidth, totbw ) ; 577 if ( sdh.freqref == "" ) 578 //sdh.freqref = MFrequency::showType( freqRef ) ; 579 sdh.freqref = "LSRK" ; 580 if ( sdh.reffreq == -1.0 ) { 581 tmpQuantCol = new ROScalarQuantColumn<Double>( spwtab, "REF_FREQUENCY" ) ; 582 Quantum<Double> qreffreq = (*tmpQuantCol)( spwId ) ; 583 delete tmpQuantCol ; 584 if ( freqRef == MFrequency::LSRK ) { 585 sdh.reffreq = qreffreq.getValue("Hz") ; 586 } 587 else { 588 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 589 sdh.reffreq = tolsr( qreffreq ).get("Hz").getValue() ; 590 } 591 } 592 Int refchan = nchan / 2 ; 593 IPosition refip( 1, refchan ) ; 594 Double refpix = 0.5*(nchan-1) ; 595 Double refval = 0.0 ; 596 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_WIDTH" ) ; 597 Double increment = (*sharedQDArrCol)( spwId )( refip ).getValue( "Hz" ) ; 598 delete sharedQDArrCol ; 599 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ; 600 sharedQDArrCol = new ROArrayQuantColumn<Double>( spwtab, "CHAN_FREQ" ) ; 601 Vector< Quantum<Double> > chanFreqs = (*sharedQDArrCol)( spwId ) ; 602 delete sharedQDArrCol ; 603 if ( freqRef == MFrequency::LSRK ) { 604 if ( even ) { 605 IPosition refip0( 1, refchan-1 ) ; 606 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 607 Double refval1 = chanFreqs(refip).getValue("Hz") ; 608 refval = 0.5 * ( refval0 + refval1 ) ; 609 } 610 else { 611 refval = chanFreqs(refip).getValue("Hz") ; 612 } 613 } 614 else { 615 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 616 if ( even ) { 617 IPosition refip0( 1, refchan-1 ) ; 618 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 619 Double refval1 = chanFreqs(refip).getValue("Hz") ; 620 refval = 0.5 * ( refval0 + refval1 ) ; 621 refval = tolsr( refval ).get("Hz").getValue() ; 622 } 623 else { 624 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ; 625 } 626 } 627 uInt freqId = table_->frequencies().addEntry( refpix, refval, increment ) ; 628 if ( ifmap.find( spwId ) == ifmap.end() ) { 457 map<Int,uInt>::iterator iter = ifmap.find( spwId ) ; 458 if ( iter == ifmap.end() ) { 459 ROScalarMeasColumn<MEpoch> *tmpMeasCol = new ROScalarMeasColumn<MEpoch>( t3, "TIME" ) ; 460 me = (*tmpMeasCol)( 0 ) ; 461 delete tmpMeasCol ; 462 Double refpix ; 463 Double refval ; 464 Double increment ; 465 spectralSetup( spwId, 466 me, 467 mp, 468 md, 469 refpix, 470 refval, 471 increment, 472 nchan, 473 sdh.freqref, 474 sdh.reffreq, 475 sdh.bandwidth, 476 tpoolr ) ; 477 freqId = table_->frequencies().addEntry( refpix, refval, increment ) ; 629 478 ifmap.insert( pair<Int, uInt>(spwId,freqId) ) ; 630 479 //os_ << "added to ifmap: (" << spwId << "," << freqId << ")" << LogIO::POST ; 631 480 } 481 else { 482 freqId = iter->second ; 483 } 632 484 633 485 // FREQ_ID 634 uintRF.attachToRecord( trec, "FREQ_ID" ) ; 635 *uintRF = freqId ; 486 *freqidRF = freqId ; 636 487 637 488 // for TSYS and TCAL … … 656 507 if ( !iswvr && sdh.poltype == "" ) sdh.poltype = getPolType( corrtype[0] ) ; 657 508 658 //time1 = gettimeofday_sec() ;659 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;509 //time1 = gettimeofday_sec() ; 510 //os_ << "end 3rd iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 660 511 // 661 512 // ITERATION: SCAN_NUMBER … … 663 514 TableIterator iter4( t3, "SCAN_NUMBER" ) ; 664 515 while( !iter4.pastEnd() ) { 665 //time0 = gettimeofday_sec() ;666 //os_ << "start 4th iteration: " << time0 << LogIO::POST ;516 //time0 = gettimeofday_sec() ; 517 //os_ << "start 4th iteration: " << time0 << LogIO::POST ; 667 518 Table t4 = iter4.table() ; 668 tcolr = tpoolr->construct( t4, "SCAN_NUMBER" ) ; 669 Int scanNum = tcolr->asInt( 0 ) ; 670 tpoolr->destroy( tcolr ) ; 519 Int scanNum = asInt( "SCAN_NUMBER", 0, t4, tpoolr ) ; 671 520 672 521 // SCANNO 673 uintRF.attachToRecord( trec, "SCANNO" ) ; 674 *uintRF = scanNum - 1 ; 675 676 // time1 = gettimeofday_sec() ; 677 // os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 522 *scanRF = scanNum - 1 ; 523 524 //time1 = gettimeofday_sec() ; 525 //os_ << "end 4th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 678 526 // 679 527 // ITERATION: STATE_ID … … 681 529 TableIterator iter5( t4, "STATE_ID" ) ; 682 530 while( !iter5.pastEnd() ) { 683 //time0 = gettimeofday_sec() ;684 //os_ << "start 5th iteration: " << time0 << LogIO::POST ;531 //time0 = gettimeofday_sec() ; 532 //os_ << "start 5th iteration: " << time0 << LogIO::POST ; 685 533 Table t5 = iter5.table() ; 686 tcolr = tpoolr->construct( t5, "STATE_ID" ) ; 687 Int stateId = tcolr->asInt( 0 ) ; 688 tpoolr->destroy( tcolr ) ; 689 tcolr = tpoolr->construct( stattab, "OBS_MODE" ) ; 690 String obstype = tcolr->asString( stateId ) ; 691 tpoolr->destroy( tcolr ) ; 534 Int stateId = asInt( "STATE_ID", 0, t5, tpoolr ) ; 535 String obstype = asString( "OBS_MODE", 0, stattab, tpoolr ) ; 692 536 if ( sdh.obstype == "" ) sdh.obstype = obstype ; 693 537 694 538 Int nrow = t5.nrow() ; 695 //time1 = gettimeofday_sec() ;696 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ;539 //time1 = gettimeofday_sec() ; 540 //os_ << "end 5th iteration init: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 697 541 698 542 uInt cycle = 0 ; … … 700 544 Cube<Float> spArr ; 701 545 Cube<Bool> flArr ; 702 if ( isFloatData_ ) { 703 ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ; 704 ROArrayColumn<Float> mFloatDataCol( t5, "FLOAT_DATA" ) ; 705 spArr = mFloatDataCol.getColumn() ; 706 flArr = mFlagCol.getColumn() ; 707 if ( sdh.fluxunit == "" ) { 708 const TableRecord &dataColKeys = mFloatDataCol.keywordSet() ; 709 if ( dataColKeys.isDefined( "UNIT" ) ) 710 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 711 } 546 reshapeSpectraAndFlagtra( spArr, 547 flArr, 548 t5, 549 npol, 550 nchan, 551 nrow, 552 corrtype ) ; 553 if ( sdh.fluxunit == "" ) { 554 String colName = "FLOAT_DATA" ; 555 if ( isData_ ) colName = "DATA" ; 556 ROTableColumn dataCol( t5, colName ) ; 557 const TableRecord &dataColKeys = dataCol.keywordSet() ; 558 if ( dataColKeys.isDefined( "UNIT" ) ) 559 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 560 else if ( dataColKeys.isDefined( "QuantumUnits" ) ) 561 sdh.fluxunit = dataColKeys.asString( "QuantumUnits" ) ; 712 562 } 713 else if ( isData_ ) { 714 spArr.resize( npol, nchan, nrow ) ; 715 flArr.resize( npol, nchan, nrow ) ; 716 ROArrayColumn<Bool> mFlagCol( t5, "FLAG" ) ; 717 ROArrayColumn<Complex> mDataCol( t5, "DATA" ) ; 718 for ( Int irow = 0 ; irow < nrow ; irow++ ) { 719 Bool crossOK = False ; 720 Matrix<Complex> sp = mDataCol( irow ) ; 721 Matrix<Bool> fl = mFlagCol( irow ) ; 722 Matrix<Float> spxy = spArr.xyPlane( irow ) ; 723 Matrix<Bool> flxy = flArr.xyPlane( irow ) ; 724 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 725 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX 726 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) { 727 if ( !crossOK ) { 728 spxy.row( ipol ) = real( sp.row( ipol ) ) ; 729 flxy.row( ipol ) = fl.row( ipol ) ; 730 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::RL ) { 731 spxy.row( ipol+1 ) = imag( sp.row( ipol ) ) ; 732 flxy.row( ipol+1 ) = fl.row( ipol ) ; 733 } 734 else { 735 spxy.row( ipol+1 ) = imag( conj( sp.row( ipol ) ) ) ; 736 flxy.row( ipol+1 ) = fl.row( ipol ) ; 737 } 738 crossOK = True ; 739 } 740 } 741 else { 742 spxy.row( ipol ) = real( sp.row( ipol ) ) ; 743 flxy.row( ipol ) = fl.row( ipol ) ; 744 } 745 } 746 } 747 if ( sdh.fluxunit == "" ) { 748 const TableRecord &dataColKeys = mDataCol.keywordSet() ; 749 if ( dataColKeys.isDefined( "UNIT" ) ) 750 sdh.fluxunit = dataColKeys.asString( "UNIT" ) ; 751 } 752 } 563 753 564 ROScalarMeasColumn<MEpoch> *mTimeCol = new ROScalarMeasColumn<MEpoch>( t5, "TIME" ) ; 754 565 Block<MEpoch> mTimeB( nrow ) ; 755 566 for ( Int irow = 0 ; irow < nrow ; irow++ ) 756 567 mTimeB[irow] = (*mTimeCol)( irow ) ; 757 ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ;758 ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ;568 // ROTableColumn *mIntervalCol = tpoolr->construct( t5, "INTERVAL" ) ; 569 // ROTableColumn *mFlagRowCol = tpoolr->construct( t5, "FLAG_ROW" ) ; 759 570 Block<Int> sysCalIdx( nrow, -1 ) ; 760 571 if ( isSysCal_ ) { … … 765 576 Int srcType = getSrcType( stateId, tpoolr ) ; 766 577 uInt diridx = 0 ; 767 MDirection::Types dirType ;768 578 uInt wid = 0 ; 769 579 Int pidx = 0 ; … … 789 599 790 600 // SRCTYPE 791 intRF.attachToRecord( trec, "SRCTYPE" ) ; 792 *intRF = srcType ; 601 *srctypeRF = srcType ; 793 602 794 603 for ( Int irow = 0 ; irow < nrow ; irow++ ) { … … 797 606 798 607 // FLAGROW 799 *flrRF = (uInt)mFlagRowCol->asBool( irow ) ; 608 // *flrRF = (uInt)mFlagRowCol->asBool( irow ) ; 609 *flrRF = (uInt)asBool( "FLAG_ROW", irow, t5, tpoolr ) ; 800 610 801 611 // SPECTRA, FLAG … … 809 619 810 620 // INTERVAL 811 *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ; 621 // *intervalRF = (Double)(mIntervalCol->asdouble( irow )) ; 622 *intervalRF = asDouble( "INTERVAL", irow, t5, tpoolr ) ; 812 623 813 624 // TSYS … … 825 636 826 637 // WEATHER_ID 827 if ( isWeather_ ) 638 if ( isWeather_ ) { 828 639 wid = getWeatherId( wid, mTimeB[irow].get("s").getValue() ) ; 829 *widRF = wid ; 830 831 832 // DIRECTION, AZEL, SCANRATE 833 if ( getPt_ ) { 834 Vector<Double> dir ; 835 Vector<Double> scanrate ; 836 String refString ; 837 diridx = getDirection( diridx, dir, scanrate, refString, pointtab, mTimeB[irow].get("s").getValue() ) ; 838 MDirection::getType( dirType, refString ) ; 839 mf.resetEpoch( mTimeB[irow] ) ; 840 mf.resetDirection( MDirection( MVDirection(dir), dirType ) ) ; 841 if ( refString == "J2000" ) { 842 *dirRF = dir ; 843 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 844 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ; 845 *azRF = (Float)azel(0) ; 846 *elRF = (Float)azel(1) ; 847 } 848 else if ( refString(0,4) == "AZEL" ) { 849 *azRF = (Float)dir(0) ; 850 *elRF = (Float)dir(1) ; 851 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 852 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ; 853 *dirRF = newdir ; 854 } 855 else { 856 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 857 Vector<Double> azel = toazel( dir ).getAngle("rad").getValue() ; 858 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 859 Vector<Double> newdir = toj2000( dir ).getAngle("rad").getValue() ; 860 *dirRF = newdir ; 861 *azRF = (Float)azel(0) ; 862 *elRF = (Float)azel(1) ; 863 } 864 if ( scanrate.size() != 0 ) { 865 *scrRF = scanrate ; 866 } 867 else { 868 *scrRF = defaultScanrate ; 869 } 640 *widRF = mwIndex_[wid] ; 870 641 } 871 642 else { 872 String ref = md.getRefString() ; 873 //Vector<Double> defaultDir = srcDir ; 874 MDirection::getType( dirType, "J2000" ) ; 875 if ( ref != "J2000" ) { 876 ROScalarMeasColumn<MEpoch> tmCol( pointtab, "TIME" ) ; 877 mf.resetEpoch( tmCol( 0 ) ) ; 878 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 879 defaultDir = toj2000( defaultDir ).getAngle("rad").getValue() ; 880 } 881 mf.resetEpoch( mTimeB[irow] ) ; 882 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 883 Vector<Double> azel = toazel( defaultDir ).getAngle("rad").getValue() ; 884 *azRF = (Float)azel(0) ; 885 *elRF = (Float)azel(1) ; 886 *dirRF = defaultDir ; 887 *scrRF = defaultScanrate ; 643 *widRF = wid ; 888 644 } 645 646 647 // DIRECTION, AZEL, SCANRATE 648 Vector<Double> dir ; 649 Vector<Double> azel ; 650 Vector<Double> scanrate ; 651 String refString ; 652 if ( getPt_ ) 653 //diridx = getDirection( diridx, dir, azel, scanrate, pointtab, mTimeB[irow], mp ) ; 654 diridx = getDirection( diridx, dir, azel, scanrate, ptcol, pdcol, mTimeB[irow], mp ) ; 655 else 656 getSourceDirection( dir, azel, scanrate, mTimeB[irow], mp, delayDir ) ; 657 *dirRF = dir ; 658 *azRF = azel[0] ; 659 *elRF = azel[1] ; 660 *scrRF = scanrate ; 889 661 890 662 // Polarization dependent things … … 893 665 *polnoRF = polnos[ipol] ; 894 666 895 //*spRF = sp.row( ipol ) ;896 //*ucarrRF = fl.row( ipol ) ;897 //*tsysRF = tsys.row( ipol ) ;898 667 spRF.define( sp.row( ipol ) ) ; 899 668 ucarrRF.define( fl.row( ipol ) ) ; … … 908 677 cycle++ ; 909 678 } 910 tpoolr->destroy( mIntervalCol ) ; 911 tpoolr->destroy( mFlagRowCol ) ; 912 913 // time1 = gettimeofday_sec() ; 914 // os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 679 680 //time1 = gettimeofday_sec() ; 681 //os_ << "end 5th iteration: " << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 915 682 916 683 iter5.next() ; 917 684 } 918 919 920 685 iter4.next() ; 921 686 } 922 923 924 687 iter3.next() ; 925 688 } 926 927 928 689 iter2.next() ; 929 690 } 930 931 932 691 iter1.next() ; 933 692 } … … 974 733 sdh.freqref = "SOURCE"; 975 734 } 735 736 if ( sdh.fluxunit == "" || sdh.fluxunit == "CNTS" ) 737 sdh.fluxunit = "K" ; 976 738 table_->setHeader( sdh ) ; 977 739 … … 1006 768 } 1007 769 1008 //double endSec = gettimeofday_sec() ;1009 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;770 //double endSec = gettimeofday_sec() ; 771 //os_ << "end MSFiller::fill() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1010 772 } 1011 773 … … 1020 782 Int MSFiller::getSrcType( Int stateId, boost::object_pool<ROTableColumn> *tpool ) 1021 783 { 1022 //double startSec = gettimeofday_sec() ;1023 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ;784 //double startSec = gettimeofday_sec() ; 785 //os_ << "start MSFiller::getSrcType() startSec=" << startSec << LogIO::POST ; 1024 786 1025 787 MSState statetab = mstable_.state() ; 1026 ROTableColumn *sharedCol ; 1027 sharedCol = tpool->construct( statetab, "OBS_MODE" ) ; 1028 String obsMode = sharedCol->asString( stateId ) ; 1029 tpool->destroy( sharedCol ) ; 1030 sharedCol = tpool->construct( statetab, "SIG" ) ; 1031 Bool sig = sharedCol->asBool( stateId ) ; 1032 tpool->destroy( sharedCol ) ; 1033 sharedCol = tpool->construct( statetab, "REF" ) ; 1034 Bool ref = sharedCol->asBool( stateId ) ; 1035 tpool->destroy( sharedCol ) ; 1036 sharedCol = tpool->construct( statetab, "CAL" ) ; 1037 Double cal = (Double)(sharedCol->asdouble( stateId )) ; 1038 tpool->destroy( sharedCol ) ; 788 String obsMode = asString( "OBS_MODE", stateId, statetab, tpool ) ; 789 Bool sig = asBool( "SIG", stateId, statetab, tpool ) ; 790 Bool ref = asBool( "REF", stateId, statetab, tpool ) ; 791 Double cal = asDouble( "CAL", stateId, statetab, tpool ) ; 1039 792 //os_ << "OBS_MODE = " << obsMode << LogIO::POST ; 1040 793 … … 1182 935 1183 936 //os_ << "srcType = " << srcType << LogIO::POST ; 1184 //double endSec = gettimeofday_sec() ;1185 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;937 //double endSec = gettimeofday_sec() ; 938 //os_ << "end MSFiller::getSrcType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1186 939 return srcType ; 1187 940 } … … 1190 943 Block<uInt> MSFiller::getPolNo( Int corrType ) 1191 944 { 1192 //double startSec = gettimeofday_sec() ;1193 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ;945 //double startSec = gettimeofday_sec() ; 946 //os_ << "start MSFiller::getPolNo() startSec=" << startSec << LogIO::POST ; 1194 947 Block<uInt> polno( 1 ) ; 1195 948 … … 1221 974 } 1222 975 //os_ << "polno = " << polno << LogIO::POST ; 1223 //double endSec = gettimeofday_sec() ;1224 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;976 //double endSec = gettimeofday_sec() ; 977 //os_ << "end MSFiller::getPolNo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1225 978 1226 979 return polno ; … … 1229 982 String MSFiller::getPolType( Int corrType ) 1230 983 { 1231 //double startSec = gettimeofday_sec() ;1232 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ;984 //double startSec = gettimeofday_sec() ; 985 //os_ << "start MSFiller::getPolType() startSec=" << startSec << LogIO::POST ; 1233 986 String poltype = "" ; 1234 987 … … 1242 995 poltype = "linpol" ; 1243 996 1244 //double endSec = gettimeofday_sec() ;1245 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;997 //double endSec = gettimeofday_sec() ; 998 //os_ << "end MSFiller::getPolType() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1246 999 return poltype ; 1247 1000 } … … 1249 1002 void MSFiller::fillWeather() 1250 1003 { 1251 //double startSec = gettimeofday_sec() ;1252 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ;1004 //double startSec = gettimeofday_sec() ; 1005 //os_ << "start MSFiller::fillWeather() startSec=" << startSec << LogIO::POST ; 1253 1006 1254 1007 if ( !isWeather_ ) { … … 1278 1031 wtab.addRow( wnrow ) ; 1279 1032 1033 Bool stationInfoExists = mWeatherSel.tableDesc().isColumn( "NS_WX_STATION_ID" ) ; 1034 Int stationId = -1 ; 1035 if ( stationInfoExists ) { 1036 // determine which station is closer 1037 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ; 1038 ROArrayColumn<Double> stationPosCol( mWeatherSel, "NS_WX_STATION_POSITION" ) ; 1039 Vector<Int> stationIds = stationCol.getColumn() ; 1040 Vector<Int> stationIdList( 0 ) ; 1041 Matrix<Double> stationPosList( 0, 3, 0.0 ) ; 1042 uInt numStation = 0 ; 1043 for ( uInt i = 0 ; i < stationIds.size() ; i++ ) { 1044 if ( !anyEQ( stationIdList, stationIds[i] ) ) { 1045 numStation++ ; 1046 stationIdList.resize( numStation, True ) ; 1047 stationIdList[numStation-1] = stationIds[i] ; 1048 stationPosList.resize( numStation, 3, True ) ; 1049 stationPosList.row( numStation-1 ) = stationPosCol( i ) ; 1050 } 1051 } 1052 //os_ << "staionIdList = " << stationIdList << endl ; 1053 Table mAntenna = mstable_.antenna() ; 1054 ROArrayColumn<Double> antposCol( mAntenna, "POSITION" ) ; 1055 Vector<Double> antpos = antposCol( antenna_ ) ; 1056 Double minDiff = -1.0 ; 1057 for ( uInt i = 0 ; i < stationIdList.size() ; i++ ) { 1058 Double diff = sum( square( antpos - stationPosList.row( i ) ) ) ; 1059 if ( minDiff < 0.0 || minDiff > diff ) { 1060 minDiff = diff ; 1061 stationId = stationIdList[i] ; 1062 } 1063 } 1064 } 1065 //os_ << "stationId = " << stationId << endl ; 1066 1280 1067 ScalarColumn<Float> *fCol ; 1281 1068 ROScalarColumn<Float> *sharedFloatCol ; … … 1322 1109 ROScalarColumn<Double> tCol( mWeatherSel, "TIME" ) ; 1323 1110 String tUnit = tqCol.getUnits() ; 1324 mwTime_= tCol.getColumn() ;1111 Vector<Double> mwTime = tCol.getColumn() ; 1325 1112 if ( tUnit == "d" ) 1326 mwTime _*= 86400.0 ;1113 mwTime *= 86400.0 ; 1327 1114 tqCol.attach( mWeatherSel, "INTERVAL" ) ; 1328 1115 tCol.attach( mWeatherSel, "INTERVAL" ) ; 1329 1116 String iUnit = tqCol.getUnits() ; 1330 mwInterval_= tCol.getColumn() ;1117 Vector<Double> mwInterval = tCol.getColumn() ; 1331 1118 if ( iUnit == "d" ) 1332 mwInterval_ *= 86400.0 ; 1119 mwInterval *= 86400.0 ; 1120 1121 if ( stationId > 0 ) { 1122 ROScalarColumn<Int> stationCol( mWeatherSel, "NS_WX_STATION_ID" ) ; 1123 Vector<Int> stationVec = stationCol.getColumn() ; 1124 uInt wsnrow = ntrue( stationVec == stationId ) ; 1125 mwTime_.resize( wsnrow ) ; 1126 mwInterval_.resize( wsnrow ) ; 1127 mwIndex_.resize( wsnrow ) ; 1128 uInt wsidx = 0 ; 1129 for ( uInt irow = 0 ; irow < wnrow ; irow++ ) { 1130 if ( stationId == stationVec[irow] ) { 1131 mwTime_[wsidx] = mwTime[irow] ; 1132 mwInterval_[wsidx] = mwInterval[irow] ; 1133 mwIndex_[wsidx] = irow ; 1134 wsidx++ ; 1135 } 1136 } 1137 } 1138 else { 1139 mwTime_ = mwTime ; 1140 mwInterval_ = mwInterval ; 1141 mwIndex_.resize( mwTime_.size() ) ; 1142 indgen( mwIndex_ ) ; 1143 } 1333 1144 //os_ << "mwTime[0] = " << mwTime_[0] << " mwInterval[0] = " << mwInterval_[0] << LogIO::POST ; 1334 //double endSec = gettimeofday_sec() ;1335 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1145 //double endSec = gettimeofday_sec() ; 1146 //os_ << "end MSFiller::fillWeather() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1336 1147 } 1337 1148 1338 1149 void MSFiller::fillFocus() 1339 1150 { 1340 //double startSec = gettimeofday_sec() ;1341 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ;1151 //double startSec = gettimeofday_sec() ; 1152 //os_ << "start MSFiller::fillFocus() startSec=" << startSec << LogIO::POST ; 1342 1153 // tentative 1343 Table tab = table_->focus().table() ; 1344 tab.addRow( 1 ) ; 1345 ScalarColumn<uInt> idCol( tab, "ID" ) ; 1346 idCol.put( 0, 0 ) ; 1347 // double endSec = gettimeofday_sec() ; 1348 // os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1154 table_->focus().addEntry( 0.0, 0.0, 0.0, 0.0 ) ; 1155 //double endSec = gettimeofday_sec() ; 1156 //os_ << "end MSFiller::fillFocus() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1349 1157 } 1350 1158 1351 1159 void MSFiller::fillTcal( boost::object_pool<ROTableColumn> *tpoolr ) 1352 1160 { 1353 //double startSec = gettimeofday_sec() ;1354 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ;1161 //double startSec = gettimeofday_sec() ; 1162 //os_ << "start MSFiller::fillTcal() startSec=" << startSec << LogIO::POST ; 1355 1163 1356 1164 if ( !isSysCal_ ) { … … 1407 1215 Table tab = table_->tcal().table() ; 1408 1216 ArrayColumn<Float> tcalCol( tab, "TCAL" ) ; 1409 ROTableColumn *sharedCol ;1410 1217 uInt oldnr = 0 ; 1411 1218 uInt newnr = 0 ; … … 1418 1225 while( !iter0.pastEnd() ) { 1419 1226 Table t0 = iter0.table() ; 1420 sharedCol = tpoolr->construct( t0, "FEED_ID" ) ; 1421 Int feedId = sharedCol->asInt( 0 ) ; 1422 tpoolr->destroy( sharedCol ) ; 1227 Int feedId = asInt( "FEED_ID", 0, t0, tpoolr ) ; 1423 1228 TableIterator iter1( t0, "SPECTRAL_WINDOW_ID" ) ; 1424 1229 while( !iter1.pastEnd() ) { 1425 1230 Table t1 = iter1.table() ; 1426 sharedCol = tpoolr->construct( t1, "SPECTRAL_WINDOW_ID" ) ; 1427 Int spwId = sharedCol->asInt( 0 ) ; 1428 tpoolr->destroy( sharedCol ) ; 1231 Int spwId = asInt( "SPECTRAL_WINDOW_ID", 0, t1, tpoolr ) ; 1429 1232 tmpTcalCol = new ROArrayColumn<Float>( t1, colTcal_ ) ; 1430 1233 ROScalarQuantColumn<Double> scTimeCol( t1, "TIME" ) ; … … 1460 1263 1461 1264 //tcalrec_.print( std::cout ) ; 1462 //double endSec = gettimeofday_sec() ;1463 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1265 //double endSec = gettimeofday_sec() ; 1266 //os_ << "end MSFiller::fillTcal() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1464 1267 } 1465 1268 1466 1269 uInt MSFiller::getWeatherId( uInt idx, Double wtime ) 1467 1270 { 1468 //double startSec = gettimeofday_sec() ;1469 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ;1271 //double startSec = gettimeofday_sec() ; 1272 //os_ << "start MSFiller::getWeatherId() startSec=" << startSec << LogIO::POST ; 1470 1273 uInt nrow = mwTime_.size() ; 1471 1274 if ( nrow < 2 ) 1472 1275 return 0 ; 1473 1276 uInt wid = nrow ; 1277 if ( idx == 0 ) { 1278 wid = 0 ; 1279 Double tStart = mwTime_[wid]-0.5*mwInterval_[wid] ; 1280 if ( wtime < tStart ) 1281 return wid ; 1282 } 1474 1283 for ( uInt i = idx ; i < nrow-1 ; i++ ) { 1475 1284 Double tStart = mwTime_[i]-0.5*mwInterval_[i] ; … … 1498 1307 //os_ << LogIO::WARN << "Couldn't find correct WEATHER_ID for time " << wtime << LogIO::POST ; 1499 1308 1500 //double endSec = gettimeofday_sec() ;1501 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1309 //double endSec = gettimeofday_sec() ; 1310 //os_ << "end MSFiller::getWeatherId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1502 1311 return wid ; 1503 1312 } … … 1505 1314 void MSFiller::getSysCalTime( Vector<MEpoch> &scTime, Vector<Double> &scInterval, Block<MEpoch> &tcol, Block<Int> &tidx ) 1506 1315 { 1507 //double startSec = gettimeofday_sec() ;1508 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ;1316 //double startSec = gettimeofday_sec() ; 1317 //os_ << "start MSFiller::getSysCalTime() startSec=" << startSec << LogIO::POST ; 1509 1318 1510 1319 if ( !isSysCal_ ) … … 1552 1361 } 1553 1362 } 1554 //double endSec = gettimeofday_sec() ;1555 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ;1363 //double endSec = gettimeofday_sec() ; 1364 //os_ << "end MSFiller::getSysCalTime() endSec=" << endSec << " (" << endSec-startSec << "sec) scnrow = " << scnrow << " tcol.nelements = " << tcol.nelements() << LogIO::POST ; 1556 1365 return ; 1557 1366 } … … 1559 1368 Block<uInt> MSFiller::getTcalId( Int fid, Int spwid, MEpoch &t ) 1560 1369 { 1561 //double startSec = gettimeofday_sec() ;1562 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ;1370 //double startSec = gettimeofday_sec() ; 1371 //os_ << "start MSFiller::getTcalId() startSec=" << startSec << LogIO::POST ; 1563 1372 //if ( table_->tcal().table().nrow() == 0 ) { 1564 1373 if ( !isSysCal_ ) { … … 1583 1392 tcalids[ipol] = ids[0] + ipol - 1 ; 1584 1393 1585 //double endSec = gettimeofday_sec() ;1586 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ;1394 //double endSec = gettimeofday_sec() ; 1395 //os_ << "end MSFiller::getTcalId() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1587 1396 return tcalids ; 1588 1397 } 1589 1398 1590 uInt MSFiller::getDirection( uInt idx, Vector<Double> &dir, Vector<Double> &srate, String &ref, MSPointing &tab, Double t ) 1591 { 1592 // double startSec = gettimeofday_sec() ; 1593 // os_ << "start MSFiller::getDirection() startSec=" << startSec << LogIO::POST ; 1399 uInt MSFiller::getDirection( uInt idx, 1400 Vector<Double> &dir, 1401 Vector<Double> &srate, 1402 String &ref, 1403 ROScalarColumn<Double> &tcol, 1404 ROArrayColumn<Double> &dcol, 1405 Double t ) 1406 { 1407 //double startSec = gettimeofday_sec() ; 1408 //os_ << "start MSFiller::getDirection1() startSec=" << startSec << LogIO::POST ; 1409 //double time0 = gettimeofday_sec() ; 1410 //os_ << "start getDirection 1st stage startSec=" << time0 << LogIO::POST ; 1594 1411 // assume that cols is sorted by TIME 1595 1412 Bool doInterp = False ; 1596 //uInt nrow = cols.nrow() ; 1597 uInt nrow = tab.nrow() ; 1413 uInt nrow = tcol.nrow() ; 1598 1414 if ( nrow == 0 ) 1599 1415 return 0 ; 1600 ROScalarMeasColumn<MEpoch> tcol( tab, "TIME" ) ; 1601 ROArrayMeasColumn<MDirection> dmcol( tab, "DIRECTION" ) ; 1602 ROArrayColumn<Double> dcol( tab, "DIRECTION" ) ; 1416 TableRecord trec = tcol.keywordSet() ; 1417 String tUnit = trec.asArrayString( "QuantumUnits" ).data()[0] ; 1418 Double factor = 1.0 ; 1419 if ( tUnit == "d" ) 1420 factor = 86400.0 ; 1603 1421 // binary search if idx == 0 1604 1422 if ( idx == 0 ) { 1605 uInt nrowb = 75000 ;1423 uInt nrowb = 1000 ; 1606 1424 if ( nrow > nrowb ) { 1607 1425 uInt nblock = nrow / nrowb + 1 ; … … 1609 1427 uInt high = min( nrowb, nrow-iblock*nrowb ) ; 1610 1428 1611 if ( tcol( high-1 ) .get( "s" ).getValue()< t ) {1429 if ( tcol( high-1 ) * factor < t ) { 1612 1430 idx = iblock * nrowb ; 1613 1431 continue ; 1614 1432 } 1615 1433 1616 Vector<MEpoch> tarr( high) ;1617 for ( uInt irow = 0 ; irow < high ; irow++ ) {1618 tarr[irow] = tcol( iblock*nrowb+irow ) ;1619 }1434 RefRows refrows( iblock*nrowb, iblock*nrowb+high, 1 ) ; 1435 Vector<Double> tarr = tcol.getColumnCells( refrows ) ; 1436 if ( tUnit == "d" ) 1437 tarr *= factor ; 1620 1438 1621 1439 uInt bidx = binarySearch( tarr, t ) ; … … 1626 1444 } 1627 1445 else { 1628 Vector<MEpoch> tarr( nrow ) ; 1629 for ( uInt irow = 0 ; irow < nrow ; irow++ ) { 1630 tarr[irow] = tcol( irow ) ; 1631 } 1446 Vector<Double> tarr = tcol.getColumn() ; 1447 if ( tUnit == "d" ) 1448 tarr *= factor ; 1632 1449 idx = binarySearch( tarr, t ) ; 1633 1450 } 1634 1451 } 1452 //double time1 = gettimeofday_sec() ; 1453 //os_ << "end getDirection 1st stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 1635 1454 // ensure that tcol(idx) < t 1636 1455 //os_ << "tcol(idx) = " << tcol(idx).get("s").getValue() << " t = " << t << " diff = " << tcol(idx).get("s").getValue()-t << endl ; 1637 while ( tcol(idx).get("s").getValue() > t && idx > 0 ) 1456 //time0 = gettimeofday_sec() ; 1457 //os_ << "start getDirection 2nd stage startSec=" << time0 << LogIO::POST ; 1458 while( tcol( idx ) * factor > t && idx > 0 ) 1638 1459 idx-- ; 1639 1460 //os_ << "idx = " << idx << LogIO::POST ; … … 1641 1462 // index search 1642 1463 for ( uInt i = idx ; i < nrow ; i++ ) { 1643 Double tref = tcol( i ).get( "s" ).getValue() ; 1464 //Double tref = tcol( i ).get( "s" ).getValue() ; 1465 Double tref = tcol( i ) * factor ; 1644 1466 if ( tref == t ) { 1645 1467 idx = i ; … … 1660 1482 } 1661 1483 } 1484 //time1 = gettimeofday_sec() ; 1485 //os_ << "end getDirection 2nd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 1662 1486 //os_ << "searched idx = " << idx << LogIO::POST ; 1663 1487 1488 //time0 = gettimeofday_sec() ; 1489 //os_ << "start getDirection 3rd stage startSec=" << time0 << LogIO::POST ; 1664 1490 //os_ << "dmcol(idx).shape() = " << dmcol(idx).shape() << LogIO::POST ; 1665 IPosition ip( dmcol(idx).shape().nelements(), 0 ) ; 1491 //IPosition ip( dmcol(idx).shape().nelements(), 0 ) ; 1492 IPosition ip( dcol(idx).shape().nelements(), 0 ) ; 1666 1493 //os_ << "ip = " << ip << LogIO::POST ; 1667 ref = dmcol(idx)(ip).getRefString() ; 1494 //ref = dmcol(idx)(ip).getRefString() ; 1495 trec = dcol.keywordSet() ; 1496 Record rec = trec.asRecord( "MEASINFO" ) ; 1497 ref = rec.asString( "Ref" ) ; 1668 1498 //os_ << "ref = " << ref << LogIO::POST ; 1669 1499 if ( doInterp ) { 1670 1500 //os_ << "do interpolation" << LogIO::POST ; 1671 1501 //os_ << "dcol(idx).shape() = " << dcol(idx).shape() << LogIO::POST ; 1672 Double tref0 = tcol(idx).get("s").getValue() ; 1673 Double tref1 = tcol(idx+1).get("s").getValue() ; 1502 // Double tref0 = tcol(idx).get("s").getValue() ; 1503 // Double tref1 = tcol(idx+1).get("s").getValue() ; 1504 Double tref0 = tcol(idx) * factor ; 1505 Double tref1 = tcol(idx+1) * factor ; 1674 1506 Matrix<Double> mdir0 = dcol( idx ) ; 1675 1507 Matrix<Double> mdir1 = dcol( idx+1 ) ; … … 1697 1529 } 1698 1530 1699 // double endSec = gettimeofday_sec() ; 1700 // os_ << "end MSFiller::getDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1531 //time1 = gettimeofday_sec() ; 1532 //os_ << "end getDirection 3rd stage endSec=" << time1 << " (" << time1-time0 << "sec)" << LogIO::POST ; 1533 //double endSec = gettimeofday_sec() ; 1534 //os_ << "end MSFiller::getDirection1() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1701 1535 return idx ; 1702 1536 } … … 1722 1556 else if ( t > target ) 1723 1557 high = idx - 1 ; 1724 else 1558 else { 1725 1559 return idx ; 1560 } 1726 1561 } 1727 1562 … … 1729 1564 1730 1565 return idx ; 1566 } 1567 1568 uInt MSFiller::binarySearch( Vector<Double> &timeList, Double target ) 1569 { 1570 Int low = 0 ; 1571 Int high = timeList.nelements() ; 1572 uInt idx = 0 ; 1573 1574 while ( low <= high ) { 1575 idx = (Int)( 0.5 * ( low + high ) ) ; 1576 Double t = timeList[idx] ; 1577 if ( t < target ) 1578 low = idx + 1 ; 1579 else if ( t > target ) 1580 high = idx - 1 ; 1581 else { 1582 return idx ; 1583 } 1584 } 1585 1586 idx = max( 0, min( low, high ) ) ; 1587 1588 return idx ; 1589 } 1590 1591 string MSFiller::getFrame() 1592 { 1593 MFrequency::Types frame = MFrequency::DEFAULT ; 1594 ROTableColumn numChanCol( mstable_.spectralWindow(), "NUM_CHAN" ) ; 1595 ROTableColumn measFreqRefCol( mstable_.spectralWindow(), "MEAS_FREQ_REF" ) ; 1596 uInt nrow = numChanCol.nrow() ; 1597 Vector<Int> measFreqRef( nrow, MFrequency::DEFAULT ) ; 1598 uInt nref = 0 ; 1599 for ( uInt irow = 0 ; irow < nrow ; irow++ ) { 1600 if ( numChanCol.asInt( irow ) != 4 ) { // exclude WVR 1601 measFreqRef[nref] = measFreqRefCol.asInt( irow ) ; 1602 nref++ ; 1603 } 1604 } 1605 if ( nref > 0 ) 1606 frame = (MFrequency::Types)measFreqRef[0] ; 1607 1608 return MFrequency::showType( frame ) ; 1609 } 1610 1611 void MSFiller::reshapeSpectraAndFlagtra( Cube<Float> &sp, 1612 Cube<Bool> &fl, 1613 Table &tab, 1614 Int &npol, 1615 Int &nchan, 1616 Int &nrow, 1617 Vector<Int> &corrtype ) 1618 { 1619 //double startSec = gettimeofday_sec() ; 1620 //os_ << "start MSFiller::reshapeSpectraAndFlagtra() startSec=" << startSec << LogIO::POST ; 1621 if ( isFloatData_ ) { 1622 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ; 1623 ROArrayColumn<Float> mFloatDataCol( tab, "FLOAT_DATA" ) ; 1624 sp = mFloatDataCol.getColumn() ; 1625 fl = mFlagCol.getColumn() ; 1626 } 1627 else if ( isData_ ) { 1628 sp.resize( npol, nchan, nrow ) ; 1629 fl.resize( npol, nchan, nrow ) ; 1630 ROArrayColumn<Bool> mFlagCol( tab, "FLAG" ) ; 1631 ROArrayColumn<Complex> mDataCol( tab, "DATA" ) ; 1632 if ( npol < 3 ) { 1633 // Cube<Complex> tmp = mDataCol.getColumn() ; 1634 // Float *sp_p = sp.data() ; 1635 // Complex *tmp_p = tmp.data() ; 1636 // for ( uInt i = 0 ; i < sp.nelements() ; i++ ) 1637 // sp_p[i] = tmp_p[i].real() ; 1638 Cube<Float> tmp = ComplexToReal( mDataCol.getColumn() ) ; 1639 IPosition start( 3, 0, 0, 0 ) ; 1640 IPosition end( 3, 2*npol-1, nchan-1, nrow-1 ) ; 1641 IPosition inc( 3, 2, 1, 1 ) ; 1642 sp = tmp( start, end, inc ) ; 1643 fl = mFlagCol.getColumn() ; 1644 } 1645 else { 1646 for ( Int irow = 0 ; irow < nrow ; irow++ ) { 1647 Bool crossOK = False ; 1648 Matrix<Complex> mSp = mDataCol( irow ) ; 1649 Matrix<Bool> mFl = mFlagCol( irow ) ; 1650 Matrix<Float> spxy = sp.xyPlane( irow ) ; 1651 Matrix<Bool> flxy = fl.xyPlane( irow ) ; 1652 for ( Int ipol = 0 ; ipol < npol ; ipol++ ) { 1653 if ( corrtype[ipol] == Stokes::XY || corrtype[ipol] == Stokes::YX 1654 || corrtype[ipol] == Stokes::RL || corrtype[ipol] == Stokes::LR ) { 1655 if ( !crossOK ) { 1656 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ; 1657 IPosition start( 1, 0 ) ; 1658 IPosition end( 1, 2*nchan-1 ) ; 1659 IPosition inc( 1, 2 ) ; 1660 spxy.row( ipol ) = tmp( start, end, inc ) ; 1661 flxy.row( ipol ) = mFl.row( ipol ) ; 1662 start = IPosition( 1, 1 ) ; 1663 spxy.row( ipol+1 ) = tmp( start, end, inc ) ; 1664 flxy.row( ipol+1 ) = mFl.row( ipol ) ; 1665 if ( corrtype[ipol] == Stokes::YX || corrtype[ipol] == Stokes::LR ) { 1666 spxy.row( ipol+1 ) = spxy.row( ipol+1 ) * (Float)-1.0 ; 1667 } 1668 crossOK = True ; 1669 } 1670 } 1671 else { 1672 Vector<Float> tmp = ComplexToReal( mSp.row( ipol ) ) ; 1673 IPosition start( 1, 0 ) ; 1674 IPosition end( 1, 2*nchan-1 ) ; 1675 IPosition inc( 1, 2 ) ; 1676 spxy.row( ipol ) = tmp( start, end, inc ) ; 1677 flxy.row( ipol ) = mFl.row( ipol ) ; 1678 } 1679 } 1680 } 1681 } 1682 } 1683 //double endSec = gettimeofday_sec() ; 1684 //os_ << "end MSFiller::reshapeSpectraAndFlagtra() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1685 } 1686 1687 uInt MSFiller::getDirection( uInt idx, 1688 Vector<Double> &dir, 1689 Vector<Double> &azel, 1690 Vector<Double> &srate, 1691 ROScalarColumn<Double> &ptcol, 1692 ROArrayColumn<Double> &pdcol, 1693 MEpoch &t, 1694 MPosition &antpos ) 1695 { 1696 //double startSec = gettimeofday_sec() ; 1697 //os_ << "start MSFiller::getDirection2() startSec=" << startSec << LogIO::POST ; 1698 String refString ; 1699 MDirection::Types dirType ; 1700 uInt diridx = getDirection( idx, dir, srate, refString, ptcol, pdcol, t.get("s").getValue() ) ; 1701 MDirection::getType( dirType, refString ) ; 1702 MeasFrame mf( t, antpos ) ; 1703 if ( refString == "J2000" ) { 1704 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 1705 azel = toazel( dir ).getAngle("rad").getValue() ; 1706 } 1707 else if ( refString(0,4) == "AZEL" ) { 1708 azel = dir.copy() ; 1709 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1710 dir = toj2000( dir ).getAngle("rad").getValue() ; 1711 } 1712 else { 1713 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ; 1714 azel = toazel( dir ).getAngle("rad").getValue() ; 1715 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1716 dir = toj2000( dir ).getAngle("rad").getValue() ; 1717 } 1718 if ( srate.size() == 0 ) { 1719 srate.resize( 2 ) ; 1720 srate = 0.0 ; 1721 } 1722 //double endSec = gettimeofday_sec() ; 1723 //os_ << "end MSFiller::getDirection2() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1724 return diridx ; 1725 } 1726 1727 void MSFiller::getSourceDirection( Vector<Double> &dir, 1728 Vector<Double> &azel, 1729 Vector<Double> &srate, 1730 MEpoch &t, 1731 MPosition &antpos, 1732 Vector<MDirection> &srcdir ) 1733 { 1734 //double startSec = gettimeofday_sec() ; 1735 //os_ << "start MSFiller::getSourceDirection() startSec=" << startSec << LogIO::POST ; 1736 Vector<Double> defaultScanrate( 2, 0.0 ) ; 1737 Vector<Double> defaultDir = srcdir[0].getAngle( "rad" ).getValue() ; 1738 if ( srcdir.nelements() > 1 ) 1739 defaultScanrate = srcdir[1].getAngle( "rad" ).getValue() ; 1740 String ref = srcdir[0].getRefString() ; 1741 MDirection::Types dirType ; 1742 MDirection::getType( dirType, ref ) ; 1743 MeasFrame mf( t, antpos ) ; 1744 if ( ref != "J2000" ) { 1745 MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ; 1746 dir = toj2000( defaultDir ).getAngle("rad").getValue() ; 1747 } 1748 else 1749 dir = defaultDir ; 1750 if ( ref != "AZELGEO" ) { 1751 MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ; 1752 azel = toazel( defaultDir ).getAngle("rad").getValue() ; 1753 } 1754 else 1755 azel = defaultDir ; 1756 srate = defaultScanrate ; 1757 //double endSec = gettimeofday_sec() ; 1758 //os_ << "end MSFiller::getSourceDirection() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1759 } 1760 1761 void MSFiller::initHeader( STHeader &header ) 1762 { 1763 header.nchan = 0 ; 1764 header.npol = 0 ; 1765 header.nif = 0 ; 1766 header.nbeam = 0 ; 1767 header.observer = "" ; 1768 header.project = "" ; 1769 header.obstype = "" ; 1770 header.antennaname = "" ; 1771 header.antennaposition.resize( 0 ) ; 1772 header.equinox = 0.0 ; 1773 header.freqref = "" ; 1774 header.reffreq = -1.0 ; 1775 header.bandwidth = 0.0 ; 1776 header.utc = 0.0 ; 1777 header.fluxunit = "" ; 1778 header.epoch = "" ; 1779 header.poltype = "" ; 1780 } 1781 1782 String MSFiller::asString( String name, 1783 uInt idx, 1784 Table tab, 1785 boost::object_pool<ROTableColumn> *pool ) 1786 { 1787 ROTableColumn *col = pool->construct( tab, name ) ; 1788 String v = col->asString( idx ) ; 1789 pool->destroy( col ) ; 1790 return v ; 1791 } 1792 1793 Bool MSFiller::asBool( String name, 1794 uInt idx, 1795 Table &tab, 1796 boost::object_pool<ROTableColumn> *pool ) 1797 { 1798 ROTableColumn *col = pool->construct( tab, name ) ; 1799 Bool v = col->asBool( idx ) ; 1800 pool->destroy( col ) ; 1801 return v ; 1802 } 1803 1804 uInt MSFiller::asuInt( String name, 1805 uInt idx, 1806 Table &tab, 1807 boost::object_pool<ROTableColumn> *pool ) 1808 { 1809 ROTableColumn *col = pool->construct( tab, name ) ; 1810 uInt v = col->asuInt( idx ) ; 1811 pool->destroy( col ) ; 1812 return v ; 1813 } 1814 1815 Int MSFiller::asInt( String name, 1816 uInt idx, 1817 Table &tab, 1818 boost::object_pool<ROTableColumn> *pool ) 1819 { 1820 ROTableColumn *col = pool->construct( tab, name ) ; 1821 Int v = col->asInt( idx ) ; 1822 pool->destroy( col ) ; 1823 return v ; 1824 } 1825 1826 Float MSFiller::asFloat( String name, 1827 uInt idx, 1828 Table &tab, 1829 boost::object_pool<ROTableColumn> *pool ) 1830 { 1831 ROTableColumn *col = pool->construct( tab, name ) ; 1832 Float v = col->asfloat( idx ) ; 1833 pool->destroy( col ) ; 1834 return v ; 1835 } 1836 1837 Double MSFiller::asDouble( String name, 1838 uInt idx, 1839 Table &tab, 1840 boost::object_pool<ROTableColumn> *pool ) 1841 { 1842 ROTableColumn *col = pool->construct( tab, name ) ; 1843 Double v = col->asdouble( idx ) ; 1844 pool->destroy( col ) ; 1845 return v ; 1846 } 1847 1848 void MSFiller::sourceInfo( Int sourceId, 1849 Int spwId, 1850 String &name, 1851 MDirection &direction, 1852 Vector<casa::Double> &properMotion, 1853 Vector<casa::Double> &restFreqs, 1854 Vector<casa::String> &transitions, 1855 Vector<casa::Double> &sysVels, 1856 boost::object_pool<ROTableColumn> *tpoolr ) 1857 { 1858 //double startSec = gettimeofday_sec() ; 1859 //os_ << "start MSFiller::sourceInfo() startSec=" << startSec << LogIO::POST ; 1860 1861 MSSource srctab = mstable_.source() ; 1862 MSSource srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == spwId ) ; 1863 if ( srctabSel.nrow() == 0 ) { 1864 srctabSel = srctab( srctab.col("SOURCE_ID") == sourceId && srctab.col("SPECTRAL_WINDOW_ID") == -1 ) ; 1865 } 1866 Int numLines = 0 ; 1867 if ( srctabSel.nrow() > 0 ) { 1868 // source name 1869 name = asString( "NAME", 0, srctabSel, tpoolr ) ; 1870 1871 // source proper motion 1872 ROArrayColumn<Double> roArrDCol( srctabSel, "PROPER_MOTION" ) ; 1873 properMotion = roArrDCol( 0 ) ; 1874 1875 // source direction as MDirection object 1876 ROScalarMeasColumn<MDirection> tmpMeasCol( srctabSel, "DIRECTION" ) ; 1877 direction = tmpMeasCol( 0 ) ; 1878 1879 // number of lines 1880 numLines = asInt( "NUM_LINES", 0, srctabSel, tpoolr ) ; 1881 } 1882 else { 1883 name = "" ; 1884 properMotion = Vector<Double>( 2, 0.0 ) ; 1885 direction = MDirection( Quantum<Double>(0.0,Unit("rad")), Quantum<Double>(0.0,Unit("rad")) ) ; 1886 } 1887 1888 restFreqs.resize( numLines ) ; 1889 transitions.resize( numLines ) ; 1890 sysVels.resize( numLines ) ; 1891 if ( numLines > 0 ) { 1892 if ( srctabSel.tableDesc().isColumn( "REST_FREQUENCY" ) ) { 1893 ROArrayQuantColumn<Double> quantArrCol( srctabSel, "REST_FREQUENCY" ) ; 1894 Array< Quantum<Double> > qRestFreqs = quantArrCol( 0 ) ; 1895 for ( int i = 0 ; i < numLines ; i++ ) { 1896 restFreqs[i] = qRestFreqs( IPosition( 1, i ) ).getValue( "Hz" ) ; 1897 } 1898 } 1899 //os_ << "restFreqs = " << restFreqs << LogIO::POST ; 1900 if ( srctabSel.tableDesc().isColumn( "TRANSITION" ) ) { 1901 ROArrayColumn<String> transitionCol( srctabSel, "TRANSITION" ) ; 1902 if ( transitionCol.isDefined( 0 ) ) 1903 transitions = transitionCol( 0 ) ; 1904 //os_ << "transitionNameCol.nrow() = " << transitionCol.nrow() << LogIO::POST ; 1905 } 1906 if ( srctabSel.tableDesc().isColumn( "SYSVEL" ) ) { 1907 ROArrayColumn<Double> roArrDCol( srctabSel, "SYSVEL" ) ; 1908 sysVels = roArrDCol( 0 ) ; 1909 } 1910 } 1731 1911 1912 //double endSec = gettimeofday_sec() ; 1913 //os_ << "end MSFiller::sourceInfo() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1914 } 1915 1916 void MSFiller::spectralSetup( Int spwId, 1917 MEpoch &me, 1918 MPosition &mp, 1919 MDirection &md, 1920 Double &refpix, 1921 Double &refval, 1922 Double &increment, 1923 Int &nchan, 1924 String &freqref, 1925 Double &reffreq, 1926 Double &bandwidth, 1927 boost::object_pool<ROTableColumn> *tpoolr ) 1928 { 1929 //double startSec = gettimeofday_sec() ; 1930 //os_ << "start MSFiller::spectralSetup() startSec=" << startSec << LogIO::POST ; 1931 1932 MSSpectralWindow spwtab = mstable_.spectralWindow() ; 1933 MeasFrame mf( me, mp, md ) ; 1934 MFrequency::Types freqRef = MFrequency::castType( (uInt)asInt( "MEAS_FREQ_REF", spwId, spwtab, tpoolr ) ) ; 1935 Bool even = False ; 1936 if ( (nchan/2)*2 == nchan ) even = True ; 1937 ROScalarQuantColumn<Double> tmpQuantCol( spwtab, "TOTAL_BANDWIDTH" ) ; 1938 Double totbw = tmpQuantCol( spwId ).getValue( "Hz" ) ; 1939 if ( nchan != 4 ) 1940 bandwidth = max( bandwidth, totbw ) ; 1941 if ( freqref == "" && nchan != 4) 1942 //sdh.freqref = MFrequency::showType( freqRef ) ; 1943 freqref = "LSRK" ; 1944 if ( reffreq == -1.0 && nchan != 4 ) { 1945 tmpQuantCol.attach( spwtab, "REF_FREQUENCY" ) ; 1946 Quantum<Double> qreffreq = tmpQuantCol( spwId ) ; 1947 if ( freqRef == MFrequency::LSRK ) { 1948 reffreq = qreffreq.getValue("Hz") ; 1949 } 1950 else { 1951 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1952 reffreq = tolsr( qreffreq ).get("Hz").getValue() ; 1953 } 1954 } 1955 Int refchan = nchan / 2 ; 1956 IPosition refip( 1, refchan ) ; 1957 refpix = 0.5*(nchan-1) ; 1958 refval = 0.0 ; 1959 ROArrayQuantColumn<Double> sharedQDArrCol( spwtab, "CHAN_WIDTH" ) ; 1960 increment = sharedQDArrCol( spwId )( refip ).getValue( "Hz" ) ; 1961 // os_ << "nchan = " << nchan << " refchan = " << refchan << "(even=" << even << ") refpix = " << refpix << LogIO::POST ; 1962 sharedQDArrCol.attach( spwtab, "CHAN_FREQ" ) ; 1963 Vector< Quantum<Double> > chanFreqs = sharedQDArrCol( spwId ) ; 1964 if ( nchan > 1 && chanFreqs[0].getValue("Hz") > chanFreqs[1].getValue("Hz") ) 1965 increment *= -1.0 ; 1966 if ( freqRef == MFrequency::LSRK ) { 1967 if ( even ) { 1968 IPosition refip0( 1, refchan-1 ) ; 1969 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1970 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1971 refval = 0.5 * ( refval0 + refval1 ) ; 1972 } 1973 else { 1974 refval = chanFreqs(refip).getValue("Hz") ; 1975 } 1976 } 1977 else { 1978 MFrequency::Convert tolsr( freqRef, MFrequency::Ref( MFrequency::LSRK, mf ) ) ; 1979 if ( even ) { 1980 IPosition refip0( 1, refchan-1 ) ; 1981 Double refval0 = chanFreqs(refip0).getValue("Hz") ; 1982 Double refval1 = chanFreqs(refip).getValue("Hz") ; 1983 refval = 0.5 * ( refval0 + refval1 ) ; 1984 refval = tolsr( refval ).get("Hz").getValue() ; 1985 } 1986 else { 1987 refval = tolsr( chanFreqs(refip) ).get("Hz").getValue() ; 1988 } 1989 } 1990 1991 //double endSec = gettimeofday_sec() ; 1992 //os_ << "end MSFiller::spectralSetup() endSec=" << endSec << " (" << endSec-startSec << "sec)" << LogIO::POST ; 1732 1993 } 1733 1994 -
branches/parallel/src/MSFiller.h
r2021 r2240 90 90 91 91 // get direction for DIRECTION, AZIMUTH, and ELEVATION columns 92 casa::uInt getDirection( casa::uInt idx, casa::Vector<casa::Double> &dir, casa::Vector<casa::Double> &srate, casa::String &ref, casa::MSPointing &tab, casa::Double t ) ; 92 casa::uInt getDirection( casa::uInt idx, 93 casa::Vector<casa::Double> &dir, 94 casa::Vector<casa::Double> &srate, 95 casa::String &ref, 96 casa::ROScalarColumn<casa::Double> &ptcol, 97 casa::ROArrayColumn<casa::Double> &pdcol, 98 casa::Double t ) ; 99 casa::uInt getDirection( casa::uInt idx, 100 casa::Vector<casa::Double> &dir, 101 casa::Vector<casa::Double> &azel, 102 casa::Vector<casa::Double> &srate, 103 casa::ROScalarColumn<casa::Double> &ptcol, 104 casa::ROArrayColumn<casa::Double> &pdcol, 105 casa::MEpoch &t, casa::MPosition &antpos ) ; 106 void getSourceDirection( casa::Vector<casa::Double> &dir, 107 casa::Vector<casa::Double> &azel, 108 casa::Vector<casa::Double> &srate, 109 casa::MEpoch &t, 110 casa::MPosition &antpos, 111 casa::Vector<casa::MDirection> &srcdir ) ; 93 112 94 113 // create key for TCAL table … … 97 116 // binary search 98 117 casa::uInt binarySearch( casa::Vector<casa::MEpoch> &timeList, casa::Double target ) ; 118 casa::uInt binarySearch( casa::Vector<casa::Double> &timeList, casa::Double target ) ; 99 119 100 120 // tool for HPC 101 121 double gettimeofday_sec() ; 102 122 123 // get frequency frame 124 std::string getFrame() ; 125 126 // reshape SPECTRA and FLAGTRA 127 void reshapeSpectraAndFlagtra( casa::Cube<casa::Float> &sp, 128 casa::Cube<casa::Bool> &fl, 129 casa::Table &tab, 130 casa::Int &npol, 131 casa::Int &nchan, 132 casa::Int &nrow, 133 casa::Vector<casa::Int> &corrtype ) ; 134 135 // initialize header 136 void initHeader( STHeader &header ) ; 137 138 // get value from table column using object pool 139 casa::String asString( casa::String name, 140 casa::uInt idx, 141 casa::Table tab, 142 boost::object_pool<casa::ROTableColumn> *pool ) ; 143 casa::Bool asBool( casa::String name, 144 casa::uInt idx, 145 casa::Table &tab, 146 boost::object_pool<casa::ROTableColumn> *pool ) ; 147 casa::Int asInt( casa::String name, 148 casa::uInt idx, 149 casa::Table &tab, 150 boost::object_pool<casa::ROTableColumn> *pool ) ; 151 casa::uInt asuInt( casa::String name, 152 casa::uInt idx, 153 casa::Table &tab, 154 boost::object_pool<casa::ROTableColumn> *pool ) ; 155 casa::Float asFloat( casa::String name, 156 casa::uInt idx, 157 casa::Table &tab, 158 boost::object_pool<casa::ROTableColumn> *pool ) ; 159 casa::Double asDouble( casa::String name, 160 casa::uInt idx, 161 casa::Table &tab, 162 boost::object_pool<casa::ROTableColumn> *pool ) ; 163 164 void sourceInfo( casa::Int sourceId, 165 casa::Int spwId, 166 casa::String &name, 167 casa::MDirection &direction, 168 casa::Vector<casa::Double> &properMotion, 169 casa::Vector<casa::Double> &restFreqs, 170 casa::Vector<casa::String> &transitions, 171 casa::Vector<casa::Double> &sysVels, 172 boost::object_pool<casa::ROTableColumn> *pool ) ; 173 174 void spectralSetup( casa::Int spwId, 175 casa::MEpoch &me, 176 casa::MPosition &mp, 177 casa::MDirection &md, 178 casa::Double &refpix, 179 casa::Double &refval, 180 casa::Double &increment, 181 casa::Int &nchan, 182 casa::String &freqref, 183 casa::Double &reffreq, 184 casa::Double &bandwidth, 185 boost::object_pool<casa::ROTableColumn> *pool ) ; 186 103 187 casa::CountedPtr<Scantable> table_ ; 104 188 casa::MeasurementSet mstable_ ; … … 126 210 casa::Vector<casa::Double> mwTime_ ; 127 211 casa::Vector<casa::Double> mwInterval_ ; 212 casa::Vector<casa::uInt> mwIndex_ ; 128 213 129 214 // Record for TCAL_ID -
branches/parallel/src/SConscript
- Property svn:mergeinfo changed (with no actual effect on merging)
Note:
See TracChangeset
for help on using the changeset viewer.