Ignore:
Timestamp:
08/24/12 10:52:04 (12 years ago)
Author:
ShinnosukeKawakami
Message:

hpc34 to trunk 24th Aug. 2012

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/external-alma

  • trunk/external-alma/atnf/PKSIO/NRODataset.cc

    r2434 r2643  
    5858NRODataset::NRODataset( string name )
    5959{
    60   LogIO os( LogOrigin( "NRODataset", "NRODataset()", WHERE ) ) ;
    61 
    6260  // memory allocation
    6361  initialize() ;
    … …  
    243241NRODataRecord *NRODataset::getRecord( int i )
    244242{
    245   LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
    246 
    247243  // DEBUG
    248244  //cout << "NRODataset::getRecord()  Start " << i << endl ;
    249245  //
    250246  if ( i < 0 || i >= rowNum_ ) {
     247    LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
    251248    //cerr << "NRODataset::getRecord()  data index out of range." << endl ;
    252249    os << LogIO::SEVERE << "data index " << i << " out of range. return NULL." << LogIO::POST ;
    … …  
    267264  }
    268265  else {
     266    LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
    269267    //cerr << "NRODataset::getRecord()  error while reading data " << i << endl ;
    270268    os << LogIO::SEVERE << "error while reading data " << i << ". return NULL." << LogIO::POST ;
    … …  
    278276int NRODataset::fillRecord( int i )
    279277{
    280   LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
    281 
    282278  int status = 0 ;
    283279
    … …  
    295291  if ( (int)fread( record_, 1, SCAN_HEADER_SIZE, fp_ ) != SCAN_HEADER_SIZE ) {
    296292    //cerr << "Failed to read scan header: " << i << endl ;
     293    LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
    297294    os << LogIO::SEVERE << "Failed to read scan header for " << i << "th row." << LogIO::POST ;
    298295    return -1 ;
    … …  
    300297  if ( (int)fread( record_->LDATA, 1, dataLen_, fp_ ) != dataLen_ ) {
    301298    //cerr << "Failed to read spectral data: " << i << endl ;
     299    LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
    302300    os << LogIO::SEVERE << "Failed to read spectral data for " << i << "th row." << LogIO::POST ;
    303301    return -1 ;
    … …  
    354352vector<double> NRODataset::getSpectrum( int i )
    355353{
    356   LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
    357  
    358354  // DEBUG
    359355  //cout << "NRODataset::getSpectrum() start process (" << i << ")" << endl ;
    … …  
    361357  // size of spectrum is not chmax_ but dataset_->getNCH() after binding
    362358  const int nchan = NUMCH ;
    363   vector<double> bspec( nchan, 0.0 ) ;  // spectrum "after" binding
     359  vector<double> spec( chmax_ ) ;  // spectrum "before" binding
    364360  // DEBUG
    365361  //cout << "NRODataset::getSpectrum()  nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
    … …  
    379375  if ( ( scale == 0.0 ) && ( offset == 0.0 ) ) {
    380376    //cerr << "NRODataset::getSpectrum()  zero spectrum (" << i << ")" << endl ;
    381     return bspec ;
     377    LogIO os( LogOrigin("NRODataset","getSpectrum",WHERE) ) ;
     378    os << LogIO::WARN << "zero spectrum for row " << i << LogIO::POST ;
     379    if ( spec.size() != nchan )
     380      spec.resize( nchan ) ;
     381    for ( vector<double>::iterator i = spec.begin() ;
     382          i != spec.end() ; i++ )
     383      *i = 0.0 ;
     384    return spec ;
    382385  }
    383386  unsigned char *cdata = (unsigned char *)record->LDATA ;
    … …  
    387390  int chmin = CHMIN ;
    388391
    389   // char -> int
    390   vector<int> ispec( chmax_, 0 ) ;
     392  // char -> int -> double
     393  vector<double>::iterator iter = spec.begin() ;
    391394
    392395  static const int shift_right[] = {
    … …  
    401404  int j = 0 ;
    402405  for ( int i = 0 ; i < chmax_ ; i++ ) {
     406    // char -> int
    403407    int ivalue = 0 ;
    404408    if ( bit == 12 ) {  // 12 bit qunatization
    405       const int idx = j + start_pos[i & 1];
     409      const int ialt = i & 1 ;
     410      const int idx = j + start_pos[ialt];
    406411      const unsigned tmp = unsigned(cdata[idx]) << 8 | cdata[idx + 1];
    407       ivalue = int((tmp >> shift_right[i & 1]) & 0xFFF);
    408       j += incr[i & 1];
    409     }
    410 
    411     ispec[i] = ivalue ;
    412     if ( ( ispec[i] < 0 ) || ( ispec[i] > 4096 ) ) {
     412      ivalue = int((tmp >> shift_right[ialt]) & 0xFFF);
     413      j += incr[ialt];
     414    }
     415
     416    if ( ( ivalue < 0 ) || ( ivalue > 4096 ) ) {
    413417      //cerr << "NRODataset::getSpectrum()  ispec[" << i << "] is out of range" << endl ;
    414       os << LogIO::SEVERE << "ispec[" << i << "] is out of range" << LogIO::EXCEPTION ;
    415       return bspec ;
     418      LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
     419      os << LogIO::SEVERE << "ivalue for row " << i << " is out of range" << LogIO::EXCEPTION ;
     420      if ( spec.size() != nchan )
     421        spec.resize( nchan ) ;
     422      for ( vector<double>::iterator i = spec.begin() ;
     423            i != spec.end() ; i++ )
     424        *i = 0.0 ;
     425      return spec ;
    416426    }
    417427    // DEBUG
    418428    //cout << "NRODataset::getSpectrum()  ispec[" << i << "] = " << ispec[i] << endl ;
    419429    //
    420   }
    421 
    422   double *const spec = new double[ chmax_ ] ;  // spectrum "before" binding
    423   // int -> double
    424   for ( int i = 0 ; i < chmax_ ; i++ ) {
    425     spec[i] = (double)( ispec[i] * scale + offset ) * dscale ;
     430
     431    // int -> double
     432    *iter = (double)( ivalue * scale + offset ) * dscale ;
    426433    // DEBUG
    427     //cout << "NRODataset::getSpectrum()  spec[" << i << "] = " << spec[i] << endl ;
     434    //cout << "NRODataset::getSpectrum()  spec[" << i << "] = " << *iter << endl ;
    428435    //
    429   }
    430 
    431   // channel binding
     436    iter++ ;
     437  }
     438
     439  // channel binding if necessary
    432440  if ( cbind != 1 ) {
    433     int k = chmin ;
    434     double sum0 = 0 ;
    435     double sum1 = 0 ;
     441    iter = spec.begin() ;
     442    advance( iter, chmin ) ;
     443    vector<double>::iterator iter2 = spec.begin() ;
    436444    for ( int i = 0 ; i < nchan ; i++ ) {
    437       for ( int j = k ; j < k + cbind ; j++ ) {
    438         sum0 += spec[k] ;
    439         sum1++ ;
     445      double sum0 = 0 ;
     446      double sum1 = 0 ;
     447      for ( int j = 0 ; j < cbind ; j++ ) {
     448        sum0 += *iter ;
     449        sum1 += 1.0 ;
     450        iter++ ;
    440451      }
    441       bspec[i] = sum0 / sum1 ;
    442       k += cbind ;
     452      *iter2 = sum0 / sum1 ;
     453      iter2++ ;
    443454      // DEBUG
    444455      //cout << "NRODataset::getSpectrum()  bspec[" << i << "] = " << bspec[i] << endl ;
    445456      //
    446457    }
    447   }
    448   else {
    449     for ( int i = 0 ; i < nchan ; i++ )
    450       bspec[i] = spec[i] ;
    451   }
    452   delete[] spec;
     458    spec.resize( nchan ) ;
     459  }
    453460
    454461  // DEBUG
    … …  
    456463  //
    457464
    458   return bspec ;
     465  return spec ;
    459466}
    460467

  • trunk/external-alma/atnf/PKSIO/NROReader.cc

    r2489 r2643  
    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
    … …  
    296300Vector<Double> NROReader::getSourceDirection()
    297301{
    298   LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
    299 
    300   Vector<Double> v ;
    301   Double srcra = Double( dataset_->getRA0() ) ;
    302   Double srcdec = Double( dataset_->getDEC0() ) ;
     302  if ( msrcdir_.nelements() == 2 )
     303    return msrcdir_ ;
     304
     305  srcdir_.resize( 2 ) ;
     306  srcdir_[0] = Double( dataset_->getRA0() ) ;
     307  srcdir_[1] = Double( dataset_->getDEC0() ) ;
    303308  char epoch[5] ;
    304309  strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
    … …  
    306311    // convert to J2000 value
    307312    MDirection result =
    308       MDirection::Convert( MDirection( Quantity( srcra, "rad" ),
    309                                        Quantity( srcdec, "rad" ),
     313      MDirection::Convert( MDirection( Quantity( srcdir_[0], "rad" ),
     314                                       Quantity( srcdir_[1], "rad" ),
    310315                                       MDirection::Ref( MDirection::B1950 ) ),
    311316                           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 ;
     317    msrcdir_ = result.getAngle().getValue() ;
     318    if ( msrcdir_[0] < 0.0 && srcdir_[0] >= 0.0 )
     319      msrcdir_[0] = 2.0 * M_PI + msrcdir_[0] ;
    316320    //cout << "NROReader::getSourceDirection()  SRCDIRECTION convert from ("
    317321    //<< srcra << "," << srcdec << ") B1950 to ("
    318322    //<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
     323    //LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
    319324    //os << LogIO::NORMAL << "SRCDIRECTION convert from ("
    320325    //   << srcra << "," << srcdec << ") B1950 to ("
    … …  
    322327  }
    323328  else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
    324     v.resize( 2 ) ;
    325     v( 0 ) = srcra ;
    326     v( 1 ) = srcdec ;
    327   }
    328    
    329   return v ;
     329    msrcdir_.reference( srcdir_ ) ;
     330  }
     331   
     332  return msrcdir_ ;
    330333}
    331334
    … …  
    333336Vector<Double> NROReader::getDirection( int i )
    334337{
    335   LogIO os( LogOrigin( "NROReader", "getDirection()", WHERE ) ) ;
    336 
    337   Vector<Double> v ;
     338  Vector<Double> v( 2 ) ;
    338339  NRODataRecord *record = dataset_->getRecord( i ) ;
    339340  char epoch[5] ;
    340341  strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
    341342  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 ;
     343  double scantime = dataset_->getScanTime( i ) ;
     344  initConvert( icoord, scantime, epoch ) ;
     345  v[0]= Double( record->SCX ) ;
     346  v[1] = Double( record->SCY ) ;
     347  if ( converter_.null() )
     348    // no conversion
     349    return v ;
     350  else {
     351    Vector<Double> v2 = (*converter_)( v ).getAngle().getValue() ;
     352    if ( v2[0] < 0.0 && v[0] >= 0.0 )
     353      v2[0] = 2.0 * M_PI + v2[0] ;
     354    return v2 ;
     355  }
     356}
     357
     358void NROReader::initConvert( int icoord, double t, char *epoch )
     359{
     360  if ( icoord == 0 && strncmp( epoch, "J2000", 5 ) == 0 )
     361    // no conversion
     362    return ;
     363
     364  if ( converter_.null() || icoord != coord_ ) {
     365    LogIO os( LogOrigin( "NROReader", "initConvert()", WHERE ) ) ;
     366    coord_ = icoord ;
     367    if ( coord_ == 0 ) {
     368      // RADEC (B1950) -> RADEC (J2000)
     369      os << "Creating converter from RADEC (B1950) to RADEC (J2000)" << LogIO::POST ;
     370      converter_ = new MDirection::Convert( MDirection::B1950,
     371                                            MDirection::J2000 ) ;
     372    }
     373    else if ( coord_ == 1 ) {
     374      // LB -> RADEC (J2000)
     375      os << "Creating converter from GALACTIC to RADEC (J2000)" << LogIO::POST ;
     376      converter_ = new MDirection::Convert( MDirection::GALACTIC,
     377                                            MDirection::J2000 ) ;
     378    }
     379    else {
     380      // coord_ == 2
     381      // AZEL -> RADEC (J2000)
     382      os << "Creating converter from AZEL to RADEC (J2000)" << LogIO::POST ;
     383      if ( mf_.null() ) {
     384        mf_ = new MeasFrame() ;
     385        vector<double> antpos = getAntennaPosition() ;
     386        Vector<Quantity> qantpos( 3 ) ;
     387        for ( int ip = 0 ; ip < 3 ; ip++ )
     388          qantpos[ip] = Quantity( antpos[ip], "m" ) ;
     389        MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
     390        mf_->set( mp ) ;
     391      }
     392      converter_ = new MDirection::Convert( MDirection::AZEL,
     393                                            MDirection::Ref(MDirection::J2000,
     394                                                            *mf_ ) ) ;
     395    }
     396  }
     397
     398  if ( coord_ == 2 ) {
     399    MEpoch me( Quantity( t, "d" ), MEpoch::UTC ) ;
     400    mf_->set( me ) ;
     401  }
    412402}
    413403
    … …  
    620610  spectra.resize( spec.size() ) ;
    621611  int index = 0 ;
    622   for ( Vector<Float>::iterator itr = spectra.begin() ; itr != spectra.end() ; itr++ ) {
    623     *itr = spec[index++] ;
    624   }
     612  Bool b ;
     613  Float *fp = spectra.getStorage( b ) ;
     614  Float *wp = fp ;
     615  for ( vector<double>::iterator i = spec.begin() ;
     616        i != spec.end() ; i++ ) {
     617    *wp = *i ;
     618    wp++ ;
     619  }
     620  spectra.putStorage( fp, b ) ;
    625621  //cout << "spec.size() = " << spec.size() << endl ;
    626622 
    627623  // flagtra
    628624  bool setValue = !( flagtra.nelements() == spectra.nelements() ) ;
    629   flagtra.resize( spectra.nelements() ) ;
    630625  if ( setValue ) {
    631626    //cout << "flagtra resized. reset values..." << endl ;
     627    flagtra.resize( spectra.nelements() ) ;
    632628    flagtra.set( 0 ) ;
    633629  }

  • trunk/external-alma/atnf/PKSIO/NROReader.h

    r2289 r2643  
    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.