Changeset 1525 for tags/asap2.3.0/src/Scantable.cpp

Timestamp:

03/13/09 11:55:57 (15 years ago)

Author:

Malte Marquarding

Message:

Tagged 2.3.0 release

Location:

tags/asap2.3.0

Files:

• . (copied) (copied from trunk)
• src/Scantable.cpp (modified) (2 diffs)

tags/asap2.3.0/src/Scantable.cpp

@@ -49 +49 @@
 #include <measures/Measures/MCDirection.h>
 //
+#include <measures/Measures/MCRadialVelocity.h>
+#include <measures/Measures/MCDoppler.h>
 #include <measures/Measures/MDirection.h>
 #include <measures/Measures/MFrequency.h>
@@ -912 +914 @@
   if ( Unit(us) == Unit("Hz") ) {
     for ( int i=0; i < nchan; ++i) {
-      Double world;
-      spc.toWorld(world, pixel[i]);
-      stlout.push_back(double(world));
+      Double worldi;
+      spc.toWorld(worldi, pixel[i]);
+      stlout.push_back(double(worldi));
     }
   } else if ( Unit(us) == Unit("km/s") ) {
-    Vector<Double> world;
-    spc.pixelToVelocity(world, pixel);
-    world.tovector(stlout);
+    
+    MeasFrame frame(mp,md,me);
+    MRadialVelocity::Types refin, refout;
+    MFrequency::Types refinfreq = freqTable_.getFrame(true);
+    MRadialVelocity::getType(refin, freqTable_.getFrameString(true));
+    MRadialVelocity::getType(refout, freqTable_.getFrameString(false));
+    SpectralCoordinate spc1 =
+      freqTable_.getSpectralCoordinate(mfreqidCol_(whichrow));
+    
+    MRadialVelocity::Ref rin(refin, frame);
+    MRadialVelocity::Ref rout(refout, frame);
+    MRadialVelocity::Convert cnv(rin, rout);
+    Vector<Double> worlds(2);  
+    Double deltav;
+    for ( int i=0; i < nchan; ++i) {      
+      
+      Double worldi;
+      if (i == 0) {
+        for (uInt k = 0; k<2;++k) {
+          spc1.toWorld(worlds[i+k], pixel[i+k]);
+          MDoppler mdop =                               \
+            MFrequency(Quantity(worlds[i+k], "Hz"), 
+                       refinfreq).toDoppler(Quantity(rf, "Hz"));
+          MRadialVelocity mrv = MRadialVelocity::fromDoppler(mdop, refin);
+          MRadialVelocity mrvo =  cnv(mrv);
+          //worlds[i+k] = cnv(mrv).get(us).getValue();
+          
+            MDoppler mdopo = MDoppler::Convert(mrvo.toDoppler(),
+            MDoppler::RADIO)();
+            MFrequency mf = MFrequency::fromDoppler(mdopo, 
+            MVFrequency(Quantity(rf, "Hz")));
+            worlds[i+k] = mf.get("GHz").getValue();
+        
+        }
+      }      
+      deltav = worlds[1]-worlds[0];
+      worldi = worlds[0] + Double(i)*deltav;
+      stlout.push_back(double(worldi));
+    }
+    cout << worlds[0] << " + " << deltav <<  "+" << spc1.increment() << endl;
   }
   return stlout;