Changeset 2761 for trunk/external-alma/atnf/PKSIO/NROReader.cc

Timestamp:

02/06/13 13:13:20 (11 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: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Added new option, freqref, to NRO filler. Posible values are:
1) 'rest' to import frequency in REST frame, which results in an exactly
same frequency label as NEWSTAR, and 2) 'vref' to import frequency
in the frame that source velocity refers, which results in the same
velocity label as NEWSTAR. The option must be given to scantable
constructor.

File:

• trunk/external-alma/atnf/PKSIO/NROReader.cc (modified) (10 diffs)

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

@@ -46 +46 @@
 #include <measures/Measures/MeasConvert.h>
 
+#include <casa/Logging/LogIO.h>
 #include <casa/IO/RegularFileIO.h>
 #include <casa/OS/File.h>
@@ -185 +186 @@
   : dataset_( NULL ),
     srcdir_( 0 ),
-    msrcdir_( 0 )
+    msrcdir_( 0 ),
+    freqRefFromVREF_( false ) // import frequency as REST
 {
   // initialization
@@ -199 +201 @@
     dataset_ = NULL ;
   }
+}
+
+// set frequency reference frame
+void NROReader::setFreqRefFromVREF( bool fromVREF )
+{
+  os_.origin( LogOrigin( "NROReader", "setFreqRefFromVREF", WHERE ) ) ;
+  os_ << ((fromVREF) ? "Take frequency reference frame from VREF" : 
+          "Use frequency reference frame REST") << LogIO::POST ;
+
+  freqRefFromVREF_ = fromVREF ; 
 }
 
@@ -387 +399 @@
         for ( int ip = 0 ; ip < 3 ; ip++ )
           qantpos[ip] = Quantity( antpos[ip], "m" ) ;
-        MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
-        mf_->set( mp ) ;
+        mp_ = MPosition( MVPosition( qantpos ), MPosition::ITRF ) ;
+        mf_->set( mp_ ) ;
       }
       converter_ = new MDirection::Convert( MDirection::AZEL,
@@ -397 +409 @@
 
   if ( coord_ == 2 ) {
-    MEpoch me( Quantity( t, "d" ), MEpoch::UTC ) ;
-    mf_->set( me ) ;
+    me_ = MEpoch( Quantity( t, "d" ), MEpoch::UTC ) ;
+    mf_->set( me_ ) ;
   }
 }
@@ -458 +470 @@
     }
   }
+  else if ( vref.compare( 0, 3, "GAL" ) == 0 ) {
+    vref = "GALACTO" ;
+  }
+  else if (vref.compare( 0, 3, "HEL" ) == 0 ) {
+    os_.origin( LogOrigin( "NROReader", "getHeaderInfo", WHERE ) ) ;
+    os_ << LogIO::WARN << "Heliocentric frame is not supported. Use Barycentric frame instead." << LogIO::POST ;
+    vref = "BARY" ;
+  }
   //freqref = vref ;
   //freqref = "LSRK" ;
-  freqref = "REST" ;
+  //freqref = "REST" ;
+  freqref = freqRefFromVREF_ ? vref : "REST" ; 
   //cout << "freqref = " << freqref << endl ;
   NRODataRecord *record = dataset_->getRecord( 0 ) ;
-  reffreq = record->FREQ0 ;
+  reffreq = record->FREQ0 ; // rest frequency
+
   //cout << "reffreq = " << reffreq << endl ;
   bw = dataset_->getBEBW()[0] ;
@@ -538 +560 @@
                             Float &winddir,      
                             Double &srcvel,
-                            Vector<Double> &propermotion,
+                            Vector<Double> &/*propermotion*/,
                             Vector<Double> &srcdir,
-                            Vector<Double> &scanrate )
+                            Vector<Double> &/*scanrate*/ )
 {
   static const IPosition oneByOne( 1, 1 );
@@ -580 +602 @@
   //cout << "freqs = [" << freqs[0] << ", " << freqs[1] << ", " << freqs[2] << "]" << endl ;
 
+  if ( freqRefFromVREF_ ) {
+    freqs = shiftFrequency( freqs, 
+                            dataset_->getURVEL(), 
+                            dataset_->getVDEF() ) ;
+  }
+
   // restfreq (for MOLECULE_ID)
   restfreq.resize( oneByOne ) ;
@@ -609 +637 @@
   vector<double> spec = dataset_->getSpectrum( irow ) ;
   spectra.resize( spec.size() ) ;
-  int index = 0 ;
+  //int index = 0 ;
   Bool b ;
   Float *fp = spectra.getStorage( b ) ;
@@ -708 +736 @@
   return dataset_->getRowNum() ;
 }
+
+vector<double> NROReader::shiftFrequency( const vector<double> &f, 
+                                          const double &v, 
+                                          const string &vdef )
+{
+  vector<double> r( f ) ;
+  double factor = v / 2.99792458e8 ;
+  if ( vdef.compare( 0, 3, "RAD" ) == 0 ) {
+    factor = 1.0 / ( 1.0 + factor ) ;
+    r[1] *= factor ;
+    r[2] *= factor ;
+  }
+  else if ( vdef.compare( 0, 3, "OPT" ) == 0 ) {
+    factor += 1.0 ;
+    r[1] *= factor ;
+    r[2] *= factor ;
+  }
+  else {
+    cout << "vdef=" << vdef << " is not supported." << endl;
+  }
+  return r ;
+}