Changeset 1868 for trunk/external-alma/atnf/PKSIO/NRODataset.cc

Timestamp:

08/10/10 12:28:15 (14 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): atnf

Description: Describe your changes here...

Sync with code/atnf/implement/PKSIO

File:

• trunk/external-alma/atnf/PKSIO/NRODataset.cc (modified) (3 diffs)

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

@@ -34 +34 @@
 #include <casa/OS/Time.h>
 #include <scimath/Mathematics/InterpolateArray1D.h>
+
+#include <measures/Measures/MeasConvert.h>
+#include <measures/Measures/MCFrequency.h>
+#include <measures/Measures/MFrequency.h>
+#include <measures/Measures/MPosition.h>
+#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MDirection.h>
 
 #include <math.h>
@@ -796 +803 @@
   }
 
+  // conversion from TOPO to LSRK
+  v[1] = toLSR( v[1], getStartIntTime( i ), record_->SCX, record_->SCY ) ;
+
+  if ( refFreq_[ib] != 0.0 ) {
+    v[1] = refFreq_[ib] ;
+  }
+  else {
+    refFreq_[ib] = v[1] ;
+  }
+
   return v ;
 }
@@ -820 +837 @@
 
 }
+
+double NRODataset::toLSR( double v, double t, double x, double y ) 
+{
+  double vlsr ;
+
+  // get epoch
+  double tcent = t + 0.5*getIPTIM()/86400.0 ;
+  MEpoch me( Quantity( tcent, "d" ), MEpoch::UTC ) ;
+
+  // get antenna position
+  MPosition mp ;
+  if ( SITE.find( "45" ) != string::npos ) {
+    // 45m telescope
+    Double posx = -3.8710235e6 ;
+    Double posy = 3.4281068e6 ;
+    Double posz = 3.7240395e6 ;
+    mp = MPosition( MVPosition( posx, posy, posz ),
+                    MPosition::ITRF ) ;
+  }
+  else {
+    // ASTE
+    Vector<Double> pos( 2 ) ;
+    pos[0] = -67.7031 ;
+    pos[1] = -22.9717 ;
+    Double sitealt = 4800.0 ;
+    mp = MPosition( MVPosition( Quantity( sitealt, "m" ),
+                                Quantum< Vector<Double> >( pos, "deg" ) ),
+                    MPosition::WGS84 ) ;
+  }
+
+  // get direction 
+  MDirection md ;
+  if ( SCNCD == 0 ) {
+    // RADEC
+    if ( EPOCH == "B1950" ) {
+      md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
+                       MDirection::B1950 ) ;
+    }
+    else {
+      md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
+                       MDirection::J2000 ) ;
+    }
+  }
+  else if ( SCNCD == 1 ) {
+    // LB
+    md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
+                     MDirection::GALACTIC ) ;
+  }
+  else {
+    // AZEL
+    md = MDirection( Quantity( Double(x), "rad" ), Quantity( Double(y), "rad" ),
+                     MDirection::AZEL ) ;
+  }
+    
+  // to LSR
+  MeasFrame mf( me, mp, md ) ;
+  MFrequency::Convert tolsr( MFrequency::TOPO, MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
+  vlsr = (double)(tolsr( Double(v) ).get( "Hz" ).getValue()) ;
+
+  return vlsr ;
+}