Changeset 2588


Ignore:
Timestamp:
07/06/12 18:30:00 (12 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: sdsave unit test

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Tuning NROReader.
Direction conversion (getDirection and getSourceDirection) is
optimized.


Location:
branches/hpc34/external-alma/atnf/PKSIO
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/hpc34/external-alma/atnf/PKSIO/NROReader.cc

    r2489 r2588  
    182182//----------------------------------------------------------------------
    183183// constructor
    184 NROReader::NROReader( string name ) {
     184NROReader::NROReader( string name )
     185  : dataset_( NULL ),
     186    srcdir_( 0 ),
     187    msrcdir_( 0 )
     188{
    185189  // initialization
    186190  filename_ = name ;
    187   dataset_ = NULL ;
     191  coord_ = -1 ;
    188192}
    189193
     
    298302  LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
    299303
    300   Vector<Double> v ;
    301   Double srcra = Double( dataset_->getRA0() ) ;
    302   Double srcdec = Double( dataset_->getDEC0() ) ;
     304  if ( msrcdir_.nelements() == 2 )
     305    return msrcdir_ ;
     306
     307  srcdir_.resize( 2 ) ;
     308  srcdir_[0] = Double( dataset_->getRA0() ) ;
     309  srcdir_[1] = Double( dataset_->getDEC0() ) ;
    303310  char epoch[5] ;
    304311  strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
     
    306313    // convert to J2000 value
    307314    MDirection result =
    308       MDirection::Convert( MDirection( Quantity( srcra, "rad" ),
    309                                        Quantity( srcdec, "rad" ),
     315      MDirection::Convert( MDirection( Quantity( srcdir_[0], "rad" ),
     316                                       Quantity( srcdir_[1], "rad" ),
    310317                                       MDirection::Ref( MDirection::B1950 ) ),
    311318                           MDirection::Ref( MDirection::J2000 ) ) () ;
    312     v = result.getAngle().getValue() ;
    313     Double srcra2 = v( 0 ) ;
    314     if ( srcra2 < 0.0 && srcra >= 0.0 )
    315       v( 0 ) = 2.0 * M_PI + srcra2 ;
     319    msrcdir_ = result.getAngle().getValue() ;
     320    if ( msrcdir_[0] < 0.0 && srcdir_[0] >= 0.0 )
     321      msrcdir_[0] = 2.0 * M_PI + msrcdir_[0] ;
    316322    //cout << "NROReader::getSourceDirection()  SRCDIRECTION convert from ("
    317323    //<< srcra << "," << srcdec << ") B1950 to ("
     
    322328  }
    323329  else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
    324     v.resize( 2 ) ;
    325     v( 0 ) = srcra ;
    326     v( 1 ) = srcdec ;
    327   }
    328    
    329   return v ;
     330    msrcdir_.reference( srcdir_ ) ;
     331  }
     332   
     333  return msrcdir_ ;
    330334}
    331335
     
    335339  LogIO os( LogOrigin( "NROReader", "getDirection()", WHERE ) ) ;
    336340
    337   Vector<Double> v ;
     341  Vector<Double> v( 2 ) ;
    338342  NRODataRecord *record = dataset_->getRecord( i ) ;
    339343  char epoch[5] ;
    340344  strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
    341345  int icoord = dataset_->getSCNCD() ;
    342   Double dirx = Double( record->SCX ) ;
    343   Double diry = Double( record->SCY ) ;
    344   if ( icoord == 1 ) {
    345     // convert from LB to RADEC
    346     MDirection result =
    347       MDirection::Convert( MDirection( Quantity( dirx, "rad" ),
    348                                        Quantity( diry, "rad" ),
    349                                        MDirection::Ref( MDirection::GALACTIC ) ),
    350                            MDirection::Ref( MDirection::J2000 ) ) () ;
    351     v = result.getAngle().getValue() ;
    352     Double dirx2 = v( 0 ) ;
    353     if ( dirx2 < 0.0 && dirx >= 0.0 )
    354       v( 0 ) = 2.0 * M_PI + dirx2 ;
    355     //cout << "NROReader::getDirection()  DIRECTION convert from ("
    356     //<< dirx << "," << diry << ") LB to ("
    357     //<< v( 0 ) << ","<< v( 1 ) << ") RADEC" << endl ;
    358     //os << LogIO::NORMAL << "DIRECTION convert from ("
    359     //   << dirx << "," << diry << ") LB to ("
    360     //   << v( 0 ) << ","<< v( 1 ) << ") RADEC" << LogIO::POST ;
    361   }
    362   else if ( icoord == 2 ) {
    363     // convert from AZEL to RADEC
    364     vector<double> antpos = getAntennaPosition() ;
    365     Vector<Quantity> qantpos( 3 ) ;
    366     for ( int ip = 0 ; ip < 3 ; ip++ )
    367       qantpos[ip] = Quantity( antpos[ip], "m" ) ;
    368     Double scantime = Double( dataset_->getScanTime( i ) ) ;
    369     MEpoch me( Quantity( scantime, "d" ), MEpoch::UTC ) ;
    370     MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
    371     MeasFrame mf( me, mp ) ;
    372     MDirection result =
    373       MDirection::Convert( MDirection( Quantity( dirx, "rad" ),
    374                                        Quantity( diry, "rad" ),
    375                                        MDirection::Ref( MDirection::AZEL ) ),
    376                            MDirection::Ref( MDirection::J2000, mf ) ) () ;
    377     v = result.getAngle().getValue() ;
    378     //cout << "NROReader::getDirection()  DIRECTION convert from ("
    379     //<< dirx << "," << diry << ") AZEL to ("
    380     //<< v( 0 ) << ","<< v( 1 ) << ") RADEC" << endl ;
    381     //os << LogIO::NORMAL << "DIRECTION convert from ("
    382     //   << dirx << "," << diry << ") AZEL to ("
    383     //   << v( 0 ) << ","<< v( 1 ) << ") RADEC" << LogIO::POST ;
    384   }
    385   else if ( icoord == 0 ) {
    386     if ( strncmp( epoch, "B1950", 5 ) == 0 ) {
    387       // convert to J2000 value
    388       MDirection result =
    389         MDirection::Convert( MDirection( Quantity( dirx, "rad" ),
    390                                          Quantity( diry, "rad" ),
    391                                          MDirection::Ref( MDirection::B1950 ) ),
    392                              MDirection::Ref( MDirection::J2000 ) ) () ;
    393       v = result.getAngle().getValue() ;
    394       Double dirx2 = v( 0 ) ;
    395       if ( dirx2 < 0.0 && dirx >= 0.0 )
    396         v( 0 ) = 2.0 * M_PI + dirx2 ;
    397       //cout << "STFiller::readNRO()  DIRECTION convert from ("
    398       //<< dirx << "," << diry << ") B1950 to ("
    399       //<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
    400       //os << LogIO::NORMAL << "DIRECTION convert from ("
    401       //   << dirx << "," << diry << ") B1950 to ("
    402       //   << v( 0 ) << ","<< v( 1 ) << ") J2000" << LogIO::POST ;
    403     }
    404     else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
    405       v.resize( 2 ) ;
    406       v( 0 ) = dirx ;
    407       v( 1 ) = diry ;
    408     }
    409   }
    410 
    411   return v ;
     346  double scantime = dataset_->getScanTime( i ) ;
     347  initConvert( icoord, scantime, epoch ) ;
     348  v[0]= Double( record->SCX ) ;
     349  v[1] = Double( record->SCY ) ;
     350  if ( converter_.null() )
     351    // no conversion
     352    return v ;
     353  else {
     354    Vector<Double> v2 = (*converter_)( v ).getAngle().getValue() ;
     355    if ( v2[0] < 0.0 && v[0] >= 0.0 )
     356      v2[0] = 2.0 * M_PI + v2[0] ;
     357    return v2 ;
     358  }
     359}
     360
     361void NROReader::initConvert( int icoord, double t, char *epoch )
     362{
     363  if ( icoord == 0 && strncmp( epoch, "J2000", 5 ) == 0 )
     364    // no conversion
     365    return ;
     366
     367  if ( converter_.null() || icoord != coord_ ) {
     368    LogIO os( LogOrigin( "NROReader", "initConvert()", WHERE ) ) ;
     369    coord_ = icoord ;
     370    if ( coord_ == 0 ) {
     371      // RADEC (B1950) -> RADEC (J2000)
     372      os << "Creating converter from RADEC (B1950) to RADEC (J2000)" << LogIO::POST ;
     373      converter_ = new MDirection::Convert( MDirection::B1950,
     374                                            MDirection::J2000 ) ;
     375    }
     376    else if ( coord_ == 1 ) {
     377      // LB -> RADEC (J2000)
     378      os << "Creating converter from GALACTIC to RADEC (J2000)" << LogIO::POST ;
     379      converter_ = new MDirection::Convert( MDirection::GALACTIC,
     380                                            MDirection::J2000 ) ;
     381    }
     382    else {
     383      // coord_ == 2
     384      // AZEL -> RADEC (J2000)
     385      os << "Creating converter from AZEL to RADEC (J2000)" << LogIO::POST ;
     386      if ( mf_.null() ) {
     387        mf_ = new MeasFrame() ;
     388        vector<double> antpos = getAntennaPosition() ;
     389        Vector<Quantity> qantpos( 3 ) ;
     390        for ( int ip = 0 ; ip < 3 ; ip++ )
     391          qantpos[ip] = Quantity( antpos[ip], "m" ) ;
     392        MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
     393        mf_->set( mp ) ;
     394      }
     395      converter_ = new MDirection::Convert( MDirection::AZEL,
     396                                            MDirection::Ref(MDirection::J2000,
     397                                                            *mf_ ) ) ;
     398    }
     399  }
     400
     401  if ( coord_ == 2 ) {
     402    MEpoch me( Quantity( t, "d" ), MEpoch::UTC ) ;
     403    mf_->set( me ) ;
     404  }
    412405}
    413406
  • branches/hpc34/external-alma/atnf/PKSIO/NROReader.h

    r2289 r2588  
    4848#include <measures/Measures/MeasFrame.h>
    4949#include <casa/Logging/LogIO.h>
     50#include <casa/Utilities/CountedPtr.h>
    5051
    5152//#include <fitsio.h>
     
    206207  // Get DIRECTION in RADEC(J2000)
    207208  virtual Vector<Double> getDirection( int i ) ;
     209  virtual void initConvert( int icoord, double t, char *epoch ) ;
    208210
    209211  // filename
     
    212214  // dataset
    213215  NRODataset *dataset_ ;
     216
     217  // source direction
     218  Vector<Double> srcdir_ ;
     219  Vector<Double> msrcdir_ ;
     220
     221  // for direction conversion
     222  CountedPtr<MDirection::Convert> converter_ ;
     223  CountedPtr<MeasFrame> mf_ ;
     224  int coord_ ;
    214225
    215226  // Logger
Note: See TracChangeset for help on using the changeset viewer.