Merge branch 'master' of file:///home/guilhem/Dropbox/gitRoot/CosmoToolbox

sample/testkd2.cpp

@@ -6,7 +6,7 @@
 #define __KD_TREE_NUMNODES
 #include "mykdtree.hpp"
 
-#define NTRY 3
+#define NTRY 100 
 #define ND 3
 
 using namespace std;
@@ -62,11 +62,12 @@ int main()
   // Check consistency
   cout << "Check consistency..." << endl;
   MyCell **ngb = new MyCell *[12];
+  double *distances = new double[12];
   
   ofstream fngb("nearest.txt");
   for (int k = 0; k < NTRY; k++) {
     cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl;
-    tree.getNearestNeighbours(xc[k], 12, ngb);
+    tree.getNearestNeighbours(xc[k], 12, ngb, distances);
 
     for (uint32_t i = 0; i < 12; i++)
       {
@@ -75,6 +76,34 @@ int main()
 	  d2 += ({double delta = xc[k][l] - ngb[i]->coord[l]; delta*delta;});
 	fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl;
       }
+	fngb << endl << endl;
+    double farther_dist = distances[11];
+    for (uint32_t i = 0; i < Ncells; i++)
+      {
+         bool found = false;
+         // If the points is not in the list, it means it is farther than the farther point
+        for (int j =0; j < 12; j++)
+          {
+             if (&cells[i] == ngb[j]) {
+               found = true;
+               break;
+             }
+          } 
+            double dist_to_seed = 0;
+            for (int l = 0; l < 3; l++)
+		{ double delta = xc[k][l]-cells[i].coord[l];
+              dist_to_seed += delta*delta; }
+        if (!found) 
+          { 
+	    if (dist_to_seed <= farther_dist)
+	      abort();
+          }
+	else
+	  {
+	    if (dist_to_seed > farther_dist)
+		abort();
+	  }
+      }
   }
 
   return 0;

src/fortran.cpp

@@ -247,7 +247,7 @@ void UnformattedWrite::endCheckpoint()
   if (checkPointAccum == 0)
     throw InvalidUnformattedAccess();
 
-  ostream::streampos curPos = f->tellp();
+  streampos curPos = f->tellp();
 
   int64_t deltaPos = curPos-checkPointRef;
 

src/loadGadget.cpp

@@ -28,7 +28,7 @@ PurePositionData *CosmoTool::loadGadgetPosition(const char *fname)
     h.npartTotal[i] = f.readInt32();
   h.flag_cooling = f.readInt32();
   h.num_files = f.readInt32();
-  h.BoxSize = f.readReal64();
+  data->BoxSize = h.BoxSize = f.readReal64();
   h.Omega0 = f.readReal64();
   h.OmegaLambda = f.readReal64();
   h.HubbleParam = f.readReal64();

src/loadRamses.cpp

@@ -1,3 +1,4 @@
+#include <sys/types.h>
 #include <regex.h>
 #include <cmath>
 #include <cstring>
@@ -436,12 +437,13 @@ CosmoTool::PhaseSpaceData *CosmoTool::loadRamsesPhase(const char *basename, int
   if (!readInfoFile(basename, outputId, info))
     return 0;
 
-  double unit_vel = info.unitLength/info.unit_t/1e5;
   double hubble = info.aexp*info.aexp/info.unit_t / (1e5/CM_IN_MPC);
-  double L0 = info.boxSize*info.unitLength/CM_IN_MPC/info.aexp;
+  double L0 = info.boxSize*info.unitLength*hubble/(100*CM_IN_MPC)/info.aexp;
+  double unit_vel = 100*L0/info.aexp;
   if (!quiet) {
     cout << "L0=" << L0 << " Mpc" << endl;
     cout << "H=" << hubble << " km/s/Mpc" << endl;
+    cout << "unit_vel=" << unit_vel << " km/s" << endl;
   }
 
   if (!quiet)

src/yorick.cpp

@@ -212,6 +212,7 @@ namespace CosmoTool {
   template<typename T>
   void loadArray(const char *fname,
 		 T*&array, uint32_t *&dimList, uint32_t& rank)
+	throw (NoSuchFileException)
   {
     NcFile f(fname, NcFile::ReadOnly);
     

src/yorick.hpp

@@ -183,9 +183,11 @@ namespace CosmoTool
   template<typename T>
   void saveArray(const char *fname,
 		 T *array, uint32_t *dimList, uint32_t rank);
+
   template<typename T>
   void loadArray(const char *fname,
-		 T*& array, uint32_t *& dimList, uint32_t& rank);
+		 T*& array, uint32_t *& dimList, uint32_t& rank)
+    throw (NoSuchFileException);
 
   ProgressiveDoubleOutput saveDoubleArrayProgressive(const char *fname, uint32_t *dimList, uint32_t rank);  
 };