Changeset 2643 for trunk/external-alma/atnf

Timestamp:

08/24/12 10:52:04 (12 years ago)

Author:

ShinnosukeKawakami

Message:

hpc34 to trunk 24th Aug. 2012

Location:

trunk

Files:

• . (modified) (1 prop)
• external-alma (modified) (1 prop)
• external-alma/atnf/PKSIO/NRODataset.cc (modified) (12 diffs)
• external-alma/atnf/PKSIO/NROReader.cc (modified) (6 diffs)
• external-alma/atnf/PKSIO/NROReader.h (modified) (3 diffs)
• external-alma/atnf/pks/pks_maths.cc (modified) (1 prop)

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

@@ -58 +58 @@
 NRODataset::NRODataset( string name ) 
 {
-  LogIO os( LogOrigin( "NRODataset", "NRODataset()", WHERE ) ) ;
-
   // memory allocation
   initialize() ;
@@ -243 +241 @@
 NRODataRecord *NRODataset::getRecord( int i )
 {
-  LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
-
   // DEBUG
   //cout << "NRODataset::getRecord()  Start " << i << endl ;
   //
   if ( i < 0 || i >= rowNum_ ) {
+    LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
     //cerr << "NRODataset::getRecord()  data index out of range." << endl ;
     os << LogIO::SEVERE << "data index " << i << " out of range. return NULL." << LogIO::POST ;
@@ -267 +264 @@
   }
   else {
+    LogIO os( LogOrigin( "NRODataset", "getRecord()", WHERE ) ) ;
     //cerr << "NRODataset::getRecord()  error while reading data " << i << endl ;
     os << LogIO::SEVERE << "error while reading data " << i << ". return NULL." << LogIO::POST ;
@@ -278 +276 @@
 int NRODataset::fillRecord( int i ) 
 {
-  LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
-
   int status = 0 ;
 
@@ -295 +291 @@
   if ( (int)fread( record_, 1, SCAN_HEADER_SIZE, fp_ ) != SCAN_HEADER_SIZE ) {
     //cerr << "Failed to read scan header: " << i << endl ;
+    LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
     os << LogIO::SEVERE << "Failed to read scan header for " << i << "th row." << LogIO::POST ;
     return -1 ;
@@ -300 +297 @@
   if ( (int)fread( record_->LDATA, 1, dataLen_, fp_ ) != dataLen_ ) {
     //cerr << "Failed to read spectral data: " << i << endl ;
+    LogIO os( LogOrigin( "NRODataset", "fillRecord()", WHERE ) ) ;
     os << LogIO::SEVERE << "Failed to read spectral data for " << i << "th row." << LogIO::POST ;
     return -1 ;
@@ -354 +352 @@
 vector<double> NRODataset::getSpectrum( int i )
 {
-  LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
-  
   // DEBUG
   //cout << "NRODataset::getSpectrum() start process (" << i << ")" << endl ;
@@ -361 +357 @@
   // size of spectrum is not chmax_ but dataset_->getNCH() after binding
   const int nchan = NUMCH ;
-  vector<double> bspec( nchan, 0.0 ) ;  // spectrum "after" binding
+  vector<double> spec( chmax_ ) ;  // spectrum "before" binding
   // DEBUG
   //cout << "NRODataset::getSpectrum()  nchan = " << nchan << " chmax_ = " << chmax_ << endl ;
@@ -379 +375 @@
   if ( ( scale == 0.0 ) && ( offset == 0.0 ) ) {
     //cerr << "NRODataset::getSpectrum()  zero spectrum (" << i << ")" << endl ;
-    return bspec ;
+    LogIO os( LogOrigin("NRODataset","getSpectrum",WHERE) ) ;
+    os << LogIO::WARN << "zero spectrum for row " << i << LogIO::POST ;
+    if ( spec.size() != nchan )
+      spec.resize( nchan ) ;
+    for ( vector<double>::iterator i = spec.begin() ;
+          i != spec.end() ; i++ )
+      *i = 0.0 ;
+    return spec ;
   }
   unsigned char *cdata = (unsigned char *)record->LDATA ;
@@ -387 +390 @@
   int chmin = CHMIN ;
 
-  // char -> int
-  vector<int> ispec( chmax_, 0 ) ;
+  // char -> int -> double
+  vector<double>::iterator iter = spec.begin() ;
 
   static const int shift_right[] = {
@@ -401 +404 @@
   int j = 0 ;
   for ( int i = 0 ; i < chmax_ ; i++ ) {
+    // char -> int
     int ivalue = 0 ;
     if ( bit == 12 ) {  // 12 bit qunatization
-      const int idx = j + start_pos[i & 1];
+      const int ialt = i & 1 ;
+      const int idx = j + start_pos[ialt];
       const unsigned tmp = unsigned(cdata[idx]) << 8 | cdata[idx + 1];
-      ivalue = int((tmp >> shift_right[i & 1]) & 0xFFF);
-      j += incr[i & 1];
-    }
-
-    ispec[i] = ivalue ;
-    if ( ( ispec[i] < 0 ) || ( ispec[i] > 4096 ) ) {
+      ivalue = int((tmp >> shift_right[ialt]) & 0xFFF);
+      j += incr[ialt];
+    }
+
+    if ( ( ivalue < 0 ) || ( ivalue > 4096 ) ) {
       //cerr << "NRODataset::getSpectrum()  ispec[" << i << "] is out of range" << endl ;
-      os << LogIO::SEVERE << "ispec[" << i << "] is out of range" << LogIO::EXCEPTION ;
-      return bspec ;
+      LogIO os( LogOrigin( "NRODataset", "getSpectrum", WHERE ) ) ;
+      os << LogIO::SEVERE << "ivalue for row " << i << " is out of range" << LogIO::EXCEPTION ;
+      if ( spec.size() != nchan )
+        spec.resize( nchan ) ;
+      for ( vector<double>::iterator i = spec.begin() ;
+            i != spec.end() ; i++ )
+        *i = 0.0 ;
+      return spec ;
     }
     // DEBUG
     //cout << "NRODataset::getSpectrum()  ispec[" << i << "] = " << ispec[i] << endl ;
     //
-  }
-
-  double *const spec = new double[ chmax_ ] ;  // spectrum "before" binding
-  // int -> double 
-  for ( int i = 0 ; i < chmax_ ; i++ ) {
-    spec[i] = (double)( ispec[i] * scale + offset ) * dscale ; 
+
+    // int -> double
+    *iter = (double)( ivalue * scale + offset ) * dscale ; 
     // DEBUG
-    //cout << "NRODataset::getSpectrum()  spec[" << i << "] = " << spec[i] << endl ;
+    //cout << "NRODataset::getSpectrum()  spec[" << i << "] = " << *iter << endl ;
     //
-  }
-
-  // channel binding
+    iter++ ;
+  }
+
+  // channel binding if necessary
   if ( cbind != 1 ) {
-    int k = chmin ;
-    double sum0 = 0 ;
-    double sum1 = 0 ;
+    iter = spec.begin() ;
+    advance( iter, chmin ) ;
+    vector<double>::iterator iter2 = spec.begin() ;
     for ( int i = 0 ; i < nchan ; i++ ) {
-      for ( int j = k ; j < k + cbind ; j++ ) {
-        sum0 += spec[k] ;
-        sum1++ ;
+      double sum0 = 0 ;
+      double sum1 = 0 ;
+      for ( int j = 0 ; j < cbind ; j++ ) {
+        sum0 += *iter ;
+        sum1 += 1.0 ;
+        iter++ ;
       }
-      bspec[i] = sum0 / sum1 ;
-      k += cbind ;
+      *iter2 = sum0 / sum1 ;
+      iter2++ ;
       // DEBUG
       //cout << "NRODataset::getSpectrum()  bspec[" << i << "] = " << bspec[i] << endl ;
       //
     }
-  }
-  else {
-    for ( int i = 0 ; i < nchan ; i++ ) 
-      bspec[i] = spec[i] ;
-  }
-  delete[] spec;
+    spec.resize( nchan ) ;
+  }
 
   // DEBUG
@@ -456 +463 @@
   //
 
-  return bspec ;
+  return spec ;
 }
 

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

@@ -182 +182 @@
 //----------------------------------------------------------------------
 // constructor
-NROReader::NROReader( string name ) {
+NROReader::NROReader( string name ) 
+  : dataset_( NULL ),
+    srcdir_( 0 ),
+    msrcdir_( 0 )
+{
   // initialization
   filename_ = name ;
-  dataset_ = NULL ;
+  coord_ = -1 ;
 }
 
@@ -296 +300 @@
 Vector<Double> NROReader::getSourceDirection()
 {
-  LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
-
-  Vector<Double> v ;
-  Double srcra = Double( dataset_->getRA0() ) ;
-  Double srcdec = Double( dataset_->getDEC0() ) ;
+  if ( msrcdir_.nelements() == 2 )
+    return msrcdir_ ;
+
+  srcdir_.resize( 2 ) ;
+  srcdir_[0] = Double( dataset_->getRA0() ) ;
+  srcdir_[1] = Double( dataset_->getDEC0() ) ;
   char epoch[5] ; 
   strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
@@ -306 +311 @@
     // convert to J2000 value
     MDirection result = 
-      MDirection::Convert( MDirection( Quantity( srcra, "rad" ),
-                                       Quantity( srcdec, "rad" ),
+      MDirection::Convert( MDirection( Quantity( srcdir_[0], "rad" ),
+                                       Quantity( srcdir_[1], "rad" ),
                                        MDirection::Ref( MDirection::B1950 ) ),
                            MDirection::Ref( MDirection::J2000 ) ) () ;
-    v = result.getAngle().getValue() ;
-    Double srcra2 = v( 0 ) ;
-    if ( srcra2 < 0.0 && srcra >= 0.0 )
-      v( 0 ) = 2.0 * M_PI + srcra2 ;
+    msrcdir_ = result.getAngle().getValue() ;
+    if ( msrcdir_[0] < 0.0 && srcdir_[0] >= 0.0 )
+      msrcdir_[0] = 2.0 * M_PI + msrcdir_[0] ;
     //cout << "NROReader::getSourceDirection()  SRCDIRECTION convert from (" 
     //<< srcra << "," << srcdec << ") B1950 to (" 
     //<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
+    //LogIO os( LogOrigin( "NROReader", "getSourceDirection()", WHERE ) ) ;
     //os << LogIO::NORMAL << "SRCDIRECTION convert from (" 
     //   << srcra << "," << srcdec << ") B1950 to (" 
@@ -322 +327 @@
   }
   else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
-    v.resize( 2 ) ;
-    v( 0 ) = srcra ;
-    v( 1 ) = srcdec ;
-  }
-    
-  return v ;
+    msrcdir_.reference( srcdir_ ) ;
+  }
+   
+  return msrcdir_ ;
 }
 
@@ -333 +336 @@
 Vector<Double> NROReader::getDirection( int i )
 {
-  LogIO os( LogOrigin( "NROReader", "getDirection()", WHERE ) ) ;
-
-  Vector<Double> v ;
+  Vector<Double> v( 2 ) ;
   NRODataRecord *record = dataset_->getRecord( i ) ;
   char epoch[5] ;
   strncpy( epoch, (dataset_->getEPOCH()).c_str(), 5 ) ;
   int icoord = dataset_->getSCNCD() ;
-  Double dirx = Double( record->SCX ) ;
-  Double diry = Double( record->SCY ) ;
-  if ( icoord == 1 ) {
-    // convert from LB to RADEC
-    MDirection result = 
-      MDirection::Convert( MDirection( Quantity( dirx, "rad" ), 
-                                       Quantity( diry, "rad" ),
-                                       MDirection::Ref( MDirection::GALACTIC ) ),
-                           MDirection::Ref( MDirection::J2000 ) ) () ;
-    v = result.getAngle().getValue() ;
-    Double dirx2 = v( 0 ) ;
-    if ( dirx2 < 0.0 && dirx >= 0.0 ) 
-      v( 0 ) = 2.0 * M_PI + dirx2 ;
-    //cout << "NROReader::getDirection()  DIRECTION convert from (" 
-    //<< dirx << "," << diry << ") LB to (" 
-    //<< v( 0 ) << ","<< v( 1 ) << ") RADEC" << endl ;
-    //os << LogIO::NORMAL << "DIRECTION convert from (" 
-    //   << dirx << "," << diry << ") LB to (" 
-    //   << v( 0 ) << ","<< v( 1 ) << ") RADEC" << LogIO::POST ;
-  }
-  else if ( icoord == 2 ) {
-    // convert from AZEL to RADEC
-    vector<double> antpos = getAntennaPosition() ;
-    Vector<Quantity> qantpos( 3 ) ;
-    for ( int ip = 0 ; ip < 3 ; ip++ )
-      qantpos[ip] = Quantity( antpos[ip], "m" ) ;
-    Double scantime = Double( dataset_->getScanTime( i ) ) ;
-    MEpoch me( Quantity( scantime, "d" ), MEpoch::UTC ) ;
-    MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
-    MeasFrame mf( me, mp ) ;
-    MDirection result = 
-      MDirection::Convert( MDirection( Quantity( dirx, "rad" ), 
-                                       Quantity( diry, "rad" ),
-                                       MDirection::Ref( MDirection::AZEL ) ),
-                           MDirection::Ref( MDirection::J2000, mf ) ) () ;
-    v = result.getAngle().getValue() ;
-    //cout << "NROReader::getDirection()  DIRECTION convert from (" 
-    //<< dirx << "," << diry << ") AZEL to (" 
-    //<< v( 0 ) << ","<< v( 1 ) << ") RADEC" << endl ;
-    //os << LogIO::NORMAL << "DIRECTION convert from (" 
-    //   << dirx << "," << diry << ") AZEL to (" 
-    //   << v( 0 ) << ","<< v( 1 ) << ") RADEC" << LogIO::POST ;
-  }
-  else if ( icoord == 0 ) {
-    if ( strncmp( epoch, "B1950", 5 ) == 0 ) {
-      // convert to J2000 value 
-      MDirection result = 
-        MDirection::Convert( MDirection( Quantity( dirx, "rad" ),
-                                         Quantity( diry, "rad" ),
-                                         MDirection::Ref( MDirection::B1950 ) ),
-                             MDirection::Ref( MDirection::J2000 ) ) () ;
-      v = result.getAngle().getValue() ;
-      Double dirx2 = v( 0 ) ;
-      if ( dirx2 < 0.0 && dirx >= 0.0 )
-        v( 0 ) = 2.0 * M_PI + dirx2 ;
-      //cout << "STFiller::readNRO()  DIRECTION convert from (" 
-      //<< dirx << "," << diry << ") B1950 to (" 
-      //<< v( 0 ) << ","<< v( 1 ) << ") J2000" << endl ;
-      //os << LogIO::NORMAL << "DIRECTION convert from (" 
-      //   << dirx << "," << diry << ") B1950 to (" 
-      //   << v( 0 ) << ","<< v( 1 ) << ") J2000" << LogIO::POST ;
-    }
-    else if ( strncmp( epoch, "J2000", 5 ) == 0 ) {
-      v.resize( 2 ) ;
-      v( 0 ) = dirx ;
-      v( 1 ) = diry ;
-    }
-  }
-
-  return v ;
+  double scantime = dataset_->getScanTime( i ) ;
+  initConvert( icoord, scantime, epoch ) ;
+  v[0]= Double( record->SCX ) ;
+  v[1] = Double( record->SCY ) ;
+  if ( converter_.null() )
+    // no conversion 
+    return v ;
+  else {
+    Vector<Double> v2 = (*converter_)( v ).getAngle().getValue() ; 
+    if ( v2[0] < 0.0 && v[0] >= 0.0 ) 
+      v2[0] = 2.0 * M_PI + v2[0] ;
+    return v2 ;
+  }
+}
+
+void NROReader::initConvert( int icoord, double t, char *epoch )
+{
+  if ( icoord == 0 && strncmp( epoch, "J2000", 5 ) == 0 )
+    // no conversion 
+    return ;
+
+  if ( converter_.null() || icoord != coord_ ) {
+    LogIO os( LogOrigin( "NROReader", "initConvert()", WHERE ) ) ;
+    coord_ = icoord ;
+    if ( coord_ == 0 ) {
+      // RADEC (B1950) -> RADEC (J2000)
+      os << "Creating converter from RADEC (B1950) to RADEC (J2000)" << LogIO::POST ;
+      converter_ = new MDirection::Convert( MDirection::B1950,
+                                            MDirection::J2000 ) ;
+    }
+    else if ( coord_ == 1 ) {
+      // LB -> RADEC (J2000)
+      os << "Creating converter from GALACTIC to RADEC (J2000)" << LogIO::POST ;
+      converter_ = new MDirection::Convert( MDirection::GALACTIC,
+                                            MDirection::J2000 ) ;
+    }
+    else {
+      // coord_ == 2
+      // AZEL -> RADEC (J2000)
+      os << "Creating converter from AZEL to RADEC (J2000)" << LogIO::POST ;
+      if ( mf_.null() ) {
+        mf_ = new MeasFrame() ;
+        vector<double> antpos = getAntennaPosition() ;
+        Vector<Quantity> qantpos( 3 ) ;
+        for ( int ip = 0 ; ip < 3 ; ip++ )
+          qantpos[ip] = Quantity( antpos[ip], "m" ) ;
+        MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
+        mf_->set( mp ) ;
+      }
+      converter_ = new MDirection::Convert( MDirection::AZEL,
+                                            MDirection::Ref(MDirection::J2000,
+                                                            *mf_ ) ) ;
+    }
+  }
+
+  if ( coord_ == 2 ) {
+    MEpoch me( Quantity( t, "d" ), MEpoch::UTC ) ;
+    mf_->set( me ) ;
+  }
 }
 
@@ -620 +610 @@
   spectra.resize( spec.size() ) ;
   int index = 0 ;
-  for ( Vector<Float>::iterator itr = spectra.begin() ; itr != spectra.end() ; itr++ ) {
-    *itr = spec[index++] ;
-  }
+  Bool b ;
+  Float *fp = spectra.getStorage( b ) ;
+  Float *wp = fp ;
+  for ( vector<double>::iterator i = spec.begin() ;
+        i != spec.end() ; i++ ) {
+    *wp = *i ;
+    wp++ ;
+  }
+  spectra.putStorage( fp, b ) ;
   //cout << "spec.size() = " << spec.size() << endl ;
   
   // flagtra
   bool setValue = !( flagtra.nelements() == spectra.nelements() ) ; 
-  flagtra.resize( spectra.nelements() ) ;
   if ( setValue ) {
     //cout << "flagtra resized. reset values..." << endl ;
+    flagtra.resize( spectra.nelements() ) ;
     flagtra.set( 0 ) ;
   }

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

@@ -48 +48 @@
 #include <measures/Measures/MeasFrame.h>
 #include <casa/Logging/LogIO.h>
+#include <casa/Utilities/CountedPtr.h>
 
 //#include <fitsio.h>
@@ -206 +207 @@
   // Get DIRECTION in RADEC(J2000)
   virtual Vector<Double> getDirection( int i ) ;
+  virtual void initConvert( int icoord, double t, char *epoch ) ;
 
   // filename 
@@ -212 +214 @@
   // dataset
   NRODataset *dataset_ ;
+
+  // source direction
+  Vector<Double> srcdir_ ;
+  Vector<Double> msrcdir_ ;
+
+  // for direction conversion
+  CountedPtr<MDirection::Convert> converter_ ;
+  CountedPtr<MeasFrame> mf_ ;
+  int coord_ ;
 
   // Logger