diff --git a/sample/testkd2.cpp b/sample/testkd2.cpp index 375f276..0226ba0 100644 --- a/sample/testkd2.cpp +++ b/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; diff --git a/src/fortran.cpp b/src/fortran.cpp index 0664c51..e6383c0 100644 --- a/src/fortran.cpp +++ b/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; diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index 763d2c6..64b4970 100644 --- a/src/loadGadget.cpp +++ b/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(); diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index 250bb3d..01cd4c3 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -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) diff --git a/src/yorick.cpp b/src/yorick.cpp index 3453390..b6ebeed 100644 --- a/src/yorick.cpp +++ b/src/yorick.cpp @@ -212,6 +212,7 @@ namespace CosmoTool { template void loadArray(const char *fname, T*&array, uint32_t *&dimList, uint32_t& rank) + throw (NoSuchFileException) { NcFile f(fname, NcFile::ReadOnly); diff --git a/src/yorick.hpp b/src/yorick.hpp index 6a2cf94..f124ae8 100644 --- a/src/yorick.hpp +++ b/src/yorick.hpp @@ -183,9 +183,11 @@ namespace CosmoTool template void saveArray(const char *fname, T *array, uint32_t *dimList, uint32_t rank); + template 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); };