Changeset 2761


Ignore:
Timestamp:
02/06/13 13:13:20 (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: 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.


Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r2643 r2761  
    4646#include <measures/Measures/MeasConvert.h>
    4747
     48#include <casa/Logging/LogIO.h>
    4849#include <casa/IO/RegularFileIO.h>
    4950#include <casa/OS/File.h>
    … …  
    185186  : dataset_( NULL ),
    186187    srcdir_( 0 ),
    187     msrcdir_( 0 )
     188    msrcdir_( 0 ),
     189    freqRefFromVREF_( false ) // import frequency as REST
    188190{
    189191  // initialization
    … …  
    199201    dataset_ = NULL ;
    200202  }
     203}
     204
     205// set frequency reference frame
     206void NROReader::setFreqRefFromVREF( bool fromVREF )
     207{
     208  os_.origin( LogOrigin( "NROReader", "setFreqRefFromVREF", WHERE ) ) ;
     209  os_ << ((fromVREF) ? "Take frequency reference frame from VREF" :
     210          "Use frequency reference frame REST") << LogIO::POST ;
     211
     212  freqRefFromVREF_ = fromVREF ;
    201213}
    202214
    … …  
    387399        for ( int ip = 0 ; ip < 3 ; ip++ )
    388400          qantpos[ip] = Quantity( antpos[ip], "m" ) ;
    389         MPosition mp( MVPosition( qantpos ), MPosition::ITRF ) ;
    390         mf_->set( mp ) ;
     401        mp_ = MPosition( MVPosition( qantpos ), MPosition::ITRF ) ;
     402        mf_->set( mp_ ) ;
    391403      }
    392404      converter_ = new MDirection::Convert( MDirection::AZEL,
    … …  
    397409
    398410  if ( coord_ == 2 ) {
    399     MEpoch me( Quantity( t, "d" ), MEpoch::UTC ) ;
    400     mf_->set( me ) ;
     411    me_ = MEpoch( Quantity( t, "d" ), MEpoch::UTC ) ;
     412    mf_->set( me_ ) ;
    401413  }
    402414}
    … …  
    458470    }
    459471  }
     472  else if ( vref.compare( 0, 3, "GAL" ) == 0 ) {
     473    vref = "GALACTO" ;
     474  }
     475  else if (vref.compare( 0, 3, "HEL" ) == 0 ) {
     476    os_.origin( LogOrigin( "NROReader", "getHeaderInfo", WHERE ) ) ;
     477    os_ << LogIO::WARN << "Heliocentric frame is not supported. Use Barycentric frame instead." << LogIO::POST ;
     478    vref = "BARY" ;
     479  }
    460480  //freqref = vref ;
    461481  //freqref = "LSRK" ;
    462   freqref = "REST" ;
     482  //freqref = "REST" ;
     483  freqref = freqRefFromVREF_ ? vref : "REST" ;
    463484  //cout << "freqref = " << freqref << endl ;
    464485  NRODataRecord *record = dataset_->getRecord( 0 ) ;
    465   reffreq = record->FREQ0 ;
     486  reffreq = record->FREQ0 ; // rest frequency
     487
    466488  //cout << "reffreq = " << reffreq << endl ;
    467489  bw = dataset_->getBEBW()[0] ;
    … …  
    538560                            Float &winddir,     
    539561                            Double &srcvel,
    540                             Vector<Double> &propermotion,
     562                            Vector<Double> &/*propermotion*/,
    541563                            Vector<Double> &srcdir,
    542                             Vector<Double> &scanrate )
     564                            Vector<Double> &/*scanrate*/ )
    543565{
    544566  static const IPosition oneByOne( 1, 1 );
    … …  
    580602  //cout << "freqs = [" << freqs[0] << ", " << freqs[1] << ", " << freqs[2] << "]" << endl ;
    581603
     604  if ( freqRefFromVREF_ ) {
     605    freqs = shiftFrequency( freqs,
     606                            dataset_->getURVEL(),
     607                            dataset_->getVDEF() ) ;
     608  }
     609
    582610  // restfreq (for MOLECULE_ID)
    583611  restfreq.resize( oneByOne ) ;
    … …  
    609637  vector<double> spec = dataset_->getSpectrum( irow ) ;
    610638  spectra.resize( spec.size() ) ;
    611   int index = 0 ;
     639  //int index = 0 ;
    612640  Bool b ;
    613641  Float *fp = spectra.getStorage( b ) ;
    … …  
    708736  return dataset_->getRowNum() ;
    709737}
     738
     739vector<double> NROReader::shiftFrequency( const vector<double> &f,
     740                                          const double &v,
     741                                          const string &vdef )
     742{
     743  vector<double> r( f ) ;
     744  double factor = v / 2.99792458e8 ;
     745  if ( vdef.compare( 0, 3, "RAD" ) == 0 ) {
     746    factor = 1.0 / ( 1.0 + factor ) ;
     747    r[1] *= factor ;
     748    r[2] *= factor ;
     749  }
     750  else if ( vdef.compare( 0, 3, "OPT" ) == 0 ) {
     751    factor += 1.0 ;
     752    r[1] *= factor ;
     753    r[2] *= factor ;
     754  }
     755  else {
     756    cout << "vdef=" << vdef << " is not supported." << endl;
     757  }
     758  return r ;
     759}

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

    r2643 r2761  
    4242#include <casa/BasicSL/String.h>
    4343#include <measures/Measures/MPosition.h>
     44#include <measures/Measures/MEpoch.h>
    4445#include <measures/Measures/MCPosition.h>
    4546#include <measures/Measures/MDirection.h>
    … …  
    5657#include <atnf/PKSIO/NRODataset.h>
    5758#include <atnf/PKSIO/NRODataRecord.h>
    58 
    59 using namespace std ;
    6059
    6160// <summary>
    … …  
    103102  // Destructor.
    104103  virtual ~NROReader() ;
     104
     105  // determine whether to import frequency as REST (frequency is same as
     106  // NEWSTAR/NOSTAR) or as is (velocity is same as NEWSTAR/NOSTAR)
     107  void setFreqRefFromVREF( bool fromVREF ) ;
    105108
    106109  // Read data header
    … …  
    133136                           uInt &beamno,
    134137                           uInt &polno,
    135                            vector<double> &freqs,   
     138                           std::vector<double> &freqs,   
    136139                           Vector<Double> &restfreq, 
    137140                           uInt &refbeamno,
    … …  
    171174
    172175  // Get IF settings
    173   virtual vector<Bool> getIFs() ;
     176  virtual std::vector<Bool> getIFs() ;
    174177
    175178  // Get Number of IFs
    … …  
    177180
    178181  // Get Beam settings
    179   virtual vector<Bool> getBeams() ;
     182  virtual std::vector<Bool> getBeams() ;
    180183
    181184  // Get Number of Beams
    … …  
    188191
    189192  // Get spectrum
    190   virtual vector< vector<double> > getSpectrum() ;
     193  virtual std::vector< std::vector<double> > getSpectrum() ;
    191194
    192195  // Get number of polarization
    … …  
    196199  virtual double getStartTime() ;
    197200  virtual double getEndTime() ;
    198   virtual vector<double> getStartIntTime() ;
     201  virtual std::vector<double> getStartIntTime() ;
    199202  //virtual double getStartIntTime( int i ) ;
    200203
    201204  // Get Antenna Position in ITRF coordinate
    202   virtual vector<double> getAntennaPosition() = 0 ;
     205  virtual std::vector<double> getAntennaPosition() = 0 ;
    203206
    204207  // Get SRCDIRECTION in RADEC(J2000)
    … …  
    209212  virtual void initConvert( int icoord, double t, char *epoch ) ;
    210213
     214  // Shift frequency by given velocity with respect to specified
     215  // velocity reference
     216  std::vector<double> shiftFrequency( const std::vector<double> &f, const double &v, const string &vref ) ;
     217
    211218  // filename
    212219  string filename_ ;
    … …  
    222229  CountedPtr<MDirection::Convert> converter_ ;
    223230  CountedPtr<MeasFrame> mf_ ;
     231  MEpoch me_ ;
     232  MPosition mp_ ;
    224233  int coord_ ;
    225234
     235  bool freqRefFromVREF_ ;
     236
    226237  // Logger
    227   //LogIO os ;
     238  LogIO os_ ;
    228239};
    229240

  • trunk/python/scantable.py

    r2754 r2761  
    255255                    self._fill([filename], unit, average, opts)
    256256                elif os.path.isfile(filename):
    257                     self._fill([filename], unit, average)
     257                    opts={'nro': {}}
     258                    nrokeys=['freqref']
     259                    for key in nrokeys:
     260                        if key in args.keys():
     261                            opts['nro'][key] = args[key]
     262                    self._fill([filename], unit, average, opts)
    258263                    # only apply to new data not "copy constructor"
    259264                    self.parallactify(parallactify)

  • trunk/src/NROFiller.cpp

    r2289 r2761  
    2323#include <casa/Quanta/MVTime.h>
    2424#include <atnf/PKSIO/SrcType.h>
     25#include <casa/Logging/LogIO.h>
    2526
    2627using namespace casa;
    … …  
    3839}
    3940
    40   bool NROFiller::open(const std::string& filename, const Record& rec)
     41bool NROFiller::open(const std::string& filename, const Record& rec)
    4142{
    4243  bool status = true ;
     44
     45  // Parse options
     46  String freqref = "DEFAULT (REST)" ;
     47  if ( rec.isDefined( "nro" ) ) {
     48    Record nrorec = rec.asRecord( "nro" ) ;
     49    if ( nrorec.isDefined( "freqref" ) ) {
     50      freqref = nrorec.asString( "freqref" ) ;
     51      freqref.upcase() ;
     52    }
     53    LogIO os( LogOrigin( "NROFiller", "open", WHERE ) ) ;
     54    os << "Parsing NRO options" << endl ;
     55    os << "   freqref = " << freqref << LogIO::POST ;
     56  }
    4357
    4458  // get appropriate reader object
    … …  
    4963    return status ;
    5064  } 
     65
     66  // Apply options
     67  if ( freqref == "REST" ) {
     68    reader_->setFreqRefFromVREF( false ) ;
     69  }
     70  else if ( freqref == "VREF" ) {
     71    reader_->setFreqRefFromVREF( true ) ;
     72  }
    5173
    5274  // get header information
Note: See TracChangeset for help on using the changeset viewer.