From 9c52937c74330dd39079f56e7afcce4225f96cc1 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem@mithril.(none)>
Date: Tue, 29 Sep 2009 17:04:10 -0500
Subject: [PATCH] Fixed velocity scaling

---
 src/loadRamses.cpp | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp
index 0eed042..d15566a 100644
--- a/src/loadRamses.cpp
+++ b/src/loadRamses.cpp
@@ -194,7 +194,7 @@ typedef struct
   double unitLength;
   double aexp;
   double boxSize;
-  double hubble;
+  double unit_t;
 } InfoData;
 
 int readInfoFile(const char *basename, int outputId, InfoData& info)
@@ -249,10 +249,7 @@ int readInfoFile(const char *basename, int outputId, InfoData& info)
   info.unitLength = infoMap["unit_l"];
   info.aexp = infoMap["aexp"];
   info.boxSize = infoMap["boxlen"];
-  
-  static const double CM_IN_MPC = 3.08e24;
-  double unit_t = infoMap["unit_t"];
-  info.hubble = info.aexp*info.aexp/unit_t / (1e5/CM_IN_MPC);
+  info.unit_t = infoMap["unit_t"];
 
   return 1;
 }
@@ -270,11 +267,12 @@ CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename,
   if (!readInfoFile(basename, outputId, info))
     return 0;
 
+  double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
   double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp;
   if (!quiet)
     {
       cout << "L0=" << L0 << " Mpc" << endl;
-      cout << "H=" << info.hubble << " km/s/Mpc" << endl;
+      cout << "H=" << hubble << " km/s/Mpc" << endl;
     }
 
   if (!quiet)
@@ -328,7 +326,7 @@ CosmoTool::PurePositionData *CosmoTool::loadRamsesPosition(const char *basename,
   assert(gd->pos != 0);
   gd->NumPart = totPart;
   gd->BoxSize = L0*1000; 
-  gd->hubble = info.hubble;
+  gd->hubble = hubble;
 
   if (!quiet)
     cout << " Total number part=" << totPart << endl
@@ -438,9 +436,13 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int
   if (!readInfoFile(basename, outputId, info))
     return 0;
 
+  double unit_vel = info.unitLength/info.unit_t;
+  double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
   double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp;
-  if (!quiet)
+  if (!quiet) {
     cout << "L0=" << L0 << " Mpc" << endl;
+    cout << "H=" << hubble << " km/s/Mpc" << endl;
+  }
 
   if (!quiet)
     cout << "Detecting number of files and particles..." << endl;
@@ -495,7 +497,8 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int
   assert(gd->vel != 0);
   gd->NumPart = totPart;
   gd->BoxSize = L0*1000; 
-
+  gd->hubble = hubble;
+  
   if (!quiet)
     cout << " Total number part=" << totPart << endl
 	 << "Loading particles ..." << endl;
@@ -551,7 +554,7 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int
 	    infile.beginCheckpoint();
 	    for (uint32_t i = 0; i < nPar; i++)
 	      {
-		vel[i][k] = infile.readReal32()*gd->BoxSize/1000*100;
+		vel[i][k] = infile.readReal32()*unit_vel;
 	      }
 	    infile.endCheckpoint();
 	  }