Ignore:
Timestamp:
07/14/11 15:02:29 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1913

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

Frequency is stored as LSRK value although raw frequency value is
in arbitrary frequency frame.
Also, direction is stored as J2000 value independently of original
direction reference code.


Location:
trunk/external-alma/asdm2ASAP
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/external-alma/asdm2ASAP/ASDMFiller.cc

    r2220 r2225  
    77#include <measures/Measures/MDirection.h>
    88#include <measures/Measures/MCDirection.h>
     9#include <measures/Measures/MFrequency.h>
     10#include <measures/Measures/MCFrequency.h>
    911#include <measures/Measures/MeasFrame.h>
    1012#include <measures/Measures/MeasConvert.h>
    … …  
    5153  antennaName_ = reader_->getAntennaName() ;
    5254
    53   //logsink->postLocally( LogMessage("antennaId_ = "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
    54   //logsink->postLocally( LogMessage("antennaName_ = "+antennaName_,LogOrigin(className_,funcName,WHERE)) ) ;
     55  //logsink_->postLocally( LogMessage("antennaId_ = "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
     56  //logsink_->postLocally( LogMessage("antennaName_ = "+antennaName_,LogOrigin(className_,funcName,WHERE)) ) ;
    5557
    5658  return status ;
    … …  
    8284  uInt numConfigDescId = configDescIdList.size() ;
    8385
    84   //logsink->postLocally( LogMessage("configDescIdList = "+String::toString(configDescIdList),LogOrigin(className_,funcName,WHERE)) ) ;
     86  //logsink_->postLocally( LogMessage("configDescIdList = "+String::toString(configDescIdList),LogOrigin(className_,funcName,WHERE)) ) ;
    8587
    8688  // get field list
    … …  
    8890  uInt numFieldId = fieldIdList.size() ;
    8991
    90   //logsink->postLocally( LogMessage("fieldIdList = "+String::toString(fieldIdList),LogOrigin(className_,funcName,WHERE)) ) ;
     92  //logsink_->postLocally( LogMessage("fieldIdList = "+String::toString(fieldIdList),LogOrigin(className_,funcName,WHERE)) ) ;
    9193
    9294  // BEAMNO is always 0 since ALMA antenna is single beam
    … …  
    187189          string fieldname = reader_->getFieldName( idata ) ;
    188190          int srctype = reader_->getSrcType( scanno, subscanno ) ;
    189           vector<double> srcDirection = reader_->getSourceDirection( idata ) ;
     191          //vector<double> srcDirection = reader_->getSourceDirection( idata ) ;
     192          vector<double> srcDirection ;
     193          string srcDirRef ;
     194          reader_->getSourceDirection( idata, srcDirection, srcDirRef ) ;
    190195          vector<double> srcProperMotion = reader_->getSourceProperMotion( idata ) ;
    191196          double sysVel = reader_->getSysVel( idata ) ;
    … …  
    343348                       hdr.poltype ) ;
    344349
     350  if ( hdr.freqref != "LSRK" ) {
     351//     String freqref ;
     352//     if (hdr.freqref == "TOPOCENT") {
     353//       freqref = "TOPO";
     354//     } else if (hdr.freqref == "GEOCENTR") {
     355//       freqref = "GEO";
     356//     } else if (hdr.freqref == "BARYCENT") {
     357//       freqref = "BARY";
     358//     } else if (hdr.freqref == "GALACTOC") {
     359//       freqref = "GALACTO";
     360//     } else if (hdr.freqref == "LOCALGRP") {
     361//       freqref = "LGROUP";
     362//     } else if (hdr.freqref == "CMBDIPOL") {
     363//       freqref = "CMB";
     364//     } else if (hdr.freqref == "SOURCE") {
     365//       freqref = "REST";
     366//     }
     367    vector<double> sdir ;
     368    string ref ;
     369    reader_->getSourceDirection( sdir, ref ) ;
     370    Vector<casa::Double> sourceDir( sdir ) ;
     371    hdr.reffreq = toLSRK( hdr.reffreq, hdr.freqref, hdr.utc, hdr.antennaposition, sdir, String(ref) ) ;
     372    hdr.freqref = "LSRK" ;
     373  }
     374
    345375  setHeader( hdr ) ;
    346376}
    … …  
    498528  azel[0] = az ;
    499529  azel[1] = el ;
    500   MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
    501   Vector<Quantity> qantpos( 3 ) ;
    502   qantpos[0] = Quantity( antpos[0], "m" ) ;
    503   qantpos[1] = Quantity( antpos[1], "m" ) ;
    504   qantpos[2] = Quantity( antpos[2], "m" ) ;
    505   MPosition mp( MVPosition( qantpos ),
    506                 MPosition::ITRF ) ;
    507 //   mp.print( os_.output() ) ;
    508   MeasFrame mf( me, mp ) ;
    509   MDirection::Convert toj2000( MDirection::AZELGEO,
    510                                MDirection::Ref( MDirection::J2000, mf ) ) ;
    511   dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
    512   //logsink->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
     530//   MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
     531//   Vector<Quantity> qantpos( 3 ) ;
     532//   qantpos[0] = Quantity( antpos[0], "m" ) ;
     533//   qantpos[1] = Quantity( antpos[1], "m" ) ;
     534//   qantpos[2] = Quantity( antpos[2], "m" ) ;
     535//   MPosition mp( MVPosition( qantpos ),
     536//                 MPosition::ITRF ) ;
     537// //   mp.print( os_.output() ) ;
     538//   MeasFrame mf( me, mp ) ;
     539//   MDirection::Convert toj2000( MDirection::AZELGEO,
     540//                                MDirection::Ref( MDirection::J2000, mf ) ) ;
     541//   dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
     542  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
     543  //logsink_->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
     544}
     545
     546Vector<casa::Double> ASDMFiller::toJ2000( Vector<casa::Double> dir,
     547                                          String dirref,
     548                                          casa::Double mjd,
     549                                          Vector<casa::Double> antpos )
     550{
     551  Vector<casa::Double> newd( dir ) ;
     552  if ( dirref != "J2000" ) {
     553    MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
     554    Vector<Quantity> qantpos( 3 ) ;
     555    qantpos[0] = Quantity( antpos[0], "m" ) ;
     556    qantpos[1] = Quantity( antpos[1], "m" ) ;
     557    qantpos[2] = Quantity( antpos[2], "m" ) ;
     558    MPosition mp( MVPosition( qantpos ),
     559                  MPosition::ITRF ) ;
     560    //   mp.print( os_.output() ) ;
     561    MeasFrame mf( me, mp ) ;
     562    MDirection::Types dirtype ;
     563    Bool b = MDirection::getType( dirtype, dirref ) ;
     564    if ( b ) {
     565      MDirection::Convert toj2000( dirtype,
     566                                   MDirection::Ref( MDirection::J2000, mf ) ) ;
     567      newd = toj2000( dir ).getAngle( "rad" ).getValue() ;
     568    }
     569  }
     570  return newd ;
    513571}
    514572
    … …  
    525583  return ftype ;
    526584}
     585
     586casa::Double ASDMFiller::toLSRK( casa::Double freq,
     587                                 String freqref,
     588                                 casa::Double utc,
     589                                 Vector<casa::Double> antpos,
     590                                 Vector<casa::Double> dir,
     591                                 String dirref )
     592{
     593  String funcName = "toLSRK" ;
     594
     595  //logsink_->postLocally( LogMessage("freqref = "+freqref,LogOrigin(className_,funcName,WHERE)) ) ;
     596  casa::Double newf = freq ;
     597  if ( freqref != "LSRK" ) {
     598    MEpoch me( Quantum<casa::Double>( utc, Unit("d") ), MEpoch::UTC ) ;
     599    Vector< Quantum<casa::Double> > antposQ( 3 ) ;
     600    for ( int i = 0 ; i < 3 ; i++ )
     601      antposQ[i] = Quantum<casa::Double>( antpos[i], Unit("m") ) ;
     602    MPosition mp( antposQ, MPosition::ITRF ) ;
     603    MDirection::Types dirtype ;
     604    Bool b = MDirection::getType( dirtype, dirref ) ;
     605    if ( !b )
     606      dirtype = MDirection::J2000 ;
     607    MDirection md( Quantum<casa::Double>( dir[0], Unit("rad") ),
     608                   Quantum<casa::Double>( dir[1], Unit("rad") ),
     609                   dirtype ) ;
     610    MeasFrame mf( me, mp, md ) ;
     611    MFrequency::Types freqtype ;
     612    b = MFrequency::getType( freqtype, freqref ) ;
     613    if ( !b )
     614      freqtype = MFrequency::TOPO ;
     615    MFrequency::Convert tolsr( freqtype,
     616                               MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
     617    newf = tolsr( Quantum<casa::Double>( freq, Unit("Hz") ) ).get( "Hz" ).getValue() ;
     618    //logsink_->postLocally( LogMessage("freq = "+String::toString(freq)+", newf = "+String::toString(newf),LogOrigin(className_,funcName,WHERE)) ) ;
     619  }
     620  return newf ;
     621}

  • trunk/external-alma/asdm2ASAP/ASDMFiller.h

    r2218 r2225  
    7474                casa::Vector<casa::Double> antpos ) ;
    7575
     76  // to J2000
     77  casa::Vector<casa::Double> toJ2000( casa::Vector<casa::Double> dir,
     78                                      casa::String dirref,
     79                                      casa::Double mjd,
     80                                      casa::Vector<casa::Double> antpos ) ;
     81
    7682  // get frequency frame enum value from string
    7783  casa::MFrequency::Types toFrameType( std::string &s ) ;
     84
     85  // to LSRK
     86  // utc must be UTC time in "d" (day)
     87  // antpos must be ITRF value in "m"
     88  casa::Double toLSRK( casa::Double freq,
     89                       casa::String freqref,
     90                       casa::Double utc,
     91                       casa::Vector<casa::Double> antpos,
     92                       casa::Vector<casa::Double> dir,
     93                       casa::String dirref ) ;
    7894
    7995  casa::CountedPtr<ASDMReader> reader_ ;

  • trunk/external-alma/asdm2ASAP/ASDMReader.cc

    r2220 r2225  
    327327  if ( spwrows[refidx]->isMeasFreqRefExists() ) {
    328328    string mfr = CFrequencyReferenceCode::name( spwrows[refidx]->getMeasFreqRef() ) ;
    329     if (mfr == "TOPO") {
    330       freqref = "TOPOCENT";
    331     } else if (mfr == "GEO") {
    332       freqref = "GEOCENTR";
    333     } else if (mfr == "BARY") {
    334       freqref = "BARYCENT";
    335     } else if (mfr == "GALACTO") {
    336       freqref = "GALACTOC";
    337     } else if (mfr == "LGROUP") {
    338       freqref = "LOCALGRP";
    339     } else if (mfr == "CMB") {
    340       freqref = "CMBDIPOL";
    341     } else if (mfr == "REST") {
    342       freqref = "SOURCE";
    343     }
    344 
     329//     if (mfr == "TOPO") {
     330//       freqref = "TOPOCENT";
     331//     } else if (mfr == "GEO") {
     332//       freqref = "GEOCENTR";
     333//     } else if (mfr == "BARY") {
     334//       freqref = "BARYCENT";
     335//     } else if (mfr == "GALACTO") {
     336//       freqref = "GALACTOC";
     337//     } else if (mfr == "LGROUP") {
     338//       freqref = "LOCALGRP";
     339//     } else if (mfr == "CMB") {
     340//       freqref = "CMBDIPOL";
     341//     } else if (mfr == "REST") {
     342//       freqref = "SOURCE";
     343//     }
     344    freqref = String( mfr ) ;
    345345  }
    346346  else {
    347347    // frequency reference is TOPOCENT by default
    348     freqref = "TOPOCENT" ;
     348    //freqref = "TOPOCENT" ;
     349    freqref = "TOPO" ;
    349350  }
    350351
    … …  
    977978  }
    978979  return dir ;
     980}
     981
     982void ASDMReader::getSourceDirection( unsigned int idx,
     983                                     vector<double> &dir,
     984                                     string &ref )
     985{
     986  dir.resize( 2 ) ;
     987  dir[0] = 0.0 ;
     988  dir[1] = 0.0 ;
     989  ref = "J2000" ;
     990  unsigned int index = dataIdList_[idx] ;
     991  //ArrayTimeInterval tint( vmsData_->v_time[index]*s2d, vmsData_->v_interval[index]*s2d ) ;
     992  double startSec = vmsData_->v_time[index] - 0.5 * vmsData_->v_interval[index] ;
     993  ArrayTimeInterval tint( startSec*s2d, vmsData_->v_interval[index]*s2d ) ;
     994  Tag ddtag( vmsData_->v_dataDescId[index], TagType::DataDescription ) ;
     995  Tag spwtag = asdm_->getDataDescription().getRowByKey(ddtag)->getSpectralWindowId() ;
     996  Tag ftag( vmsData_->v_fieldId[index], TagType::Field ) ;
     997  FieldRow *frow = asdm_->getField().getRowByKey( ftag ) ;
     998  string srcname ;
     999  if ( frow->isSourceIdExists() ) {
     1000    //logsink_->postLocally( LogMessage("sourceId exists",LogOrigin(className_,funcName,WHERE)) ) ;
     1001    int sid = frow->getSourceId() ;
     1002    SourceRow *srow = asdm_->getSource().getRowByKey( sid, tint, spwtag ) ;
     1003    vector<Angle> srcdir = srow->getDirection() ;
     1004    if ( srow->isDirectionCodeExists() ) {
     1005      ref = CDirectionReferenceCode::toString( srow->getDirectionCode() ) ;
     1006    }
     1007    dir[0] = srcdir[0].get() ;
     1008    dir[1] = srcdir[1].get() ;
     1009  }
     1010}
     1011
     1012void ASDMReader::getSourceDirection( vector<double> &dir, string &ref )
     1013{
     1014  dir.resize( 2 ) ;
     1015  ref = "J2000" ; // default is J2000
     1016  SourceTable &tab = asdm_->getSource() ;
     1017  SourceRow *row = tab.get()[0] ;
     1018  vector<Angle> dirA = row->getDirection() ;
     1019  dir[0] = dirA[0].get() ;
     1020  dir[1] = dirA[1].get() ;
     1021  if ( row->isDirectionCodeExists() ) {
     1022    ref = CDirectionReferenceCode::toString( row->getDirectionCode() ) ;
     1023  }
    9791024}
    9801025
    … …  
    16181663  String funcName = "toJ2000" ;
    16191664
    1620   casa::Vector<casa::Double> azel( 2 ) ;
     1665//   casa::Vector<casa::Double> azel( 2 ) ;
     1666//   azel[0] = az ;
     1667//   azel[1] = el ;
     1668//   casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
     1669//   casa::Vector<casa::Quantity> qantpos( 3 ) ;
     1670//   qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
     1671//   qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
     1672//   qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
     1673//   casa::MPosition mp( casa::MVPosition( qantpos ),
     1674//                       casa::MPosition::ITRF ) ;
     1675//   //ostringstream oss ;
     1676//   //mp.print( oss ) ;
     1677//   //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
     1678
     1679//   casa::MeasFrame mf( me, mp ) ;
     1680//   casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
     1681//                                      casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
     1682//   casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
     1683//   //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
     1684//   dir.resize(2) ;
     1685//   dir[0] = (double)(cdir[0]) ;
     1686//   dir[1] = (double)(cdir[1]) ;
     1687  vector<double> azel( 2 ) ;
    16211688  azel[0] = az ;
    16221689  azel[1] = el ;
    1623   casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
    1624   casa::Vector<casa::Quantity> qantpos( 3 ) ;
    1625   qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
    1626   qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
    1627   qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
    1628   casa::MPosition mp( casa::MVPosition( qantpos ),
    1629                       casa::MPosition::ITRF ) ;
    1630   //ostringstream oss ;
    1631   //mp.print( oss ) ;
    1632   //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
    1633 
    1634   casa::MeasFrame mf( me, mp ) ;
    1635   casa::MDirection::Convert toj2000( casa::MDirection::AZELGEO,
    1636   //MDirection::Convert toj2000( MDirection::AZEL,
    1637   //MDirection::Convert toj2000( MDirection::AZELSW,
    1638   //MDirection::Convert toj2000( MDirection::AZELSWGEO,
    1639   //MDirection::Convert toj2000( MDirection::AZELNE,
    1640   //MDirection::Convert toj2000( MDirection::AZELNEGEO,
    1641                                      casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
    1642   casa::Vector<casa::Double> cdir = toj2000( azel ).getAngle( "rad" ).getValue() ;
    1643   //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
    1644   dir.resize(2) ;
    1645   dir[0] = (double)(cdir[0]) ;
    1646   dir[1] = (double)(cdir[1]) ;
     1690  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
     1691}
     1692
     1693vector<double> ASDMReader::toJ2000( vector<double> dir,
     1694                                    casa::String dirref,
     1695                                    double mjd,
     1696                                    casa::Vector<casa::Double> antpos )
     1697{
     1698  casa::String funcname = "toJ2000" ;
     1699
     1700  vector<double> newd( dir ) ;
     1701  if ( dirref != "J2000" ) {
     1702    casa::MEpoch me( casa::Quantity( (casa::Double)mjd, "d" ), casa::MEpoch::UTC ) ;
     1703    casa::Vector<casa::Quantity> qantpos( 3 ) ;
     1704    qantpos[0] = casa::Quantity( antpos[0], "m" ) ;
     1705    qantpos[1] = casa::Quantity( antpos[1], "m" ) ;
     1706    qantpos[2] = casa::Quantity( antpos[2], "m" ) ;
     1707    casa::MPosition mp( casa::MVPosition( qantpos ),
     1708                        casa::MPosition::ITRF ) ;
     1709    //ostringstream oss ;
     1710    //mp.print( oss ) ;
     1711    //logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
     1712   
     1713    casa::MeasFrame mf( me, mp ) ;
     1714    casa::MDirection::Types dirtype ;
     1715    casa::Bool b = casa::MDirection::getType( dirtype, dirref ) ;
     1716    if ( b ) {
     1717      casa::MDirection::Convert toj2000( dirtype,
     1718                                         casa::MDirection::Ref( casa::MDirection::J2000, mf ) ) ;
     1719      casa::Vector<casa::Double> cdir = toj2000( dir ).getAngle( "rad" ).getValue() ;
     1720      //logsink_->postLocally( LogMessage("cdir = "+String::toString(cdir),LogOrigin(className_,funcName,WHERE)) ) ;
     1721      newd[0] = (double)(cdir[0]) ;
     1722      newd[1] = (double)(cdir[1]) ;
     1723    }
     1724  }
     1725  return newd ;
    16471726}
    16481727

  • trunk/external-alma/asdm2ASAP/ASDMReader.h

    r2220 r2225  
    237237   **/
    238238  std::vector<double> getSourceDirection( unsigned int idx ) ;
     239  void getSourceDirection( unsigned int idx,
     240                           std::vector<double> &dir,
     241                           std::string &ref  ) ;
     242
     243  /**
     244   * get source direction with reference
     245   *
     246   * @param dir source direction
     247   * @param reference frame
     248   * @param idx for Source table
     249   **/
     250  void getSourceDirection( std::vector<double> &dir, std::string &ref ) ;
    239251 
    240252  /**
    … …  
    496508                casa::Vector<casa::Double> antpos ) ;
    497509
     510  /**
     511  * to J2000
     512  *
     513  * @param dir pointing direction
     514  * @param dirref direction reference
     515  * @param mjd reference time
     516  * @param antpos antenna position vector
     517  * @return new direction
     518  **/
     519  std::vector<double> toJ2000( std::vector<double> dir,
     520                               casa::String dirref,
     521                               double mjd,
     522                               casa::Vector<casa::Double> antpos ) ;
    498523  /**
    499524   * get nIF
Note: See TracChangeset for help on using the changeset viewer.