renamed src to source

c_source/CMakeLists.txt

@@ -0,0 +1,47 @@
+
+include(FindPkgConfig)
+
+pkg_check_modules(CAIROMM cairomm-1.0)
+pkg_check_modules(SIGC sigc++-2.0)
+pkg_check_modules(CAIRO cairo)
+pkg_check_modules(FREETYPE freetype2)
+
+find_path(CAIROMMCONFIG_INCLUDE_PATH NAMES cairommconfig.h HINTS ${CAIROMM_INCLUDE_DIRS})
+find_path(CAIROMM_INCLUDE_PATH NAMES cairomm/cairomm.h HINTS ${CAIROMM_INCLUDE_DIRS})
+find_path(SIGC_INCLUDE_PATH NAMES sigc++/slot.h HINTS ${SIGC_INCLUDE_DIRS})
+find_path(SIGCCONFIG_INCLUDE_PATH NAMES sigc++config.h HINTS ${SIGC_INCLUDE_DIRS})
+find_path(CAIRO_INCLUDE_PATH NAMES cairo.h HINTS ${CAIRO_INCLUDE_DIRS})
+find_path(FREETYPE_INCLUDE_PATH NAMES freetype/config/ftheader.h HINTS ${FREETYPE_INCLUDE_DIRS})
+
+find_library(CAIROMM_LIB NAMES ${CAIROMM_LIBRARIES} HINTS CAIROMM_LIBRARY_DIRS)
+
+IF (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB)
+
+  SET(CAIRO_FOUND 1)
+  SET(ALL_CAIROMM_LIBS ${CAIROMM_LIB})
+  SET(CAIRO_HEADERS ${CAIROMM_INCLUDE_PATH} 
+    ${CAIROMMCONFIG_INCLUDE_PATH} ${SIGC_INCLUDE_PATH} 
+    ${SIGCCONFIG_INCLUDE_PATH} ${CAIRO_INCLUDE_PATH}
+    ${FREETYPE_INCLUDE_PATH})
+
+ELSE (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB)
+
+  SET(CAIRO_FOUND 0)
+
+ENDIF (CAIROMMCONFIG_INCLUDE_PATH AND CAIROMM_INCLUDE_PATH AND SIGC_INCLUDE_PATH AND SIGCCONFIG_INCLUDE_PATH AND CAIRO_INCLUDE_PATH AND FREETYPE_INCLUDE_PATH AND CAIROMM_LIB)
+
+
+include_directories(
+  ${QHULL_INCLUDES}
+  ${HEALPIX_INCLUDE_PATH}
+  ${CMAKE_CURRENT_BINARY_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/util)
+
+SET(ZOB_LIBS zobovTool)
+
+set(computeAverageDistortion_SRCS computeAverageDistortion.cpp)
+add_genopt(computeAverageDistortion_SRCS computeAverageDistortion.ggo computeAverageDistortion_conf STRUCTNAME Params)
+#add_executable(computeAverageDistortion ${computeAverageDistortion_SRCS})
+#target_link_libraries(computeAverageDistortion ${ZOB_LIBS})
+
+subdirs(util prep pruning jozov2 zobov)

c_source/jozov2/CMakeLists.txt

@@ -0,0 +1,26 @@
+function(omp_add_flags TARGET LANG)
+   if(OPENMP_FOUND)
+     if(NOT LANG)
+       get_target_property(LANG ${TARGET} LINKER_LANGUAGE)
+     endif()
+     if(LANG MATCHES "C(XX)?")
+       set_property(TARGET ${TARGET} APPEND
+                    PROPERTY COMPILE_FLAGS ${OpenMP_${LANG}_FLAGS})
+       set_property(TARGET ${TARGET} APPEND
+                    PROPERTY LINK_FLAGS ${OpenMP_${LANG}_FLAGS})
+     else()
+       message(WARNING "omp_add_flags: target '${TARGET}'"
+                      " link language '${LANG}' is not supported")
+    endif()
+  endif()
+endfunction()
+
+add_executable(jozov2 jozov2.cpp jozov2_io.cpp jozov2_zones.cpp jozov2_watershed.cpp findrtop.c)
+
+if (ENABLE_OPENMP)
+  Find_Package(OpenMP)
+  omp_add_flags(jozov2 CXX)
+  add_definitions(-DOPENMP)
+ENDIF(ENABLE_OPENMP)
+
+install(TARGETS jozov2 DESTINATION ${VIDE_BIN})

c_source/jozov2/findrtop.c

@@ -0,0 +1,81 @@
+#include <math.h>
+#include <stdlib.h>
+
+/*------------------------------------------------------------------------------
+  Find nb elements of Real array a having the largest value.
+  Returns index iord of these elements, ordered so iord[0] corresponds
+  to element a[iord[0]] having the largest value.
+  If nb > na, then the last nb-na elements of iord are undefined.
+  Elements of a that are equal are left in their original order.
+
+  Courtesy Andrew Hamilton.
+*/
+void findrtop(double *a, int na, int *iord, int nb)
+{
+#undef	ORDER
+#define	ORDER(ia, ja)	a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja)
+
+  int i, ia, ib, it, n;
+  
+  n = (na <= nb)? na : nb;
+  if (n <= 0) return;
+  
+  /* heap first n elements, so smallest element is at top of heap */
+  for (ib = (n >> 1); ib < n; ib++) {
+    iord[ib] = ib;
+  }
+  for (ia = (n >> 1) - 1; ia >= 0; ia--) {
+    i = ia;
+    for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+      if (ib+1 < n) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(ia, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = ia;
+  }
+  
+  /* now compare rest of elements of array to heap */
+  for (ia = n; ia < na; ia++) {
+    /* if new element is greater than smallest, sift it into heap */
+    i = 0;
+    if (ORDER(ia, iord[i])) {
+      for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+	if (ib+1 < n) {
+	  if (ORDER(iord[ib], iord[ib+1])) ib++;
+	}
+	if (ORDER(ia, iord[ib])) {
+	  iord[i] = iord[ib];
+	  i = ib;
+	} else {
+	  break;
+	}
+      }
+      iord[i] = ia;
+    }
+  }
+  
+  /* unheap iord so largest element is at top */
+  for (ia = n - 1; ia > 0; ia--) {
+    it = iord[ia];
+    i = 0;
+    iord[ia] = iord[i];
+    for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) {
+      if (ib+1 < ia) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(it, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = it;
+  }
+}

c_source/jozov2/jozov2.cpp

@@ -0,0 +1,174 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#include <boost/format.hpp>
+#include <exception>
+#include <iostream>
+#include <fstream>
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <assert.h>
+/* jovoz.c by Mark Neyrinck */
+#include "jozov2.hpp"
+#include "zobov.hpp"
+
+using namespace std;
+using boost::format;
+
+
+
+int main(int argc,char **argv)
+{
+  ofstream txt;
+  PARTICLE *p;
+  pid_t np;
+
+  ZONE *z;
+  int nzones;
+
+  string adjfile, volfile, zonfile, zonfile2, txtfile;
+  int *jumped;
+  int *zonenum;
+  float voltol, prob;
+
+  float maxvol, minvol;
+  double *sorter;
+  int *iord;
+
+  int mockIndex;
+
+  if (argc != 8) {
+    printf("Wrong number of arguments.\n");
+    printf("arg1: adjacency file\n");
+    printf("arg2: volume file\n");
+    printf("arg3: output file containing particles in each zone\n");
+    printf("arg4: output file containing zones in each void\n");
+    printf("arg5: output text file\n");
+    printf("arg6: Density threshold (0 for no threshold)\n");
+    printf("arg7: Beginning index of mock galaxies\n\n");
+    exit(0);
+  }
+  adjfile = argv[1];
+  volfile = argv[2];
+  zonfile = argv[3];
+  zonfile2 = argv[4];
+  txtfile = argv[5];
+  if (sscanf(argv[6],"%f",&voltol) == 0) {
+    printf("Bad density threshold.\n");
+    exit(0);
+  }
+  if (sscanf(argv[7],"%d",&mockIndex) == 0) {
+    printf("Bad mock galaxy index.\n");
+    exit(0);
+  }
+  printf("TOLERANCE: %f\n", voltol);
+  if (voltol <= 0.) {
+    printf("Proceeding without a density threshold.\n");
+    voltol = 1e30;
+  }
+
+  try
+    {
+      readAdjacencyFile(adjfile, p, np);
+    }
+  catch(const FileError& e)
+    {
+      return 1;
+    }
+  if (mockIndex < 0)
+    mockIndex = np;
+
+  /* Check that we got all the pairs */
+  for (int i = 0; i < np; i++)
+    {
+      if (p[i].ncnt != p[i].nadj && i < mockIndex) {
+	cout
+	  << format("We didn't get all of %d's adj's; %d != %d.") 
+	  % i % p[i].ncnt % p[i].nadj
+	  << endl;
+	p[i].nadj = p[i].ncnt;
+      }
+    }
+
+  /* Volumes */
+  try
+    {
+      readVolumeFile(volfile, p, np, mockIndex);
+    }
+  catch (const FileError& e)
+    {
+      return 2;
+    }
+
+  buildZones(p, np, jumped, z, nzones, zonenum);
+  writeZoneFile(zonfile, p, np, z, nzones, zonenum, jumped);
+
+  maxvol = 0.;
+  minvol = BIGFLT;
+  for(pid_t i = 0; i < np; i++){
+    if (p[i].dens > maxvol)
+      maxvol = p[i].dens;
+    if (p[i].dens < minvol) 
+      minvol = p[i].dens;
+  }
+  cout << format("Densities range from %e to %e.") % minvol % maxvol << endl;
+  FF;
+
+  doWatershed(p, np, z, nzones, maxvol, voltol);
+  writeVoidFile(zonfile2, z, nzones);
+
+  sorter = new double[nzones+1];
+  /* Assign sorter by probability (could use volume instead) */
+  for (int h = 0; h < nzones; h++)
+    sorter[h] = (double)z[h].denscontrast;
+    
+  /* Text output file */
+
+  printf("about to sort ...\n");FF;
+
+  iord = new int[nzones];
+
+  findrtop(sorter, nzones, iord, nzones);
+  delete[] sorter;
+
+  txt.open(txtfile.c_str());
+  txt << format("%d particles, %d voids.") % np % nzones << endl;
+  txt << "Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb" << endl;
+  for (int h=0; h<nzones; h++)
+    {
+      int i = iord[h];
+      prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
+      if (z[i].np == 1) 
+	continue;
+      
+      txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e")
+	   % (h+1) % i % z[i].core % p[z[i].core].dens % z[i].vol
+           % z[i].np % z[i].nhl % z[i].voljoin % z[i].npjoin
+	   % z[i].denscontrast % prob << endl;
+
+    } /* h+1 to start from 1, not zero */
+  txt.close();
+
+  delete[] iord;
+  delete[] z;
+  delete[] p;
+
+  cout << "Done!" << endl;
+  return(0);
+}

c_source/jozov2/jozov2.hpp

@@ -0,0 +1,65 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#ifndef __JOZOV2_HPP
+#define __JOZOV2_HPP
+
+#include <string>
+#include <exception>
+#include "zobov.hpp"
+
+#define BIGFLT 1e30 /* Biggest possible floating-point number */
+#define NLINKS 1000 /* Number of possible links with the same rho_sl */
+#define FF cout.flush()
+
+class FileError: virtual std::exception
+{
+};
+
+void readAdjacencyFile(const std::string& adjfile, PARTICLE*& p, pid_t& np)
+  throw(FileError);
+
+void readVolumeFile(const std::string& volfile, PARTICLE *p, pid_t np, 
+                    pid_t mockIndex)
+  throw(FileError);
+
+void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped, 
+                       pid_t *numinh, pid_t& numZones);
+
+void buildZoneAdjacencies(PARTICLE *p, pid_t np,
+                          ZONE *z, ZONET *zt,
+                          int numZones,
+                          pid_t *jumped,
+                          int *zonenum,
+                          int *numinh);
+
+void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped,
+                ZONE*& z, int& nzones,
+                int*& zonenum);
+
+void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol);
+
+void writeZoneFile(const std::string& zonfile, PARTICLE* p, pid_t np,
+                   ZONE *z, int numZones, int* zonenum, int *jumped);
+
+void writeVoidFile(const std::string& zonfile2, ZONE *z, int numZones);
+
+
+extern "C" void findrtop(double *a, int na, int *iord, int nb);
+
+#endif

c_source/jozov2/jozov2_io.cpp

@@ -0,0 +1,199 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_io.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <boost/format.hpp>
+#include "jozov2.hpp"
+#include "zobov.hpp"
+
+using namespace std;
+using boost::format;
+
+void readAdjacencyFile(const string& adjfile, PARTICLE*& p, pid_t& np)
+ throw(FileError)
+{
+  ifstream adj(adjfile.c_str());
+
+  if (!adj) {
+    cout << format("Unable to open %s") % adjfile << endl;
+    throw FileError();
+  }
+  adj.read((char *)&np, sizeof(pid_t));
+  
+  cout << format("adj: %d particles") % np << endl;
+  FF;
+
+  p = new PARTICLE[np];
+  
+  /* Adjacencies*/
+  for (pid_t i=0; i < np; i++) {
+    adj.read((char*)&p[i].nadj, sizeof(pid_t));
+    if (!adj)
+      throw FileError();
+
+    /* The number of adjacencies per particle */
+    if (p[i].nadj > 0)
+      p[i].adj = new pid_t[p[i].nadj];
+    else 
+      p[i].adj = 0;
+    p[i].ncnt = 0; /* Temporarily, it's an adj counter */
+  }
+
+  for (pid_t i = 0; i < np; i++)
+    {
+      pid_t nin;
+      adj.read((char*)&nin, sizeof(pid_t));
+      if (!adj)
+        throw FileError();
+
+      if (nin > 0)
+        {
+          for (int k=0; k < nin; k++) {
+            pid_t j;
+            
+            adj.read((char *)&j,sizeof(pid_t));
+            if (j < np)
+              {
+                /* Set both halves of the pair */
+                assert(i < j);
+                if (p[i].ncnt == p[i].nadj)
+                  {
+		    cout << format("OVERFLOW for particle %d (pending %d). List of accepted:") % i % j << endl;
+                    for (int q = 0; q < p[i].nadj; q++)
+                      cout << format("  %d") % p[i].adj[q] << endl;
+//                    throw FileError();
+                  }
+                else               
+                if (p[j].ncnt == p[j].nadj)
+                  {
+		    cout << format("OVERFLOW for particle %d (pending %d). List of accepted:") % j % i << endl;
+                    for (int q = 0; q < p[j].nadj; q++)
+                      cout << format("  %d\n") % p[j].adj[q] << endl;
+//                    throw FileError();
+                  }
+                else{
+                p[i].adj[p[i].ncnt] = j;
+                p[j].adj[p[j].ncnt] = i;
+                p[i].ncnt++; 
+                p[j].ncnt++; }
+              }
+            else
+              {
+                cout << format("%d: adj = %d") % i % j << endl;
+              }
+          }
+        }
+    }
+}
+
+void readVolumeFile(const std::string& volfile, PARTICLE *p, pid_t np, 
+                    pid_t mockIndex)
+  throw(FileError)
+{
+  ifstream vol(volfile.c_str());
+  pid_t np2;
+  
+  if (!vol)
+    {
+      cout << "Unable to open volume file " << volfile << endl;
+      throw FileError();
+    }
+  vol.read((char *)&np2, sizeof(pid_t));
+  if (np != np2) {
+    cout << format("Number of particles doesn't match! %d != %d") % np %np2 << endl;
+    throw FileError();
+  }
+  cout << format("%d particles") % np << endl;
+  FF;
+
+  for (pid_t i = 0; i < np; i++) {
+    vol.read((char*)&p[i].dens, sizeof(float));
+    if (((p[i].dens < 1e-30) || (p[i].dens > 1e30)) && (i < mockIndex)) {
+      cout << format("Whacked-out volume found, of particle %d: %f") % i % p[i].dens << endl;
+      p[i].dens = 1.;
+    }
+    p[i].dens = 1./p[i].dens; /* Get density from volume */
+  }
+}
+
+void writeZoneFile(const std::string& zonfile, PARTICLE* p, pid_t np,
+                   ZONE *z, int numZones, int* zonenum, int *jumped)
+{
+  pid_t **m = new pid_t *[numZones];
+  pid_t *nm = new pid_t[numZones];
+
+  /* Not in the zone struct since it'll be freed up (contiguously, we hope)
+     soon */
+  for (int h=0; h < numZones; h++)
+    {
+      m[h] = new pid_t[z[h].np];
+      nm[h] = 0;
+    }
+
+  for (pid_t i = 0; i < np; i++)
+    {
+      int h = zonenum[jumped[i]];
+      if (z[h].core == i)
+        {
+          m[h][nm[h]] = m[h][0];
+          m[h][0] = i; /* Put the core particle at the top of the list */
+        } 
+      else
+        {
+          m[h][nm[h]] = i;
+        }
+      nm[h] ++;
+    }
+  delete[] nm;
+
+  ofstream zon(zonfile.c_str());
+  if (!zon) 
+    {
+      cout << format("Problem opening zonefile %s.") % zonfile << endl;
+      throw FileError();
+    }
+  
+  zon.write((char *)&np, sizeof(pid_t));
+  zon.write((char *)&numZones, sizeof(int));
+  for (int h = 0; h < numZones; h++)
+    {
+      zon.write((char *)&(z[h].np), sizeof(pid_t));
+      zon.write((char *)m[h],z[h].np * sizeof(pid_t));
+      delete[] m[h];
+    }
+  delete[] m;
+}
+
+void writeVoidFile(const string& zonfile2, ZONE *z, int numZones)
+{
+  ofstream zon2(zonfile2.c_str());
+  if (!zon2)
+    {
+      cout << format("Problem opening zonefile %s.)") % zonfile2 << endl;
+      throw FileError();
+    }
+  zon2.write((char *)&numZones,sizeof(int));
+  cout << "Writing void/zone relations..." << endl;
+  for (int h = 0; h < numZones; h++)
+    {
+      zon2.write((char *)&z[h].nhl, sizeof(int));
+      zon2.write((char *)z[h].zonelist, z[h].nhl*sizeof(int));
+    }
+}

c_source/jozov2/jozov2_watershed.cpp

@@ -0,0 +1,294 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_watershed.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#ifdef OPENMP
+#include <omp.h>
+#endif
+
+#include <queue>
+#include <set>
+#include <iostream>
+#include <fstream>
+#include <boost/format.hpp>
+#include <string>
+#include "jozov2.hpp"
+#include "zobov.hpp"
+
+using namespace std;
+using boost::format;
+
+struct ZoneDensityPair
+{
+  int h;
+  double density;
+  double core;
+
+  bool operator<(const ZoneDensityPair& p2) const
+  {
+    return (density > p2.density) || (density==p2.density && core > p2.core);
+  }
+};
+
+typedef priority_queue<ZoneDensityPair> ZoneQueue;
+
+static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, char *inyet, int h)
+{
+  ZoneDensityPair zdp;
+  ZONE& z_h = z[h];
+  bool interior = true;
+
+  assert(inyet[h] == 1);
+  for (int za = 0; za < z_h.nadj; za++)
+   {
+     zdp.h = z_h.adj[za];
+     zdp.density = z_h.slv[za];
+     zdp.core = p[z[zdp.h].core].dens;
+     if (inyet[zdp.h] == 0)
+       {
+	 q.push(zdp);
+	 interior = false;
+       }
+   }
+  if (interior)
+    inyet[h] = 2;
+}
+
+void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol)
+{
+  /* Text output file */
+
+#pragma omp parallel
+  {
+    char *inyet, *inyet2;
+    int *zonelist, *zonelist2;
+    int nhl;
+    int prev_ii = -1;
+
+    inyet = new char[numZones];
+    inyet2 = new char[numZones];
+    zonelist = new int[numZones];
+    zonelist2 = new int[numZones];
+
+    fill(inyet, inyet + numZones, 0);
+    fill(inyet2, inyet2 + numZones, 0);
+
+    nhl = 0;
+#pragma omp for schedule(dynamic,1)
+    for (int h = 0; h < numZones; h++)
+      {
+        int nhlcount = 0;
+        float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens;
+        bool beaten;
+        priority_queue<ZoneDensityPair> to_process;
+	int link0;
+        int save_nhl;
+        pid_t save_npjoin;
+
+        for (int hl = 0; hl < nhl; hl++)
+          inyet[zonelist[hl]] = 0;
+        
+        zonelist[0] = h;
+        inyet[h] = 1;
+        save_nhl = nhl = 1;
+        save_npjoin = z[h].npjoin = z[h].np;
+        z_cur_core_dens = p[z[h].core].dens;
+
+        build_process_queue(to_process, z, p, inyet, h);
+	
+        do {
+          /* Find the lowest-volume (highest-density) adjacency */
+          beaten = false;
+          
+          if (to_process.empty())
+            {
+              beaten = true;
+              z[h].leak = maxvol;
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+              continue;
+            }
+
+	  do
+	    {
+	      lowvol = to_process.top().density;
+	      link0 = to_process.top().h; 
+	      to_process.pop();
+	    }
+	  while ((inyet[link0] != 0) && (!to_process.empty()));
+
+	  if (to_process.empty())
+	    {
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+	      beaten = true;
+	      z[h].leak = maxvol;
+	      continue;
+	    }
+
+          /* See if there's a beater */
+          if (previous_lowvol != lowvol)
+            {
+              save_npjoin = z[h].npjoin;
+              save_nhl = nhl;
+              previous_lowvol = lowvol;
+            }
+
+	  if (lowvol > voltol)
+            {
+              beaten = true;
+              z[h].leak = lowvol;
+              continue;
+            }
+          
+          if (p[z[link0].core].dens < z_cur_core_dens)
+            {
+	      beaten = true;
+              z[h].leak = lowvol;
+              continue;
+            }
+          
+          /* Add everything linked to the link(s) */
+          int nhl2 = 0;
+	  zonelist2[0] = link0;
+	  inyet2[link0] = 1;
+	  nhl2=1;
+	  
+	  bool added = true;
+	  while (added && !beaten)
+	    {
+	      added = false;
+	      for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
+		{
+		  int h2 = zonelist2[hl];
+                  
+		  if (inyet2[h2] == 1) {
+		    bool interior = true; /* Guilty until proven innocent */
+                    
+		    for (int za = 0; za < z[h2].nadj; za ++) {
+		      int link2 = z[h2].adj[za];
+		      
+		      if ((inyet[link2]+inyet2[link2]) == 0) {
+			interior = false;
+			if (z[h2].slv[za] <= lowvol) {
+			  if (p[z[link2].core].dens < z_cur_core_dens) {
+			    beaten = true;
+			    break;
+			  }
+			  zonelist2[nhl2] = link2;
+			  inyet2[link2] = 1;
+			  nhl2++;
+			  added = true;
+			}
+		      }
+		    }
+		    
+		    if (interior)
+		      inyet2[h2] = 2;
+		  }
+		}
+	    }
+	  
+	  /* See if there's a beater */
+	  if (beaten) {
+	    z[h].leak = lowvol;
+	  } else {
+
+	    for (int h2 = 0; h2 < nhl2; h2++) {
+	      int new_h = zonelist2[h2];
+	      
+              zonelist[nhl] = new_h;
+              assert(inyet[new_h] == 0);
+              z[h].npjoin += z[new_h].np;
+              inyet[new_h] = 1;
+              if (inyet2[new_h] != 2)
+                build_process_queue(to_process, z, p, inyet, new_h);
+              nhl++;
+            }
+          }
+	  
+          for (int hl = 0; hl < nhl2; hl++)
+            inyet2[zonelist2[hl]] = 0;
+          if (nhl/10000 > nhlcount) {
+            if (nhlcount == 0)
+              (cout << format("Zone %d: %d") % h % nhl).flush(); 
+            else
+              (cout << format(" %d [%d]") % nhl % to_process.size()).flush();
+            nhlcount = nhl/10000;
+          }
+	}
+	while((lowvol < BIGFLT) && (!beaten));
+
+        if (!beaten)
+          { 
+            save_npjoin = z[h].npjoin;
+            save_nhl = nhl;
+          }
+	
+	z[h].denscontrast = z[h].leak/p[z[h].core].dens;
+        if (z[h].denscontrast < 1.) 
+	  z[h].denscontrast = 1.;
+        
+        
+        /* Don't sort; want the core zone to be first */
+        
+        if (nhlcount > 0) { /* Outputs the number of zones in large voids */
+          printf(" h%d:%d\n",h,nhl);
+          FF;
+        }
+        /* Calculate volume */
+        z[h].npjoin = save_npjoin;
+        z[h].voljoin = 0.;
+        z[h].zonelist = new int[save_nhl];
+        z[h].numzones = save_nhl;
+        for (int q = 0; q < save_nhl; q++) {
+          z[h].voljoin += z[zonelist[q]].vol;
+          z[h].zonelist[q] = zonelist[q];
+        }
+        
+        z[h].nhl = save_nhl;
+      }
+    delete[] zonelist;
+    delete[] zonelist2;
+    delete[] inyet;
+    delete[] inyet2;
+    
+  }
+  
+  double maxdenscontrast = 0;
+#pragma omp parallel shared(maxdenscontrast)
+  {
+    double maxdenscontrast_local = 0;
+
+#pragma omp for schedule(static)
+    for (int h = 0; h < numZones; h++)
+      {
+        /* find biggest denscontrast */
+        if (z[h].denscontrast > maxdenscontrast_local) {
+          maxdenscontrast_local = (double)z[h].denscontrast;
+        }
+      }
+
+#pragma omp critical
+    {
+      if (maxdenscontrast_local > maxdenscontrast)
+        maxdenscontrast = maxdenscontrast_local;
+    }
+  }
+
+  cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl;
+}

c_source/jozov2/jozov2_zones.cpp

@@ -0,0 +1,243 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/jozov2_zones.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#include <iostream>
+#include <fstream>
+#include "jozov2.hpp"
+#include "zobov.hpp"
+#include <boost/format.hpp>
+
+using namespace std;
+using boost::format;
+
+void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped, 
+                      pid_t *numinh, pid_t& numZones)
+{  
+  pid_t *jumper = new pid_t[np];
+  float minvol;
+
+  /* find jumper */
+  for (pid_t i = 0; i < np; i++) {
+    PARTICLE& cur_p = p[i];
+
+    minvol = cur_p.dens;
+    jumper[i] = -1;
+    for (int j = 0; j < cur_p.nadj; j++) {
+      if (p[cur_p.adj[j]].dens < minvol) {
+	jumper[i] = cur_p.adj[j];
+	minvol = p[jumper[i]].dens;
+      }
+    }
+    numinh[i] = 0;
+  }
+
+  (cout << "About to jump ... " << endl).flush();
+  
+  /* Jump */
+  for (pid_t i = 0; i < np; i++) {
+    jumped[i] = i;
+    while (jumper[jumped[i]] > -1)
+      jumped[i] = jumper[jumped[i]];
+    numinh[jumped[i]]++;
+  }
+  (cout << "Post-jump ..." << endl).flush();
+  
+  pid_t loc_NumZones = 0;
+#pragma omp parallel for schedule(static) reduction(+:loc_NumZones)
+  for (pid_t i = 0; i < np; i++)
+    if (numinh[i] > 0) 
+      loc_NumZones++;
+
+  numZones = loc_NumZones;
+  cout << format("%d initial zones found") % numZones << endl;
+
+  delete[] jumper;
+}
+
+void buildZoneAdjacencies(PARTICLE *p, pid_t np,
+                          ZONE *z, ZONET *zt,
+                          int numZones,
+                          pid_t *jumped,
+                          int *zonenum,
+                          int *numinh)
+{
+  /* Count border particles */
+  for (pid_t i = 0; i < np; i++)
+    for (int j = 0; j < p[i].nadj; j++) {
+      pid_t testpart = p[i].adj[j];
+      if (jumped[i] != jumped[testpart])
+	zt[zonenum[jumped[i]]].nadj++;
+    }
+
+  size_t nadjAlloced = 0;
+  cout << "Total numZones = " << numZones << endl;
+  try
+    {
+      for (int h = 0; h < numZones; h++) {
+        if (zt[h].nadj > 0) {
+          zt[h].adj = new pid_t[zt[h].nadj];
+          zt[h].slv = new float[zt[h].nadj];
+          nadjAlloced += zt[h].nadj;
+          zt[h].nadj = 0;
+        } else {
+          zt[h].adj = 0;
+          zt[h].slv = 0;
+        }
+      }
+    }
+  catch(const std::bad_alloc& a)
+    {
+      cout << "Could not allocate memory for zone adjacencies (nadj so far: " << nadjAlloced << ", memory needed: " << (nadjAlloced*(sizeof(pid_t)+sizeof(float))) << ")" << endl;
+      throw a;
+    }
+
+  /* Find "weakest links" */
+  for (pid_t i = 0; i < np; i++)
+    {
+      int h = zonenum[jumped[i]];
+      ZONET& zt_h = zt[h];
+    
+      for (int j = 0; j < p[i].nadj; j++)
+        {
+          pid_t testpart = p[i].adj[j];
+          bool already;
+
+          if (h != zonenum[jumped[testpart]]) {
+            if (p[testpart].dens > p[i].dens) {
+              /* there could be a weakest link through testpart */
+              already = false;
+              for (int za = 0; za < zt_h.nadj; za++)
+                if (zt_h.adj[za] == zonenum[jumped[testpart]]) {
+                  already = true;
+                  if (p[testpart].dens < zt_h.slv[za]) {
+                    zt_h.slv[za] = p[testpart].dens;
+                  }
+                }
+              if (!already) {
+                zt_h.adj[zt_h.nadj] = zonenum[jumped[testpart]];
+                zt_h.slv[zt_h.nadj] = p[testpart].dens;
+                zt_h.nadj++;
+              }
+            } else { /* There could be a weakest link through i */
+              already = false;
+              for (int za = 0; za < zt_h.nadj; za++)
+                if (zt_h.adj[za] == zonenum[jumped[testpart]]) {
+                  already = true;
+                  if (p[i].dens < zt_h.slv[za]) {
+                    zt_h.slv[za] = p[i].dens;
+                  }
+                }
+              if (!already) {
+                zt_h.adj[zt_h.nadj] = zonenum[jumped[testpart]];
+                zt_h.slv[zt_h.nadj] = p[i].dens;
+                zt_h.nadj++;
+              }
+            }
+          }
+        }
+    }
+  (cout <<"Found zone adjacencies" << endl).flush();
+
+  /* Free particle adjacencies */
+  for (pid_t i = 0; i < np; i++)
+    {
+      if (p[i].adj != 0)
+        delete[] p[i].adj;
+    }
+
+  /* Use z instead of zt */
+  for (int h = 0; h < numZones; h++) {
+    z[h].nadj = zt[h].nadj;
+    if (zt[h].nadj > 0) {
+      z[h].adj = new pid_t[zt[h].nadj];
+      z[h].slv = new float[zt[h].nadj];
+      for (int za = 0; za < zt[h].nadj; za++) {
+        z[h].adj[za] = zt[h].adj[za];
+        z[h].slv[za] = zt[h].slv[za];
+      }
+      delete[] zt[h].adj;
+      delete[] zt[h].slv;
+    } else {
+      z[h].adj = 0;
+      z[h].slv = 0;
+    }
+    z[h].np = numinh[z[h].core];
+  }
+}
+
+void buildZones(PARTICLE *p, pid_t np, pid_t *&jumped, 
+                ZONE*& z, int& nzones,
+                int*& zonenum)
+{
+  pid_t *numinh;
+  ZONET *zt;
+
+  jumped = new pid_t[np];
+  numinh = new pid_t[np];
+
+  buildInitialZones(p, np, jumped, numinh, nzones);
+
+  try
+    {
+      z = new ZONE[nzones];
+      zt = new ZONET[nzones];
+      zonenum = new int[np];
+    }
+  catch (const std::bad_alloc& e)
+    {
+      cout << "Unable to do zone allocations" << endl;
+      throw e;
+    }
+
+  int h = 0;
+  for (pid_t i = 0; i < np; i++)
+    {
+      if (numinh[i] > 0) {
+        z[h].core = i;
+        z[h].vol = 0;
+        z[h].np = z[h].npjoin = z[h].nadj = z[h].nhl = 0;
+        z[h].leak = 0;
+        z[h].adj = 0;
+        z[h].slv = 0;
+        z[h].denscontrast = 0;
+        z[h].vol = z[h].voljoin = 0;
+        zonenum[i] = h;
+        h++;
+      } else {
+        zonenum[i] = -1;
+      }
+    }
+  for (size_t z = 0; z < nzones; z++) {
+    zt[z].nadj = 0;
+    zt[z].adj = 0;
+    zt[z].slv = 0;
+  }
+
+  buildZoneAdjacencies(p, np, z, zt,
+                       h, jumped, zonenum, numinh);
+
+  delete[] zt;
+  delete[] numinh;
+
+  for (pid_t i=0; i < np; i++) {
+    int h = zonenum[jumped[i]];
+    z[h].vol += 1.0/(double)p[i].dens;
+    z[h].numzones = 0;
+    z[h].zonelist = 0;
+  }
+}

c_source/jozov2/zobov.hpp

@@ -0,0 +1,52 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/zobov2/jozov2/zobov.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+#ifndef __ZOBOV_HPP
+#define __ZOBOV_HPP
+
+typedef struct Particle {
+  float dens;
+  int nadj;
+  int ncnt;
+  int *adj;
+} PARTICLE;
+
+typedef struct Zone {
+  int core; /* Identity of peak particle */
+  int np; /* Number of particles in zone */
+  int npjoin; /* Number of particles in the joined void */
+  int nadj; /* Number of adjacent zones */
+  int nhl; /* Number of zones in final joined void */
+  float leak; /* Volume of last leak zone*/
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+  float denscontrast; /* density contrast */
+  double vol; /* Total volume of all particles in the zone */
+  double voljoin; /* Total volume of all particles in the joined void */
+
+  int *zonelist; /* Zones bound to the void. */
+  int numzones; /* Number of zones bound. */
+} ZONE;
+
+typedef struct ZoneT {
+  int nadj; /* Number of zones on border */
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+} ZONET;
+
+#endif

c_source/prep/CMakeLists.txt

@@ -0,0 +1,22 @@
+include_directories(${CMAKE_CURRENT_BINARY_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/loaders)
+
+IF(SDF_SUPPORT)
+  add_definitions(-DSDF_SUPPORT)
+ENDIF(SDF_SUPPORT)
+
+SET(prepSimulation_SRCS prepSimulation.cpp )
+add_genopt(prepSimulation_SRCS prepSimulation.ggo prepSimulation_conf STRUCTNAME prepSimulation_info)
+add_executable(prepSimulation ${prepSimulation_SRCS})
+target_link_libraries(prepSimulation simu_loaders ${ZOB_LIBS} ${LIBSDF_LIBRARY} ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES} ${HDF5_LIBRARIES} ${DL_LIBRARY})
+
+SET(prepObservation_SRCS prepObservation.cpp)
+add_genopt(prepObservation_SRCS prepObservation.ggo prepObservation_conf STRUCTNAME prepObservation_info)
+add_executable(prepObservation ${prepObservation_SRCS})
+target_link_libraries(prepObservation ${ZOB_LIBS} ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES} ${HDF5_LIBRARIES} ${HEALPIX_LIBRARIES} ${DL_LIBRARY})
+
+set_target_properties(prepObservation PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
+
+subdirs(loaders)
+
+install(TARGETS prepSimulation prepObservation DESTINATION ${VIDE_BIN})
+

c_source/prep/generateMock.cpp

@@ -0,0 +1,68 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/generateTestMock.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <iostream>
+#include <cstdlib>
+#include <CosmoTool/fortran.hpp>
+
+using namespace CosmoTool;
+using namespace std;
+
+#define LX 1.0
+#define LY 1.0
+#define LZ 1.0
+
+#define NUMPART (16*16*16)
+
+int main(int argc, char **argv)
+{
+  UnformattedWrite f("particles.bin");
+  
+  f.beginCheckpoint();
+  f.writeInt32(NUMPART);
+  f.endCheckpoint();
+
+  cout << "Writing X components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LX);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing Y components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LY);
+    }
+  f.endCheckpoint();
+  
+  cout << "Writing Z components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < NUMPART; i++)
+    {
+      f.writeReal32(drand48()*LZ);
+    }
+  f.endCheckpoint();
+
+  return 0;
+}

c_source/prep/loaders/CMakeLists.txt

@@ -0,0 +1,11 @@
+SET(SIMU_LOADERS_SRCS 
+   simulation_loader.cpp ramses_loader.cpp 
+   gadget_loader.cpp flash_loader.cpp)
+
+IF (SDF_SUPPORT)
+ add_definitions(-DSDF_SUPPORT)
+ SET(SIMU_LOADERS_SRCS ${SIMU_LOADERS_SRCS} sdf_loader.cpp multidark_loader.cpp)
+ENDIF (SDF_SUPPORT)
+
+
+add_library(simu_loaders ${SIMU_LOADERS_SRCS})

c_source/prep/loaders/basic_loader.cpp

@@ -0,0 +1,149 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/basic_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <iostream>
+#include <boost/format.hpp>
+#include <fstream>
+#include <CosmoTool/yorick.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+using boost::format;
+using boost::str;
+
+class BasicGroupLoader: public SimulationLoader
+{
+private:
+  string storage;
+  SimuData *header;
+  int flags, numFiles;
+public:
+  BasicGroupLoader(const string& storage_path, SimuData *header, int flags, int numfiles)
+  {
+    this->header = header;
+    this->storage = storage_path;
+    this->flags = flags;
+    this->numFiles = numfiles;
+  }
+  
+  SimuData *getHeader() {
+    return header;
+  }
+  
+  int num_files() {
+    return numFiles;
+  }
+  
+  SimuData *loadFile(int id) {    
+    if (id < 0 || id >= numFiles)
+      return 0;
+
+    SimuData *simu = new SimuData;
+    uint32_t *dimlist, dimrank;
+    string fname;
+
+    simu->time = header->time;
+    simu->TotalNumPart = header->NumPart;
+    simu->Omega_M = header->Omega_M;
+    simu->Omega_Lambda = header->Omega_Lambda;
+    simu->NumPart = -1;
+
+    if (flags & NEED_POSITION)
+      {
+	loadArray(str(format("%s/x_%d.nc") % storage % id), 
+		  simu->Pos[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	simu->NumPart = dimlist[0];
+	
+	loadArray(str(format("%s/y_%d.nc") % storage % id), 
+		  simu->Pos[1], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/z_%d.nc") % storage % id), 
+		  simu->Pos[2], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+      }
+
+    if (flags & NEED_VELOCITY)
+      {
+	loadArray(str(format("%s/vx_%d.nc") % storage % id), 
+		  simu->Vel[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	if (simu->NumPart < 0)
+	  simu->NumPart = dimlist[0];
+
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/vy_%d.nc") % storage % id), 
+		  simu->Vel[0], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+	
+	loadArray(str(format("%s/vz_%d.nc") % storage % id), 
+		  simu->Vel[2], dimlist, dimrank);
+	assert(dimrank == 1);
+	assert(simu->NumPart == dimlist[0]);
+      }
+
+    if (flags & NEED_GADGET_ID)
+      {
+	loadArray(str(format("%s/id_%d.nc") % storage % id), 
+		  simu->Id, dimlist, dimrank);
+	assert(dimrank == 1);
+	if (simu->NumPart < 0)
+	  simu->NumPart = dimlist[0];
+
+	assert(simu->NumPart == dimlist[0]);
+      }
+    
+    return simu;
+  }
+
+  ~BasicGroupLoader()
+  {
+    delete header;
+  }
+};
+
+SimulationLoader *basicGroupLoader(const std::string& simupath, 
+				   int flags)
+{
+  SimuData *header;
+  ifstream f;
+  string header_path = simupath + "/header.txt";
+  int numFiles;
+
+  header = new SimuData;
+  f.open(header_path.c_str());
+
+  f >> header->time
+    >> header->Omega_M
+    >> header->Omega_Lambda 
+    >> header->NumPart
+    >> numFiles;
+
+  return new BasicGroupLoader(simupath, header, flags, numFiles);
+}

c_source/prep/loaders/flash_loader.cpp

@@ -0,0 +1,134 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/flash_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadFlash.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class FlashLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  SimuData *gadget_header;
+  string snapshot_name;
+  SimulationPreprocessor *preproc;
+public:
+  FlashLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), gadget_header(header), preproc(p)
+  {
+  }
+  
+  ~FlashLoader()
+  {
+    delete gadget_header;
+  }
+  
+  SimuData *getHeader() {
+    return gadget_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+    if (onefile && id > 0)
+      return 0;
+    if (id >= _num_files)
+      return 0;
+
+    if (onefile)
+      d = loadFlashMulti(snapshot_name.c_str(), -1, load_flags);
+    else
+      d = loadFlashMulti(snapshot_name.c_str(), id, load_flags);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    applyTransformations(d);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *flashLoader(const std::string& snapshot, int flags, SimulationPreprocessor *p)
+{
+  bool singleFile;
+  int num_files;
+  SimuData *d;
+
+  try
+    {
+      d = loadFlashMulti(snapshot.c_str(), -1, 0);
+      singleFile = true;
+      num_files = 1;
+    }
+  catch (const NoSuchFileException& e)
+    {
+      try
+        {
+          d = loadFlashMulti(snapshot.c_str(), 0, 0);
+	  num_files = 0;
+        }
+      catch(const NoSuchFileException& e)
+        {
+          return 0;
+        }
+    }
+    
+  assert(d != 0);
+  SimuData *header = d;
+    
+  if (!singleFile)
+    {
+      try
+	{
+	  while ((d = loadFlashMulti(snapshot.c_str(), num_files, 0)) != 0)
+	    {
+	      num_files++;
+	      delete d;
+	    }
+	}
+      catch(const NoSuchFileException& e)
+	{
+	}
+    }
+    
+  return new FlashLoader(snapshot, header, flags, singleFile, num_files, p);
+}

c_source/prep/loaders/gadget_loader.cpp

@@ -0,0 +1,151 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/gadget_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <vector>
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadGadget.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+
+class GadgetLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  double unitMpc;
+  int gadgetFormat;
+  SimuData *gadget_header;
+  string snapshot_name;
+  SimulationPreprocessor *preproc;
+public:
+  GadgetLoader(const string& basename, SimuData *header, int flags, bool singleFile, int _num, double unit, int gadgetFormat, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), onefile(singleFile), _num_files(_num), unitMpc(1/unit), gadget_header(header), gadgetFormat(gadgetFormat), preproc(p)
+  {
+  }
+  
+  ~GadgetLoader()
+  {
+    delete gadget_header;
+  }
+  
+  SimuData *getHeader() {
+    return gadget_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+    if (onefile && id > 0)
+      return 0;
+    if (id >= _num_files)
+      return 0;
+
+    if (onefile)
+      d = loadGadgetMulti(snapshot_name.c_str(), -1, load_flags, gadgetFormat);
+    else
+      d = loadGadgetMulti(snapshot_name.c_str(), id, load_flags, gadgetFormat);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    for (int k = 0; k < 3; k++)
+      {
+        if (d->Pos[k] != 0)
+          {
+            for (long i = 0; i < d->NumPart; i++)
+              d->Pos[k][i] *= unitMpc;
+          }
+      }
+    d->BoxSize *= unitMpc;
+
+    applyTransformations(d);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+
+SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags, int gadgetFormat, SimulationPreprocessor *p)
+{
+  bool singleFile = false;
+  int num_files;
+  SimuData *d;
+
+  cout << " File to load is:" << snapshot.c_str() << endl;
+
+  try
+    {
+      d = loadGadgetMulti(snapshot.c_str(), -1, 0, gadgetFormat);
+      singleFile = true;
+      num_files = 1;
+    }
+  catch (const NoSuchFileException& e)
+    {
+      try
+        {
+          d = loadGadgetMulti(snapshot.c_str(), 0, 0, gadgetFormat);
+	  num_files = 0;
+        }
+      catch(const NoSuchFileException& e)
+        {
+          return 0;
+        }
+    }
+    
+  assert(d != 0);
+  SimuData *header = d;
+
+  header->BoxSize /= Mpc_unitLength;
+    
+  if (!singleFile)
+    {
+      try
+	{
+	  while ((d = loadGadgetMulti(snapshot.c_str(), num_files, 0, gadgetFormat)) != 0)
+	    {
+	      num_files++;
+	      delete d;
+	    }
+	}
+      catch(const NoSuchFileException& e)
+	{
+	}
+    }
+    
+  return new GadgetLoader(snapshot, header, flags, singleFile, num_files, Mpc_unitLength, gadgetFormat, p);
+}

c_source/prep/loaders/multidark_loader.cpp

@@ -0,0 +1,177 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/multidark_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <boost/format.hpp>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+using boost::format;
+using namespace std;
+using namespace CosmoTool;
+
+static const double one_kpc = 3.08567802e16; /* km */
+static const double one_Gyr = 3.1558149984e16; /* sec */
+
+class MultiDarkLoader: public SimulationLoader
+{
+protected:
+  SimuData *header;
+  string darkname;
+  SimulationPreprocessor *preproc;
+public:
+  MultiDarkLoader(const std::string& name, SimuData *h, SimulationPreprocessor *p)
+    : preproc(p), darkname(name), header(h)
+  {
+  }
+
+  ~MultiDarkLoader()
+  {
+    delete header;
+  }
+
+  int num_files()
+  {
+    return 1;
+  }
+
+  SimuData *getHeader()
+  {
+    return header;
+  }
+
+  SimuData *loadFile(int id)
+  {
+    if (id != 0)
+      return 0;
+
+    ifstream fp(darkname.c_str());
+    SimuData *simu = new SimuData;
+
+    fp >> simu->BoxSize >> simu->Omega_M >> simu->Hubble >> simu->time >> simu->NumPart;
+    simu->time = 1./(1.+simu->time); // convert to scale factor
+    simu->TotalNumPart = simu->NumPart;
+    simu->Omega_Lambda = 1.0 - simu->Omega_M;
+
+    long estimated = (preproc == 0) ? simu->NumPart : preproc->getEstimatedPostprocessed(simu->NumPart);
+    long allocated = estimated; 
+
+    for (int k = 0; k < 3; k++) {
+      simu->Pos[k] = new float[allocated];
+      simu->Vel[k] = new float[allocated];
+    }
+    simu->Id = new int64_t[allocated];
+    long *uniqueID = new long[allocated];
+    long *index = new long[allocated];
+
+    double tempData;
+
+    double shift = 0.5*simu->BoxSize;
+    double rescale_position = simu->Hubble*1.e-5*100./simu->time;
+    double rescale_velocity = one_kpc/one_Gyr;
+#define INFINITY std::numeric_limits<float>::max()
+    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
+
+    cout << "loading multidark particles" << endl;
+    long actualNumPart = 0;
+
+    for (long i = 0; i < simu->NumPart; i++) {
+      SingleParticle p;
+
+      fp >> p.ID >> p.Pos[0] >> p.Pos[1]
+         >> p.Pos[2] >> p.Vel[2] >> p.Vel[1] >> p.Vel[0] >> tempData;
+
+      if (p.ID == -99 && 
+          p.Pos[0] == -99 && p.Pos[1] == -99 && 
+          p.Pos[2] == -99 && p.Vel[2] == -99)
+        break;
+
+      //p.Pos[0] = p.Pos[0]*rescale_position + shift;
+      //p.Pos[1] = p.Pos[1]*rescale_position + shift;
+      //p.Pos[2] = p.Pos[2]*rescale_position + shift;
+      //p.Vel[2] = p.Vel[2]*rescale_velocity;
+
+      // enforce box size in case of roundoff error
+      for (int k = 0; k < 3; k++) {
+        if (p.Pos[k] < 0) p.Pos[k] += simu->BoxSize;
+        if (p.Pos[k] >= simu->BoxSize) p.Pos[k] -= simu->BoxSize;
+      }
+
+      if (preproc != 0 && !preproc->accept(p))
+        continue;
+
+      copyParticleToSimu(p, simu, actualNumPart);
+      uniqueID[actualNumPart]= p.ID;
+      index[actualNumPart] = i;
+      actualNumPart++;
+      if (actualNumPart == allocated)
+        {
+          allocated += (estimated+9)/10;
+          reallocSimu(simu, allocated); 
+          reallocArray(uniqueID, allocated, actualNumPart);
+          reallocArray(index, allocated, actualNumPart);
+        }
+
+       for (int k = 0; k < 3; k++) {
+         min_pos[k] = std::min(min_pos[k], p.Pos[k]);
+         max_pos[k] = std::max(max_pos[k], p.Pos[k]);
+       }
+    }
+    for (int k = 0; k < 3; k++) cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
+
+    applyTransformations(simu);
+    simu->NumPart = actualNumPart;
+    simu->TotalNumPart = actualNumPart;
+    simu->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+    simu->new_attribute("index", index, delete_adaptor<long>);
+    return simu;
+  }
+};
+
+SimulationLoader *multidarkLoader(const string& multidarkname, SimulationPreprocessor *p)
+{
+  SimuData *header;
+  int actualNumPart;
+  ifstream fp(multidarkname.c_str());
+
+  cout << "opening multidark file " << multidarkname << endl;
+  if (!fp)
+    {
+      cout << "could not open file!" << endl;
+      return 0;
+    }
+
+  header = new SimuData();
+  fp >> header->BoxSize >> header->Omega_M >> header->Hubble >> header->time >> header->NumPart;
+  
+  header->time = 1./(1.+header->time); // convert to scale factor
+  header->TotalNumPart = header->NumPart;
+  header->Omega_Lambda = 1.0 - header->Omega_M;
+
+  return new MultiDarkLoader(multidarkname, header, p);
+}
+
+  
+

c_source/prep/loaders/ramses_loader.cpp

@@ -0,0 +1,183 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/ramses_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <CosmoTool/loadRamses.hpp>
+#include <CosmoTool/fortran.hpp>
+#include "simulation_loader.hpp"
+
+// ben edit
+
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <sys/types.h>
+#include <regex.h>
+#include <cstring>
+#include <cstdlib>
+#include <cassert>
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <map>
+
+
+// end ben edit
+
+
+using namespace std;
+using namespace CosmoTool;
+
+class RamsesLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  int _num_files;
+  int baseid;
+  bool double_precision;
+  SimuData *ramses_header;
+  string snapshot_name;
+  SimulationPreprocessor *preproc;
+public:
+  RamsesLoader(const string& basename, int baseid, bool dp, SimuData *header, int flags, int _num, SimulationPreprocessor *p)
+    : snapshot_name(basename), load_flags(flags), _num_files(_num), double_precision(dp),
+      ramses_header(header), preproc(p)
+  {
+  }
+  
+  ~RamsesLoader()
+  {
+    delete ramses_header;
+  }
+  
+  SimuData *getHeader() {
+    return ramses_header;
+  }
+  
+  int num_files() {
+    return _num_files;
+  }
+  
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    
+
+// cludgy hack until I can get baseid working in this function... the problem is with   SimuData *loadFile(int id) ... just need to learn a bit more C++ ~ Ben
+
+    string baseidstr = snapshot_name.c_str();
+    unsigned found = baseidstr.find_last_of("/");
+    baseidstr = baseidstr.substr(found-5,found);
+    baseidstr = baseidstr.substr(0,5);  // remove trailing slash
+    baseid = atoi(baseidstr.c_str());
+    
+    if (id >= _num_files)
+      return 0;
+
+    d = loadRamsesSimu(snapshot_name.c_str(), baseid, id, double_precision, load_flags);
+    assert(d != 0);
+
+    if (d->Id != 0)
+      {
+	long *uniqueID = new long[d->NumPart];
+	for (long i = 0; i < d->NumPart; i++)
+	  {
+	    uniqueID[i] = d->Id[i];
+	  }
+	d->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+      }
+
+    applyTransformations(d);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags, SimulationPreprocessor *p)
+{
+  SimuData *d, *header;
+  int num_files = 0; // how many particle files are there?
+
+  header = loadRamsesSimu(snapshot.c_str(), baseid, 0, double_precision, 0);
+
+  if (header == 0)
+    return 0;
+
+
+// count number of CPU's for this output... kinda copy/paste of loadRamses.cpp in cosmotool/src
+
+  ostringstream ss_fname;  
+  ss_fname << snapshot.c_str() << "/info_" << setfill('0') << setw(5) << baseid << ".txt";
+  cout << "Opening info file " << ss_fname.str() << " to find cpu number" << endl;
+  ifstream infile(ss_fname.str().c_str());
+  if (!infile)
+    return 0;
+
+
+  int err;
+  regex_t unit_l_rx;
+
+  //  const char *pattern = "^unit_l[ ]*=[ ]*([0-9\\.E+\\-]+)";
+  const char *pattern = "^([A-Za-z_]+)[ ]*=[ ]*([0-9\\.E+\\-]+)";
+
+  err = regcomp (&unit_l_rx, pattern, REG_EXTENDED);
+  cout << unit_l_rx.re_nsub << endl;
+  if (err)
+    {
+      char errString[255];
+      regerror(err, &unit_l_rx, errString, sizeof(errString));
+      cout << errString << endl;
+      abort();
+    }
+
+  map<string,double> infoMap;
+  string line;
+  while (getline(infile, line))
+    {
+      regmatch_t allMatch[4];
+      if (!regexec(&unit_l_rx, line.c_str(), 4, allMatch, 0))
+	{
+	  uint32_t start0 = allMatch[1].rm_so, end0 = allMatch[1].rm_eo;
+	  uint32_t start1 = allMatch[2].rm_so, end1 = allMatch[2].rm_eo;
+
+	  string keyword = line.substr(start0, end0-start0);
+	  istringstream iss(line.substr(start1, end1-start1));
+	  double unitLength;
+
+	  iss >> unitLength;
+
+	  infoMap[keyword] = unitLength;
+	}
+    }
+  
+  regfree(&unit_l_rx);
+
+
+  num_files = infoMap["ncpu"];
+
+
+  return new RamsesLoader(snapshot, baseid, double_precision, header, flags, num_files, p);
+}
+

c_source/prep/loaders/sdf_loader.cpp

@@ -0,0 +1,308 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/sdf_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <iostream>
+#include <boost/format.hpp>
+#include <vector>
+#include <cassert>
+#include <string>
+#include "sdfloader_internal.hpp"
+#include "simulation_loader.hpp"
+#include <libSDF/SDF.h>
+
+#undef SDFgetfloatOrDie
+#undef SDFgetdoubleOrDie
+#undef SDFgetintOrDie
+
+using boost::format;
+
+using namespace std;
+using namespace CosmoTool;
+
+static const double one_kpc = 3.08567802e16; /* km */
+static const double one_Gyr = 3.1558149984e16; /* sec */
+
+static void SDFgetfloatOrDie(SDF *sdfp, const char *name, float *v)
+{
+   if( SDFgetfloat(sdfp, name, v) ) {
+     cerr << format("SDFgetfloat(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+static void SDFgetdoubleOrDie(SDF *sdfp, const char *name, double *v)
+{
+   if( SDFgetdouble(sdfp, name, v) ) {
+     cerr << format("SDFgetdouble(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+static void SDFgetintOrDie(SDF *sdfp, const char *name, int *v)
+{
+   if( SDFgetint(sdfp, name, v) ) {
+     cerr << format("SDFgetint(%s) failed") % name << endl;
+     abort();
+   }
+}
+
+class SDFLoader: public SimulationLoader
+{
+private:
+  int load_flags;
+  bool onefile;
+  int _num_files;
+  double unitMpc;
+  SimuData *sdf_header;
+  SDF *sdfp;
+  SimulationPreprocessor *preproc;
+  int num_splitting;
+public:
+  SDFLoader(SDF *sdf_file, SimuData *header, int flags, int num_splitting,
+            SimulationPreprocessor *p)
+    : sdfp(sdf_file), load_flags(flags),
+      sdf_header(header), preproc(p)
+  {
+    this->num_splitting = num_splitting;
+  }
+  
+  ~SDFLoader()
+  {
+    delete sdf_header;
+  }
+  
+  SimuData *getHeader() {
+    return sdf_header;
+  }
+  
+  int num_files() {
+    return num_splitting;
+  }
+
+  int64_t getStart(int id)
+  {
+    return sdf_header->TotalNumPart * int64_t(id) / num_splitting;
+  }
+
+  int64_t getNumberInSplit(int id)
+  {
+    return getStart(id+1)-getStart(id);
+  }
+
+  void rescaleParticles(SimuData *d,
+                        double rescale_position,
+                        double rescale_velocity)
+  {
+#define INFINITY std::numeric_limits<float>::max()
+    float shift = 0.5*d->BoxSize;
+    rescale_position /= d->time;
+    float min_pos[3] = {INFINITY,INFINITY, INFINITY}, max_pos[3] = {-INFINITY,-INFINITY,-INFINITY};
+
+    if (d->Pos[0] != 0)
+      {
+        for (int k = 0; k < 3; k++)
+          {
+            for (int64_t i = 0; i < d->NumPart; i++)
+              {
+                d->Pos[k][i] = (d->Pos[k][i]*rescale_position + shift);
+                min_pos[k] = std::min(min_pos[k], d->Pos[k][i]);
+                max_pos[k] = std::max(max_pos[k], d->Pos[k][i]);
+              }
+            cout << boost::format("min[%d] = %g, max[%d] = %g") % k % min_pos[k] % k %max_pos[k] << endl;
+          }
+      }
+    if (d->Vel[0] != 0)
+      {
+        for (int k = 0; k < 3; k++)
+          {
+            for (int64_t i = 0; i < d->NumPart; i++)
+              {
+                d->Vel[k][i] *= rescale_velocity;
+              }
+          }
+      }
+  }
+
+
+  void enforceBoxSize(SimuData *d)
+  {
+    for (int k = 0; k < 3; k++)
+      {
+        for (int64_t i = 0; i < d->NumPart; i++)
+         {
+           while (d->Pos[k][i]<0)
+             d->Pos[k][i] += d->BoxSize;
+           while (d->Pos[k][i]>=d->BoxSize)
+             d->Pos[k][i] -= d->BoxSize;
+         }
+       }
+  }
+
+  SimuData *loadFile(int id) {
+    SimuData *d;
+    int64_t starts[7];
+    int nobjs[7];
+    int strides[7];
+    void *addrs[7];
+    const char *names[7];
+    int ret;
+
+    if (id >= num_splitting)
+      return 0;
+
+    d = new SimuData;
+    d->BoxSize = sdf_header->BoxSize;
+    d->time = sdf_header->time;
+    d->Hubble = sdf_header->Hubble;
+    d->Omega_M = sdf_header->Omega_M;
+    d->Omega_Lambda = sdf_header->Omega_Lambda;
+    
+    d->NumPart = getNumberInSplit(id);
+
+    int64_t numPartToLoad = getNumberInSplit(id);
+    int64_t start = getStart(id);
+    int k = 0;
+#define SETUP_READ(name, vec) \
+    names[k] = name, starts[k] = start, nobjs[k] = numPartToLoad, strides[k] = sizeof(vec[0]), addrs[k] = vec, k++;
+
+    if (load_flags & NEED_POSITION)
+      {
+        const char *name_template[3] = { "x", "y", "z" };
+        for (int q = 0; q < 3; q++)
+          {
+            d->Pos[q] = new float[numPartToLoad];
+            SETUP_READ(name_template[q], d->Pos[q]);
+          }
+      }
+    if (load_flags & NEED_VELOCITY)
+      {
+        const char *name_template[3] = { "vx", "vy", "vz" };
+        for (int q = 0; q < 3; q++)
+          {
+            d->Vel[q] = new float[numPartToLoad];
+            SETUP_READ(name_template[q], d->Vel[q]);
+          }
+      }
+    if (load_flags & NEED_GADGET_ID)
+      {
+        d->Id = new int64_t[numPartToLoad];
+        SETUP_READ("ident", d->Id);
+      }
+#undef SETUP_READ
+
+    ret = SDFseekrdvecsarr(sdfp, k, (char**)names, starts, nobjs, addrs, strides);
+    if (ret != 0)
+      {
+        cerr << format("Splitting %d/%d, SDF read failure: %s") % id % num_splitting % SDFerrstring << endl;
+        abort();
+      }
+
+    if (load_flags & (NEED_POSITION | NEED_VELOCITY))
+      rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr);
+
+    enforceBoxSize(d);
+    applyTransformations(d);
+    basicPreprocessing(d, preproc);
+
+    return d;
+  }
+};
+
+SimulationLoader *sdfLoader(const std::string& snapshot, int flags, 
+                            int num_splitting,
+                            SimulationPreprocessor *p)
+{
+  SDF *sdfp;
+  int fileversion;
+  SimuData *hdr;
+  int64_t gnobj;
+
+  sdfp = SDFopen(0, snapshot.c_str());
+  if (sdfp == 0)
+    {
+      return 0;
+    }
+   cout << "Loading SDF with artificial splitting " << num_splitting << endl;
+  hdr = new SimuData;
+  
+  SDFgetintOrDefault(sdfp, "version", &fileversion, 1);
+  double h0;
+  if (fileversion == 1)
+    {
+      SDFgetfloatOrDie(sdfp, "Omega0", &hdr->Omega_M);
+      SDFgetfloatOrDie(sdfp, "Lambda_prime", &hdr->Omega_Lambda);
+      SDFgetfloatOrDie(sdfp, "H0", &hdr->Hubble);
+      hdr->Hubble *= 1000.*(one_kpc/one_Gyr);
+      h0 = hdr->Hubble / 100.;
+    }
+  else
+    {
+      float Or, Of;
+      SDFgetfloatOrDie(sdfp, "Omega0_m", &hdr->Omega_M);
+      SDFgetfloatOrDie(sdfp, "Omega0_lambda", &hdr->Omega_Lambda);
+      SDFgetfloatOrDie(sdfp, "Omega0_r", &Or);
+      SDFgetfloatOrDie(sdfp, "Omega0_fld", &Of);
+
+      hdr->Omega_M += Or;
+      hdr->Omega_Lambda += Of;
+      SDFgetfloatOrDie(sdfp, "h_100", &hdr->Hubble);
+      hdr->Hubble *= 100.;
+      h0 = hdr->Hubble / 100.; 
+    }
+
+  if (SDFhasname("R0", sdfp))
+    {
+      double R0;
+
+      SDFgetdoubleOrDie(sdfp, "R0", &R0);
+      hdr->BoxSize = 2.0*0.001*R0*h0;
+    }
+
+  if (SDFhasname("Rx", sdfp))
+    {
+      double R0;
+
+      SDFgetdoubleOrDie(sdfp, "Rx", &R0);
+      hdr->BoxSize = 2.0*0.001*R0*h0;
+    }
+
+  if (SDFhasname("redshift", sdfp))
+    {
+      double redshift;
+
+      SDFgetdoubleOrDie(sdfp, "redshift", &redshift);
+      hdr->time = 1/(1+redshift);
+    }
+   
+  if (SDFgetint64(sdfp, "npart", &gnobj))
+    {
+      gnobj = SDFnrecs("x", sdfp);
+      cerr << format("[Warning] No 'npart' found in SDF file '%s', guessing npart=%ld from SDFnrecs") % snapshot % gnobj << endl;
+    }
+  hdr->NumPart = hdr->TotalNumPart = gnobj;
+
+  return new SDFLoader(sdfp, hdr, flags, num_splitting, p);
+}
+
+
+void sdfLoaderInit(int& argc, char **& argv)
+{
+}

c_source/prep/loaders/sdfloader_internal.hpp

@@ -0,0 +1,27 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/sdfloader_internal.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef __VOID_SDF_LOADER_INTERNAL_HPP
+#define __VOID_SDF_LOADER_INTERNAL_HPP
+
+void sdfLoaderInit(int& argc, char **&argv);
+
+#endif

c_source/prep/loaders/simulation_loader.cpp

@@ -0,0 +1,110 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/simulation_loader.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cmath>
+#include <CosmoTool/loadSimu.hpp>
+#include "simulation_loader.hpp"
+#ifdef SDF_SUPPORT
+#include "sdfloader_internal.hpp"
+#endif
+
+using std::min;
+using namespace CosmoTool;
+
+template<typename T> void reallocArray(T *& a, long newSize, long toCopy)
+{
+  T *b = new T[newSize];
+  if (a != 0)
+    {
+      memcpy(b, a, sizeof(T)*toCopy);
+      delete[] a;
+    }
+  a = b;
+}
+
+void SimulationLoader::reallocSimu(SimuData *s, long newNumPart)
+{
+  long to_copy = min(newNumPart, s->NumPart);
+  
+  for (int j = 0; j < 3; j++)
+    {
+      reallocArray(s->Pos[j], newNumPart, to_copy);
+      reallocArray(s->Vel[j], newNumPart, to_copy);
+    }
+  reallocArray(s->Id, newNumPart, to_copy);
+  reallocArray(s->type, newNumPart, to_copy);
+}
+
+void SimulationLoader::applyTransformations(SimuData *s)
+{
+  float redshift_gravity = do_redshift ? 1.0 : 0.0;
+
+  for (int i = 0; i < s->NumPart; i++)
+    {
+      s->Pos[redshift_axis][i] += 
+	redshift_gravity*s->Vel[redshift_axis][i]/100.;
+    }
+}
+
+
+void SimulationLoader::basicPreprocessing(SimuData *d, 
+                                          SimulationPreprocessor *preproc)
+{
+  if (preproc == 0)
+    return;
+
+  long numAccepted = 0;
+  bool *accepted = new bool[d->NumPart];
+  for (long i = 0; i < d->NumPart; i++)
+    {
+      SingleParticle p;
+      
+      for (int k = 0; k < 3; k++)
+        {
+          p.Pos[k] = (d->Pos[k]) ? 0 : d->Pos[k][i];
+          p.Vel[k] = (d->Vel[k]) ? 0 : d->Vel[k][i];
+        }
+      p.ID = (d->Id) ? 0 : d->Id[i];
+      
+      accepted[i] = preproc->accept(p);
+      numAccepted += accepted[i];
+    }
+  
+  for (int k = 0; k < 3; k++)
+    {
+      filteredCopy(d->Pos[k], accepted, d->NumPart);
+      filteredCopy(d->Vel[k], accepted, d->NumPart);
+    }
+  filteredCopy(d->Id, accepted, d->NumPart);
+  filteredCopy(d->type, accepted, d->NumPart);
+  
+  filterAttribute<long>(d, "uniqueID", accepted, d->NumPart);
+  
+  d->NumPart = numAccepted;
+  delete[] accepted;
+}
+
+void simulationLoadersInit(int& argc, char **& argv)
+{
+#ifdef SDF_SUPPORT
+  sdfLoaderInit(argc, argv);
+#endif
+}

c_source/prep/loaders/simulation_loader.hpp

@@ -0,0 +1,156 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/loaders/simulation_loader.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef _MOCK_SIMULATION_LOADER_HPP
+#define _MOCK_SIMULATION_LOADER_HPP
+
+#include <string>
+#include <algorithm>
+#include <CosmoTool/loadSimu.hpp>
+
+struct SingleParticle
+{
+  float Pos[3];
+  float Vel[3];
+  long Index;
+  long ID;
+};
+
+class SimulationPreprocessor
+{
+public:
+  SimulationPreprocessor() {}
+  virtual ~SimulationPreprocessor() {}
+
+  virtual long getEstimatedPostprocessed(long numParticles) { return numParticles; };
+  virtual bool accept(const SingleParticle& p) { return true; }
+  virtual void reset() {}
+};
+
+class SimulationLoader
+{
+protected:
+  bool do_redshift;
+  int redshift_axis;
+
+  SimulationLoader() 
+  {
+    do_redshift = false;
+    redshift_axis = 2;
+  }
+
+  template<typename T> void reallocArray(T *& a, long newSize, long toCopy)
+  {
+    T *b = new T[newSize];
+    if (a != 0)
+     {
+       std::copy(a, a+toCopy, b);
+       delete[] a;
+     }
+   a = b;
+  } 
+
+
+
+  void reallocSimu(CosmoTool::SimuData *s, long newNumPart);
+ 
+
+  void basicPreprocessing(CosmoTool::SimuData *d, SimulationPreprocessor *preproc);
+  void applyTransformations(CosmoTool::SimuData *s);
+
+  void copyParticleToSimu(const SingleParticle& p, CosmoTool::SimuData *s, long index)
+  {
+    s->Pos[0][index] = p.Pos[0];
+    s->Pos[1][index] = p.Pos[1];
+    s->Pos[2][index] = p.Pos[2];
+    if (s->Vel[0])
+      s->Vel[0][index] = p.Vel[0];
+    if (s->Vel[1])
+      s->Vel[1][index] = p.Vel[1];
+    if (s->Vel[2])
+      s->Vel[2][index] = p.Vel[2];
+    s->Id[index] = p.ID;
+  }
+public:
+  virtual ~SimulationLoader() {}
+  
+  void doRedshift(bool set = true) { do_redshift = set; }
+  void setVelAxis(int axis) { redshift_axis = axis; }
+
+  virtual CosmoTool::SimuData *getHeader() = 0;
+  virtual int num_files() = 0;
+  virtual CosmoTool::SimuData* loadFile(int id) = 0;
+
+  
+  template<typename T>
+  void filteredCopy(T *a, bool *accepted, long N)
+  {
+    long i = 0, j = 0;
+
+    if (a == 0)
+      return;
+
+    while (i < N)
+      {
+        if (!accepted[i])
+          {
+            i++;
+            continue;
+          }
+
+        a[j] = a[i];
+        j++;
+        i++;
+      }
+  }
+  
+  template<typename T>
+  void filterAttribute(CosmoTool::SimuData *d, const std::string& attrname, bool *accepted, long NumPart)
+  {
+    if (d->attributes.find(attrname) == d->attributes.end())
+      return;
+ 
+    long *l = d->as<T>(attrname);
+    filteredCopy<T>(l, accepted, NumPart);
+  }
+};
+
+
+template<typename T>
+void delete_adaptor(void *ptr)
+{
+  T *ptr_T = reinterpret_cast<T *>(ptr);
+
+  delete[] ptr_T;
+}
+
+
+// Unit length is the size of one Mpc in the simulation units
+SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags, int gadgetFormat, SimulationPreprocessor *p);
+SimulationLoader *flashLoader(const std::string& snapshot, int flags, SimulationPreprocessor *p);
+SimulationLoader *multidarkLoader(const std::string& snapshot, SimulationPreprocessor *p);
+SimulationLoader *ramsesLoader(const std::string& snapshot, int baseid, bool double_precision, int flags, SimulationPreprocessor *p);
+SimulationLoader *sdfLoader(const std::string& snapshot, int flags, int num_splitting, SimulationPreprocessor *p);
+
+/* This is required for some parallel I/O handler (thus MPI beneath it) */
+void simulationLoadersInit(int& argc, char **& argv);
+
+#endif

c_source/prep/prepObservation.cpp

@@ -0,0 +1,651 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/prepObservation.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <healpix_map.h>
+#include <healpix_map_fitsio.h>
+#include <boost/format.hpp>
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <string>
+#include "prepObservation_conf.h"
+#include "contour_pixels.hpp"
+#include <netcdf>
+#include <CosmoTool/fortran.hpp>
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_integration.h>
+
+#define LIGHT_SPEED 299792.458
+
+using namespace std;
+using boost::format;
+using namespace CosmoTool;
+using namespace netCDF;
+
+struct NYU_Data
+{
+  int index;
+  int sector;
+  int region;
+  double ra, dec;
+  double cz;
+  double fgotten;
+  double phi_z;
+  long uniqueID;
+};
+
+struct Position
+{
+  double xyz[3];
+};
+
+struct ParticleData
+{
+  vector<int> id_gal;
+  vector<double> ra;
+  vector<double> dec;
+  vector<double> redshift;
+  vector<double> catalogID;
+  vector<long> uniqueID;
+  int id_mask;
+  // PMS
+  int mask_index;
+  // END PMS
+  vector<Position> pos;
+  double box[3][2];
+  double Lmax;
+};
+
+typedef vector<NYU_Data> NYU_VData;
+
+// this defines the expansion function that we will integrate
+// Laveaux & Wandelt (2012) Eq. 24
+struct my_expan_params { double Om; double w0; double wa; };
+double expanFun (double z, void * p) {
+  struct my_expan_params * params = (struct my_expan_params *)p;
+  double Om = (params->Om);
+  double w0 = (params->w0);
+  double wa = (params->wa);
+
+  //const double h0 = 1.0;
+  const double h0 = 0.71;
+  double ez;
+
+  double wz = w0 + wa*z/(1+z);
+
+  ez = Om*pow(1+z,3) + (1.-Om);
+  //ez = Om*pow(1+z,3) + pow(h0,2) * (1.-Om)*pow(1+z,3+3*wz);
+
+  ez = sqrt(ez);
+  //ez = sqrt(ez)/h0;
+
+  ez = 1./ez;
+
+  return  ez;
+}
+
+void loadData(const string& fname, NYU_VData & data)
+{
+  ifstream f(fname.c_str());
+   
+  while (!f.eof())
+    {
+      NYU_Data d;
+      f >> d.index >> d.sector >> d.region >> d.ra >> d.dec >> d.cz >> d.fgotten >> d.phi_z;
+      // Maubert - avoid double counting of last particle in array if EOF is after a blank line
+      if (!f) 
+      { 
+      continue; 
+      } 
+      // End Maubert
+      d.uniqueID = d.index;
+      data.push_back(d);
+    }
+}
+
+void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data, 
+                            bool useComoving, double omegaM)
+{
+  double d2r = M_PI/180;
+
+  gsl_function expanF;
+  expanF.function = &expanFun;
+  struct my_expan_params expanParams;
+  double maxZ = 2.0, z, result, error, *dL, *redshifts;
+  int numZ = 1000, iZ;
+  size_t nEval;
+
+  expanParams.Om = omegaM;
+  expanParams.w0 = -1.0;
+  expanParams.wa = 0.0;
+  expanF.params = &expanParams;
+
+  dL = (double *) malloc(numZ * sizeof(double));
+  redshifts = (double *) malloc(numZ * sizeof(double));
+
+  for (iZ = 0; iZ < numZ; iZ++) {
+    z = iZ * maxZ/numZ;
+    //gsl_integration_qng(&expanF, 0.0, z, 1.e-6, 1.e-6, &result, &error, &nEval);
+    dL[iZ] = result;
+    redshifts[iZ] = z;
+  }
+
+  gsl_interp *interp = gsl_interp_alloc(gsl_interp_linear, numZ); 
+  gsl_interp_init(interp, redshifts, dL, numZ);
+  gsl_interp_accel *acc = gsl_interp_accel_alloc();
+
+  output_data.pos.resize(data.size());
+  output_data.id_gal.resize(data.size());
+  output_data.ra.resize(data.size());
+  output_data.dec.resize(data.size());
+  output_data.redshift.resize(data.size());
+  output_data.uniqueID.resize(data.size());
+
+  for (int j = 0; j < 3; j++)
+    {
+      output_data.box[j][0] = -INFINITY;
+      output_data.box[j][1] = INFINITY;
+    }
+
+  for (int i = 0; i < data.size(); i++)
+    {  
+      double ra = data[i].ra*d2r, dec = data[i].dec*d2r;
+      Position& p = output_data.pos[i];
+
+      if (useComoving) {
+        //double pos = gsl_interp_eval(interp, redshifts, dL, data[i].cz, acc);
+        // Maubert - Lower boundary in redshift set to 0 to be consistent with pruneVoids (was 1.e-6 before).
+        gsl_integration_qng(&expanF, 0.0, data[i].cz/LIGHT_SPEED,
+                                        1.e-6, 
+                                        1.e-6, &result, &error, &nEval);
+        double Dc = result*LIGHT_SPEED;
+        p.xyz[0] = Dc*cos(ra)*cos(dec);
+        p.xyz[1] = Dc*sin(ra)*cos(dec);
+        p.xyz[2] = Dc*sin(dec);
+      } else {
+        p.xyz[0] = data[i].cz*cos(ra)*cos(dec);
+        p.xyz[1] = data[i].cz*sin(ra)*cos(dec);
+        p.xyz[2] = data[i].cz*sin(dec);
+      }
+//printf("CREATE %e %e\n", data[i].cz, sqrt(p.xyz[0]*p.xyz[0] + p.xyz[1]*p.xyz[1] + p.xyz[2]*p.xyz[2]));
+      output_data.id_gal[i] = data[i].index;
+      output_data.ra[i] = ra;
+      output_data.dec[i] = dec;
+      output_data.redshift[i] = data[i].cz;
+      output_data.uniqueID[i] = data[i].uniqueID;
+
+      for (int j = 0; j < 3; j++)
+	{
+	  if (p.xyz[j] > output_data.box[j][0])
+	    output_data.box[j][0] = p.xyz[j];
+	  if (p.xyz[j] < output_data.box[j][1])
+	    output_data.box[j][1] = p.xyz[j];	  
+	}
+//printf("INSERT GAL %d %e %e %e\n", output_data.id_gal[i], p.xyz[0], p.xyz[1], p.xyz[2]);
+    }
+
+  // normalize box
+  float left = 1.e99;
+  float right = -1.e99;
+  for (int j = 0; j < 3; j++) {
+    if (output_data.box[j][1] < left) left = output_data.box[j][1];
+    if (output_data.box[j][0] > right) right = output_data.box[j][0];
+  }
+  for (int j = 0; j < 3; j++) {
+    output_data.box[j][1] = left;
+    output_data.box[j][0] = right;
+  }
+
+  cout << format("Galaxy position generated: %d galaxies") % output_data.pos.size() << endl;
+  cout << format("box is %g < x < %g;  %g < y < %g;   %g < z < %g") 
+    % (1e-2*output_data.box[0][1]) % (1e-2*output_data.box[0][0])
+    % (1e-2*output_data.box[1][1]) % (1e-2*output_data.box[1][0])
+    % (1e-2*output_data.box[2][1]) % (1e-2*output_data.box[2][0]) << endl;
+ 
+  gsl_interp_free(interp); 
+}
+
+void generateSurfaceMask(prepObservation_info& args ,
+			 Healpix_Map<float>& mask, 
+			 vector<int>& pixel_list, 
+			 vector<int>& full_mask_list, 
+			 NYU_VData& data, 
+			 ParticleData& output_data,
+			 bool useComoving, 
+			 double omegaM)
+{
+
+  //Maubert - Needed for comobile distance in mock_sphere
+  gsl_function expanF;
+  expanF.function = &expanFun;
+  struct my_expan_params expanParams;
+  expanParams.Om = omegaM;
+  expanParams.w0 = -1.0;
+  expanParams.wa = 0.0;
+  expanF.params = &expanParams;
+  
+  double result, error ;
+  size_t nEval;
+  //End Maubert - Needed for comobile distance in mock_sphere
+
+  // Find the first free index
+  int idx = -1;
+  int insertion = 0;
+  double volume = pixel_list.size()*1.0/mask.Npix()*4*M_PI;
+  int numToInsert;
+
+  for (int i = 0; i < output_data.id_gal.size(); i++)
+    {
+      if (idx < output_data.id_gal[i])
+	idx = output_data.id_gal[i]+1;
+    }
+
+  output_data.id_mask = idx;
+ 
+// PMS
+  output_data.mask_index = output_data.id_gal.size();
+// END PMS 
+  cout << "Generate surface mask..." << endl;
+
+  double Rmax = -1;
+  for (int j = 0; j < 3; j++)
+    {
+      Rmax = max(Rmax, max(output_data.box[j][0], -output_data.box[j][1]));
+    }
+
+  output_data.Lmax = Rmax;
+
+// PMS - write a small text file with galaxy position
+  FILE *fp;
+  fp = fopen("galaxies.txt", "w");
+  for (int i = 0; i < data.size(); i++) {  
+    Position& p = output_data.pos[i];
+    fprintf(fp, "%e %e %e\n", 
+            (p.xyz[0]), 
+            (p.xyz[1]), 
+            (p.xyz[2]));
+  }
+  fclose(fp);
+// END PMS
+
+  cout << format("Rmax is %g, surface volume is %g") % (Rmax/100) % (volume/(4*M_PI)) << endl;
+  volume *= Rmax*Rmax*Rmax/3/1e6;
+  numToInsert = (int)floor(volume*args.density_fake_arg);
+  cout << format("3d volume to fill: %g (Mpc/h)^3") % volume << endl;
+
+  cout << format("Will insert %d particles") % numToInsert << endl;
+
+  fp = fopen("mock_galaxies.txt", "w");
+
+  double pct = 0;
+  for (int i = 0; i < numToInsert; i++) {
+    double new_pct = i*100./numToInsert;      
+
+    if (new_pct-pct > 5.) {
+	    pct = new_pct;
+	    cout << format(" .. %3.0f %%") % pct << endl;
+ 	  }
+
+    Position p;
+    bool stop_here;
+
+    do {
+	    int p0 = (int)floor(drand48()*pixel_list.size());
+	    vec3 v = mask.pix2vec(pixel_list[p0]);
+	    double r = Rmax*pow(drand48(),1./3);
+	  
+	    p.xyz[0] = v.x * r;
+	    p.xyz[1] = v.y * r;
+	    p.xyz[2] = v.z * r;
+
+	    stop_here = true;
+	    for (int j = 0; j < 3; j++) {
+	      if (p.xyz[j] > output_data.box[j][0] ||
+		        p.xyz[j] < output_data.box[j][1])
+		      stop_here = false;
+	    }
+	  }
+    while (!stop_here);
+      
+// PMS : write mock galaxies to a small file for plotting
+    fprintf(fp, "%e %e %e\n", 
+            (p.xyz[0]), 
+            (p.xyz[1]), 
+            (p.xyz[2]));
+// END PMS
+    output_data.pos.push_back(p);
+    output_data.id_gal.push_back(idx);
+    output_data.ra.push_back(-1);
+    output_data.dec.push_back(-1);
+    output_data.redshift.push_back(-1);
+    output_data.uniqueID.push_back(-1);
+//printf("INSERT MOCK %d %e %e %e\n", idx, p.xyz[0], p.xyz[1], p.xyz[2]);
+    insertion++;
+  }
+
+  fclose(fp);
+
+  // PMS
+  // TEST - insert mock galaxies along box edge
+  fp = fopen("mock_boundary.txt", "w");
+  double dx[3];
+  dx[0] = output_data.box[0][1] - output_data.box[0][0];
+  dx[1] = output_data.box[1][1] - output_data.box[1][0];
+  dx[2] = output_data.box[2][1] - output_data.box[2][0];
+
+  int nPart = 100;
+// TEST
+  for (int iDir = 0; iDir < 0; iDir++) {
+  for (int iFace = 0; iFace < 0; iFace++) {
+  //for (int iDir = 0; iDir < 3; iDir++) {
+  //for (int iFace = 0; iFace < 2; iFace++) {
+
+    int iy = (iDir + 1) % 3;
+    int iz = (iDir + 2) % 3;
+
+    for (int i = 0; i < nPart; i++) {
+    for (int j = 0; j < nPart; j++) {
+      Position p;
+
+      p.xyz[iDir] = output_data.box[iDir][iFace];
+      p.xyz[iy] = i * dx[iy]/nPart + output_data.box[iy][0];
+      p.xyz[iz] = j * dx[iz]/nPart + output_data.box[iz][0];
+
+      output_data.pos.push_back(p);
+      output_data.id_gal.push_back(idx);
+      output_data.ra.push_back(-1);
+      output_data.dec.push_back(-1);
+      output_data.redshift.push_back(-1);
+      output_data.uniqueID.push_back(-1);
+      insertion++;
+
+      fprintf(fp, "%e %e %e\n", 
+              (p.xyz[0]),
+              (p.xyz[1]), 
+              (p.xyz[2]));
+    } 
+    }
+  } 
+  } 
+  fclose(fp);
+  // END PMS
+
+   // PMS
+  // TEST - insert mock galaxies along spheres of survey redshift boundaries
+  fp = fopen("mock_sphere.txt", "w");
+  //Maubert - insert mock galaxies according to useComoving specification
+  // Added & compute rmin & rmax out of loop
+  double rmin ;
+  double rmax ; 
+  if (useComoving) {
+        gsl_integration_qng(&expanF, 0.0, args.zMin_arg,
+                                        1.e-6, 
+                                        1.e-6, &result, &error, &nEval);
+        rmin = result* LIGHT_SPEED;
+        
+        gsl_integration_qng(&expanF, 0.0, args.zMax_arg,
+                                        1.e-6, 
+                                        1.e-6, &result, &error, &nEval);
+        rmax = result* LIGHT_SPEED;
+  } else {
+         rmin = args.zMin_arg * LIGHT_SPEED;
+         rmax = args.zMax_arg * LIGHT_SPEED;
+  }
+    
+  //for (int q = 0; q < 0; q++) {
+  for (int q = 0; q < full_mask_list.size(); q++) {
+    vec3 v = mask.pix2vec(full_mask_list[q]);
+	  
+    Position p;
+
+
+    if (rmin > 0.) {
+      p.xyz[0] = v.x * rmin;
+      p.xyz[1] = v.y * rmin;
+      p.xyz[2] = v.z * rmin;
+      output_data.pos.push_back(p);
+      output_data.id_gal.push_back(idx);
+      output_data.ra.push_back(-1);
+      output_data.dec.push_back(-1);
+      output_data.redshift.push_back(-1);
+      output_data.uniqueID.push_back(-1);
+      insertion++;
+      fprintf(fp, "%e %e %e\n", 
+              (p.xyz[0]),
+              (p.xyz[1]), 
+              (p.xyz[2]));
+    }
+
+
+    p.xyz[0] = v.x * rmax;
+    p.xyz[1] = v.y * rmax;
+    p.xyz[2] = v.z * rmax;
+    output_data.pos.push_back(p);
+    output_data.id_gal.push_back(idx);
+    output_data.ra.push_back(-1);
+    output_data.dec.push_back(-1);
+    output_data.redshift.push_back(-1);
+    output_data.uniqueID.push_back(-1);
+    insertion++;
+    fprintf(fp, "%e %e %e\n", 
+            (p.xyz[0]),
+            (p.xyz[1]), 
+            (p.xyz[2]));
+  }
+  fclose(fp);
+  // END PMS
+
+  cout << format("Done. Inserted %d particles.") % insertion << endl;
+}
+
+void saveData(ParticleData& pdata)
+{
+  NcFile f("particles.nc", NcFile::replace);
+
+  NcDim d = f.addDim("space", 3);
+  NcDim p = f.addDim("Np", pdata.pos.size());
+  NcVar v = f.addVar("particles", ncDouble, {d, p});
+  double *x = new double[pdata.pos.size()];
+
+  for (int j = 0; j < 3; j++)
+    {
+
+      for (int i = 0; i < pdata.pos.size(); i++)
+	x[i] = pdata.pos[i].xyz[j];
+
+      v.putVar({size_t(j), 0}, {1, pdata.pos.size()}, x);
+    }
+
+  v = f.addVar("id_gal", ncInt, std::vector<NcDim>({p}));
+  v.putVar(&pdata.id_gal[0]);
+
+  delete[] x;
+  
+}
+
+void saveForZobov(ParticleData& pdata, const string& fname, const string& paramname)
+{
+  UnformattedWrite f(fname);
+  static const char axis[] = { 'X', 'Y', 'Z' };
+  double Lmax = pdata.Lmax;
+  double r2d = 180./M_PI;
+
+  f.beginCheckpoint();
+  f.writeInt32(pdata.pos.size());
+  f.endCheckpoint();
+
+
+  for (int j = 0; j < 3; j++)
+    {
+      cout << format("Writing %c components...") % axis[j] << endl;
+      f.beginCheckpoint();
+      for (uint32_t i = 0; i < pdata.pos.size(); i++)
+	{
+	  f.writeReal32((pdata.pos[i].xyz[j]+Lmax)/(2*Lmax));
+	}
+      f.endCheckpoint();
+    }
+
+  cout << format("Writing RA...")  << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < pdata.pos.size(); i++) {
+	  f.writeReal32(pdata.ra[i]*r2d);
+	}
+  f.endCheckpoint();
+   
+  cout << format("Writing Dec...")  << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < pdata.pos.size(); i++) {
+	  f.writeReal32(pdata.dec[i]*r2d);
+	}
+  f.endCheckpoint();
+    
+  cout << format("Writing Redshift...")  << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < pdata.pos.size(); i++) {
+	  f.writeReal32(pdata.redshift[i]);
+	}
+  f.endCheckpoint();
+   
+  cout << format("Writing Unique ID...")  << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < pdata.pos.size(); i++) {
+	  f.writeInt64(pdata.uniqueID[i]);
+	}
+  f.endCheckpoint();
+   
+  NcFile fp(paramname.c_str(), NcFile::replace);
+
+  fp.putAtt("range_x_min", ncDouble, -Lmax/100.);
+  fp.putAtt("range_x_max", ncDouble, Lmax/100.);
+  fp.putAtt("range_y_min", ncDouble, -Lmax/100.);
+  fp.putAtt("range_y_max", ncDouble, Lmax/100.);
+  fp.putAtt("range_z_min", ncDouble, -Lmax/100.);
+  fp.putAtt("range_z_max", ncDouble, Lmax/100.);
+  fp.putAtt("mask_index", ncInt, pdata.mask_index); // PMS
+  fp.putAtt("is_observation", ncInt, 1); // PMS
+
+  int nOutputPart = pdata.mask_index;
+  //int nOutputPart = pdata.pos.size();
+
+  NcDim NumPart_dim = fp.addDim("numpart_dim", nOutputPart);
+  NcVar v = fp.addVar("particle_ids", ncInt, NumPart_dim);
+  //NcVar v2 = fp.addVar("expansion", ncDouble, NumPart_dim);
+
+  //double *expansion_fac = new double[pdata.pos.size()];
+
+  //for (int i = 0; i <  pdata.pos.size(); i++)
+  //  expansion_fac[i] = 1.0;
+
+  v.putVar({0}, {size_t(nOutputPart)}, &pdata.id_gal[0]);
+  //v2->put(expansion_fac, pdata.pos.size());
+
+  //delete[] expansion_fac;
+
+/*
+  FILE *infoFile = fopen("sample_info.txt", "w");
+  fprintf(infoFile, "x_min = %f\n", -Lmax/100.);  
+  fprintf(infoFile, "x_max = %f\n", Lmax/100.);
+  fprintf(infoFile, "y_min = %f\n", -Lmax/100.);  
+  fprintf(infoFile, "y_max = %f\n", Lmax/100.);  
+  fprintf(infoFile, "z_min = %f\n", -Lmax/100.);
+  fprintf(infoFile, "z_max = %f\n", Lmax/100.);
+  fprintf(infoFile, "mask_index = %d\n", pdata.mask_index);
+  fprintf(infoFile, "total_particles = %d\n", pdata.pos.size());
+  fclose(infoFile);
+*/
+
+}
+
+int main(int argc, char **argv)
+{
+  prepObservation_info args_info;
+  prepObservation_conf_params args_params;
+ 
+  prepObservation_conf_init(&args_info);
+  prepObservation_conf_params_init(&args_params);
+  
+  args_params.check_required = 0;
+  if (prepObservation_conf_ext (argc, argv, &args_info, &args_params))
+    return 1;
+  
+  if (!args_info.configFile_given)
+    {
+      if (prepObservation_conf_required (&args_info, PREPOBSERVATION_CONF_PACKAGE))
+        return 1;
+    }
+  else
+    {
+      args_params.check_required = 1;
+      args_params.initialize = 0;
+      if (prepObservation_conf_config_file (args_info.configFile_arg,
+					 &args_info,
+					 &args_params))
+	return 1;
+    }
+  
+  prepObservation_conf_print_version();
+
+  cout << "Loading data " << args_info.catalog_arg << "..." << endl;
+  vector<NYU_Data> data;
+  Healpix_Map<float> o_mask;
+  vector<int> pixel_list;
+  vector<int> full_mask_list;
+  ParticleData output_data;
+
+  loadData(args_info.catalog_arg, data);
+  
+
+  cout << "Loading mask..." << endl;
+  read_Healpix_map_from_fits(args_info.mask_arg, o_mask);
+
+  Healpix_Map<float> mask;
+
+  mask.SetNside(128, RING);
+  mask.Import(o_mask);
+
+  computeContourPixels(mask,pixel_list);
+  computeMaskPixels(mask,full_mask_list);
+
+  // We compute a cube holding all the galaxies + the survey surface mask
+
+  cout << "Placing galaxies..." << endl;
+  placeGalaxiesInCube(data, output_data, args_info.useComoving_flag, 
+                         args_info.omegaM_arg);
+  generateSurfaceMask(args_info, mask, pixel_list, full_mask_list,
+                      data, output_data,args_info.useComoving_flag, 
+                         args_info.omegaM_arg);
+  
+  saveForZobov(output_data, args_info.output_arg, args_info.params_arg);
+  //  saveData(output_data);
+  
+  // PMS
+  FILE *fp = fopen("mask_index.txt", "w");
+  fprintf(fp, "%d", output_data.mask_index);
+  fclose(fp);
+
+  fp = fopen("total_particles.txt", "w");
+  fprintf(fp, "%d", output_data.pos.size());
+  fclose(fp);
+  printf("Done!\n");
+  // END PMS
+  return 0;
+}

c_source/prep/prepObservation.ggo

@@ -0,0 +1,19 @@
+package "prepObservation"
+version "alpha"
+
+option "configFile" - "Configuration filename" string optional
+
+option "catalog" - "Input NYU-VAGC catalog" string required
+option "mask" - "Healpix mask of unobserved data (in Equatorial coordinates)" string required
+option "density_fake" - "Number density of boundary fake tracers (1 h^3/ Mpc^3)"  double optional default="1"
+
+option "zMin" - "Minimum redshift of data" double required
+
+option "zMax" - "Maximum redshift of data" double required
+
+option "output"	- "Filename of particle datafile" string required
+option "params" - "Output parameters of the datacube" string required
+
+option "useComoving" - "Convert to real space using LCDM cosmology" flag off
+
+option "omegaM" - "Omega Matter for fiducial cosmology" double optional default="0.27"

c_source/prep/prepSimulation.cpp

@@ -0,0 +1,791 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/mock/prepSimulation.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <gsl/gsl_rng.h>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/format.hpp>
+#include <cmath>
+#include <cassert>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <algorithm>
+#include <CosmoTool/loadSimu.hpp>
+#include <CosmoTool/interpolate.hpp>
+#include <CosmoTool/fortran.hpp>
+#include <CosmoTool/algo.hpp>
+#include "prepSimulation_conf.h"
+#include "gslIntegrate.hpp"
+#include <netcdf>
+#include "simulation_loader.hpp"
+
+using namespace std;
+using namespace CosmoTool;
+using boost::format;
+using namespace netCDF;
+
+#define LIGHT_SPEED 299792.458
+
+typedef boost::function2<void, SimuData*, double*> MetricFunctor;
+
+struct TotalExpansion
+{
+  double Omega_M, Omega_L;
+
+  double operator()(double z)
+  {
+    return 1/sqrt(Omega_M*cube(1+z) + Omega_L);
+  }
+};
+
+Interpolate make_cosmological_redshift(double OM, double OL, double z0, double z1, int N = 5000)
+{
+  TotalExpansion e_computer;
+  double D_tilde, Q, Qprime;
+  InterpolatePairs pairs;
+  
+  e_computer.Omega_M = OM;
+  e_computer.Omega_L = OL;
+  
+  pairs.resize(N);
+  ofstream f("comoving_distance.txt");
+
+  for (int i = 0; i < N; i++)
+    {      
+      double z = z0 + (z1-z0)/N*i;
+
+      pairs[i].second = z;
+      pairs[i].first = gslIntegrate(e_computer, 0, z, 1e-3);
+      f << z << " " <<  pairs[i].first << endl;
+    }
+
+  return buildFromVector(pairs);
+}
+
+void metricTransform(SimuData *data, int axis, bool reshift, bool pecvel, double* expfact, bool cosmo_flag)
+{
+  int x0, x1, x2;
+
+  switch (axis) {
+  case 0:
+    x0 = 1; x1 = 2; x2 = 0;
+    break;
+  case 1:
+    x0 = 0; x1 = 2; x2 = 1;
+    break;
+  case 2:
+    x0 = 0; x1 = 1; x2 = 2;
+    break;
+  default:
+    abort();
+  }
+
+  double z0 = 1/data->time - 1;
+  Interpolate z_vs_D = 
+    make_cosmological_redshift(data->Omega_M, data->Omega_Lambda, 
+			       0., z0+8*data->BoxSize*100/LIGHT_SPEED);
+  // Redshift 2*z0 should be sufficient ? This is fragile. 
+  //A proper solver is needed here.
+  double z_base = reshift ? z0 : 0;
+  TotalExpansion e_computer;
+  double baseComovingDistance;
+
+  cout << "Using base redshift z=" << z0 << " " << z0+8*data->BoxSize*100/LIGHT_SPEED << endl;
+
+  e_computer.Omega_M = data->Omega_M;
+  e_computer.Omega_L = data->Omega_Lambda;
+  baseComovingDistance = LIGHT_SPEED/100.* gslIntegrate(e_computer, 0, z0, 1e-3);
+  cout << "Comoving distance = " << baseComovingDistance << " Mpc/h" << endl;
+
+  if (cosmo_flag) cout << "Will place particles on a lightcone..." << endl;
+
+  float minZ = 1.e99;
+  float maxZ = 0;
+
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      float& x = data->Pos[x0][i];
+      float& y = data->Pos[x1][i];
+      float& z = data->Pos[x2][i];
+      float& v = data->Vel[2][i];
+      float z_old = z;
+
+      double reduced_red = (z + baseComovingDistance)*100./LIGHT_SPEED;
+      double reduced_base = (reshift ? (baseComovingDistance*100./LIGHT_SPEED) : 0);
+      try
+        {
+
+          // Distorted redshift
+          if (reduced_red == 0)
+            z = 0;
+          else if (cosmo_flag)
+	    z = (z_vs_D.compute(reduced_red)-z_base)*LIGHT_SPEED/100.;
+          else
+            z = (reduced_red-reduced_base)*LIGHT_SPEED/100.0;
+
+          if (expfact)
+            expfact[i] = z / z_old; 
+	  // Add peculiar velocity
+	  if (pecvel)
+	    z += v/(100*e_computer(z0));
+       }
+      catch(const InvalidRangeException& e) {
+       cout << "Trying to interpolate out of the tabulated range." << endl;
+       cout << "The offending value is z=" << reduced_red << endl;
+       abort();
+      }
+      if (z > maxZ) maxZ = z;
+      if (z < minZ) minZ = z;
+    }
+
+    printf("Range of z: %.2f - %.2f\n", minZ, maxZ);
+}
+
+// slightly perturb particle positions
+void joggleParticles(SimuData *data) {
+  cout << "Joggling particle positions..." << endl;
+  gsl_rng *myRng = gsl_rng_alloc(gsl_rng_taus);
+  int seed = 314159;
+  gsl_rng_set(myRng, seed);
+  for (uint32_t i = 0; i < data->NumPart; i++) {
+    data->Pos[0][i] += 1.e-3*gsl_rng_uniform(myRng);
+    data->Pos[1][i] += 1.e-3*gsl_rng_uniform(myRng); 
+    data->Pos[2][i] += 1.e-3*gsl_rng_uniform(myRng);
+    data->Pos[0][i] -= 1.e-3*gsl_rng_uniform(myRng);
+    data->Pos[1][i] -= 1.e-3*gsl_rng_uniform(myRng); 
+    data->Pos[2][i] -= 1.e-3*gsl_rng_uniform(myRng);
+  }
+} // end joggleParticles
+
+void generateOutput(SimuData *data, int axis, 
+		    const std::string& fname)
+{
+  UnformattedWrite f(fname);
+
+  cout << "Generating output particles to " << fname << endl;
+ 
+  int x0, x1, x2;
+  
+  switch (axis) {
+  case 0:
+    x0 = 1; x1 = 2; x2 = 0;
+    break;
+  case 1:
+    x0 = 0; x1 = 2; x2 = 1;
+    break;
+  case 2:
+    x0 = 0; x1 = 1; x2 = 2;
+    break;
+  default:
+    abort();
+  }
+
+  f.beginCheckpoint();
+  f.writeInt32(data->NumPart);
+  f.endCheckpoint();
+
+  cout << "Writing X components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x0][i]);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing Y components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x1][i]);
+    }
+  f.endCheckpoint();
+  
+  cout << "Writing Z components..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x2][i]);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing RA..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x0][i]);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing Dec..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x1][i]);
+    }
+  f.endCheckpoint();
+
+  cout << "Writing redshift..." << endl;
+  f.beginCheckpoint();
+  for (uint32_t i = 0; i < data->NumPart; i++)
+    {
+      f.writeReal32(data->Pos[x2][i]*LIGHT_SPEED);
+    }
+  f.endCheckpoint();
+
+  long *uniqueID = data->as<long>("uniqueID");
+  if (uniqueID != 0)
+    {
+      cout << "Writing unique ID..." << endl;
+      f.beginCheckpoint();
+      for (uint32_t i = 0; i < data->NumPart; i++)
+	{
+	  f.writeInt64(uniqueID[i]);
+	}
+      f.endCheckpoint();
+    }
+}
+
+// This function prepares the list of targets for the specified snapshot. The target list is not
+// cleared. That way new particles can be appended if this is a multi-file snapshot.
+void selectBox(SimuData *simu, std::vector<long>& targets, prepSimulation_info& args_info, SimulationPreprocessor *preselect)
+{
+  double ranges[3][2] = {
+    { args_info.rangeX_min_arg, args_info.rangeX_max_arg },
+    { args_info.rangeY_min_arg, args_info.rangeY_max_arg },
+    { args_info.rangeZ_min_arg, args_info.rangeZ_max_arg }
+  };
+  long numAccepted =0;
+
+  for (uint32_t i = 0; i < simu->NumPart; i++)
+    {
+      bool acceptance = true;
+      SingleParticle p;
+
+      for (int j = 0; j < 3; j++) {
+	acceptance = 
+	  acceptance &&
+	  (simu->Pos[j][i] > ranges[j][0]) && 
+	  (simu->Pos[j][i] <= ranges[j][1]);	
+        p.Pos[j] = simu->Pos[j][i];
+        p.Vel[j] = (simu->Vel[j] != 0) ? simu->Vel[j][i] : 0;
+      }
+      p.ID = (simu->Id != 0) ? simu->Id[i] : -1;
+
+      if (preselect != 0)
+        acceptance = acceptance && preselect->accept(p);
+      
+      if (acceptance) {
+	targets.push_back(i);
+        numAccepted++; 
+      }
+    }
+  cout << "SELECTBOX: Num accepted here = " << numAccepted <<  " / input = " << simu->NumPart << " (after resubsampling)" << endl; 
+}
+
+class PreselectParticles: public SimulationPreprocessor
+{
+private:
+  gsl_rng *rng;
+  double subsample;
+  int seed;
+public:
+  PreselectParticles(double s, int seed_value)
+    : subsample(s), seed(seed_value), rng(gsl_rng_alloc(gsl_rng_default))
+  {
+    gsl_rng_set(rng, seed);
+  }
+
+  virtual ~PreselectParticles()
+  {
+    gsl_rng_free(rng);
+  } 
+
+  bool accept(const SingleParticle& p)
+  {
+    return gsl_rng_uniform(rng) < subsample;
+  }
+
+  void reset()
+  {
+    gsl_rng_set(rng, seed);
+  }
+};
+
+void createBox(SimuData *simu, vector<long>& targets, vector<long>& snapshot_split, SimuData *& boxed, prepSimulation_info& args_info)
+{
+  double *ranges = new double[6];
+  double *mul = new double[3];
+  long *simu_uniqueID = simu->as<long>("uniqueID");
+
+  ranges[0] = args_info.rangeX_min_arg;
+  ranges[1] = args_info.rangeX_max_arg;
+  ranges[2] = args_info.rangeY_min_arg;
+  ranges[3] = args_info.rangeY_max_arg;
+  ranges[4] = args_info.rangeZ_min_arg;
+  ranges[5] = args_info.rangeZ_max_arg;
+
+  boxed = new SimuData;
+  boxed->Hubble = simu->Hubble;
+  boxed->Omega_M = simu->Omega_M;
+  boxed->Omega_Lambda = simu->Omega_Lambda;
+  boxed->time = simu->time;
+  boxed->BoxSize = simu->BoxSize;
+  boxed->NumPart = targets.size();
+
+  for (int j = 0; j < 3; j++)
+    {
+      boxed->Pos[j] = new float[boxed->NumPart];
+      boxed->Vel[j] = new float[boxed->NumPart];
+      //boxed->Vel[j] = 0;
+      mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]);
+    }
+  cout << "Min range = " << ranges[0] << " " << ranges[2] << " " << ranges[4] << endl;
+  cout << "Max range = " << ranges[1] << " " << ranges[3] << " " << ranges[5] << endl;
+  cout << "Number of accepted particles: " << boxed->NumPart << endl; 
+  cout << "Rescaling factors = " << mul[0] << " " << mul[1] << " " << mul[2] << endl;
+
+  // PMS
+  FILE *fp = fopen("mask_index.txt", "w");
+  fprintf(fp, "%ld", boxed->NumPart);
+  fclose(fp);
+
+  fp = fopen("total_particles.txt", "w");
+  fprintf(fp, "%ld", boxed->NumPart);
+  fclose(fp);
+  printf("Done!\n");
+  // END PMS
+
+
+  long *uniqueID = new long[boxed->NumPart];
+  long *particle_id = new long[boxed->NumPart];
+  double *expansion_fac = new double[boxed->NumPart];
+  long *snap_split = new long[snapshot_split.size()];
+  int *numsnap_info = new int[1];
+
+  copy(targets.begin(), targets.end(), particle_id);
+  copy(snapshot_split.begin(), snapshot_split.end(), snap_split);
+  *numsnap_info = snapshot_split.size();
+
+  boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
+  boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
+  boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>); 
+  boxed->new_attribute("mul", mul, delete_adaptor<double>);
+  boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
+  boxed->new_attribute("snapshot_split", snap_split, delete_adaptor<long>);
+  boxed->new_attribute("num_snapshots", numsnap_info, delete_adaptor<int>);
+}
+
+void buildBox(SimuData *simu, long num_targets, long loaded, 
+	      SimuData *boxed, double *efac)
+{
+  uint32_t k = 0;
+  long *uniqueID = boxed->as<long>("uniqueID");
+  long *simu_uniqueID = simu->as<long>("uniqueID");
+  double *expansion_fac = boxed->as<double>("expansion_fac");
+  double *mul = boxed->as<double>("mul");
+  double *ranges = boxed->as<double>("ranges");
+  long *particle_id = boxed->as<long>("particle_id");
+
+  for (uint32_t i = 0; i < num_targets; i++, loaded++)
+    {
+      long pid = particle_id[loaded];
+      //assert(pid < simu->NumPart);
+      assert(loaded < boxed->NumPart);
+      
+      for (int j = 0; j < 3; j++)
+	{
+	  boxed->Pos[j][loaded] = max(min((simu->Pos[j][pid]-ranges[j*2])*mul[j], double(1)), double(0));
+	  boxed->Vel[j][loaded] = simu->Vel[j][pid];
+	  assert(boxed->Pos[j][loaded] >= 0);
+	  assert(boxed->Pos[j][loaded] <= 1);
+	}
+      uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[pid] : 0;
+      expansion_fac[loaded] = efac[pid];
+    }
+}
+
+void saveBox(SimuData *&boxed, const std::string& outbox, prepSimulation_info& args_info)
+{
+  double *ranges = boxed->as<double>("ranges");
+  NcFile f(outbox.c_str(), NcFile::replace);
+  long *particle_id = boxed->as<long>("particle_id");
+  double *expansion_fac = boxed->as<double>("expansion_fac");
+  long *snapshot_split = boxed->as<long>("snapshot_split");
+  int num_snapshots = *boxed->as<int>("num_snapshots");
+  long *uniqueID = boxed->as<long>("uniqueID");
+  float *velX = boxed->Vel[0];
+  float *velY = boxed->Vel[1];
+  float *velZ = boxed->Vel[2];
+
+  f.putAtt("range_x_min", ncDouble, ranges[0]);
+  f.putAtt("range_x_max", ncDouble, ranges[1]);
+  f.putAtt("range_y_min", ncDouble,  ranges[2]);
+  f.putAtt("range_y_max", ncDouble, ranges[3]);
+  f.putAtt("range_z_min", ncDouble, ranges[4]);
+  f.putAtt("range_z_max", ncDouble, ranges[5]);
+  f.putAtt("mask_index", ncInt, -1);
+  f.putAtt("is_observation", ncInt, 0);
+  f.putAtt("data_subsampling", ncInt, args_info.subsample_arg);
+
+  NcDim NumPart_dim = f.addDim("numpart_dim", boxed->NumPart);
+  NcDim NumSnap_dim = f.addDim("numsnap_dim", num_snapshots);
+  NcVar v = f.addVar("particle_ids", ncInt64, NumPart_dim);
+  NcVar v2 = f.addVar("expansion", ncDouble, NumPart_dim);
+  NcVar v3 = f.addVar("snapshot_split", ncInt64, NumSnap_dim);
+
+  v.putVar({0}, {size_t(boxed->NumPart)}, particle_id);
+  v2.putVar({0}, {size_t(boxed->NumPart)}, expansion_fac);  
+  v3.putVar({0}, {size_t(num_snapshots)}, snapshot_split);
+  if (uniqueID != 0)
+    {
+      NcVar v4 = f.addVar("unique_ids_lsb", ncInt, NumPart_dim);
+      NcVar v5 = f.addVar("unique_ids_msb", ncInt, NumPart_dim);
+
+      nclong *tmp_int = new nclong[boxed->NumPart];
+      for (long i = 0; i < boxed->NumPart; i++)
+         tmp_int[i] = (nclong)(((unsigned long)uniqueID[i]) & 0xffffffff);
+      v4.putVar({0}, {size_t(boxed->NumPart)}, tmp_int);
+      for (long i = 0; i < boxed->NumPart; i++)
+         tmp_int[i] = (nclong)((((unsigned long)uniqueID[i]) & 0xffffffff) >> 32);
+      v5.putVar({0}, {size_t(boxed->NumPart)}, tmp_int);
+      delete[] tmp_int;
+    }
+
+  NcVar v6 = f.addVar("vel_x", ncFloat, NumPart_dim);
+  NcVar v7 = f.addVar("vel_y", ncFloat, NumPart_dim);
+  NcVar v8 = f.addVar("vel_z", ncFloat, NumPart_dim);
+  v6.putVar({0}, {size_t(boxed->NumPart)}, velX);
+  v7.putVar({0}, {size_t(boxed->NumPart)}, velY);
+  v8.putVar({0}, {size_t(boxed->NumPart)}, velZ);
+}
+
+void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, prepSimulation_info& args_info)
+{
+  NcFile f(args_info.inputParameter_arg, NcFile::read);
+  NcVar *v;
+  long *particle_id;
+  double *expansion_fac;
+  long *uniqueID;
+  long *snapshot_split;
+  int *num_snapshots = new int[1];
+
+  boxed = new SimuData;
+  boxed->Hubble = simu->Hubble;
+  boxed->Omega_M = simu->Omega_M;
+  boxed->Omega_Lambda = simu->Omega_Lambda;
+  boxed->time = simu->time;
+  boxed->BoxSize = simu->BoxSize;
+
+  NcGroupAtt d_sub = f.getAtt("data_subsampling");
+  auto checkAtt = [&args_info](NcGroupAtt a) {
+    if (a.isNull())
+      return true;
+    double subsampling;
+    a.getValues(&subsampling);
+    return subsampling/args_info.subsample_arg - 1 > 1e-5;
+  };
+  if (checkAtt(d_sub))
+   {
+     cerr << "Parameter file was not generated with the same simulation subsampling argument. Particles will be different. Stop here." <<endl;
+     exit(1);
+   }
+
+  NcVar v_id = f.getVar("particle_ids");
+  NcVar v_snap = f.getVar("snapshot_split");
+  double *ranges;
+  double *mul;
+
+  std::vector<NcDim> edges1 = v_id.getDims();
+  std::vector<NcDim> dim_snap = v_snap.getDims();
+  assert(edges1.size()==1);
+  assert(dim_snap.size()==1);
+
+  boxed->NumPart = edges1[0].getSize();
+  *num_snapshots = dim_snap[0].getSize();
+
+  particle_id = new long[boxed->NumPart];
+  uniqueID = new long[boxed->NumPart];
+  mul = new double[3];
+  ranges = new double[6];
+  snapshot_split = new long[*num_snapshots];
+  expansion_fac = new double[boxed->NumPart];
+
+
+  boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
+  boxed->new_attribute("mul", mul, delete_adaptor<double>);
+  boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
+  boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
+  boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor<int>);
+  boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor<long>);
+  boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
+
+  v_id.getVar(particle_id);
+  v_snap.getVar(snapshot_split);
+
+  f.getAtt("range_x_min").getValues(&ranges[0]);
+  f.getAtt("range_x_max").getValues(&ranges[1]);
+  f.getAtt("range_y_min").getValues(&ranges[2]);
+  f.getAtt("range_y_max").getValues(&ranges[3]);
+  f.getAtt("range_z_min").getValues(&ranges[4]);
+  f.getAtt("range_z_max").getValues(&ranges[5]);
+
+  for (int j = 0; j < 3; j++)
+    {
+      boxed->Pos[j] = new float[boxed->NumPart];
+      boxed->Vel[j] = new float[boxed->NumPart];
+      //boxed->Vel[j] = 0;
+      mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]);
+    }
+  
+  uint32_t k = 0;
+  NcVar v_uniq_lsb = f.getVar("unique_ids_lsb");
+  NcVar v_uniq_msb = f.getVar("unique_ids_lsb");
+  nclong *tmp_int;
+
+  tmp_int = new nclong[boxed->NumPart];
+
+  v_uniq_lsb.getVar(tmp_int);
+  for (long i = 0; i < boxed->NumPart; i++)
+    uniqueID[i] = tmp_int[i];
+
+  v_uniq_msb.getVar(tmp_int);
+  for (long i = 0; i < boxed->NumPart; i++)
+    uniqueID[i] |= (unsigned long)(tmp_int[i]) << 32;
+
+  delete[] tmp_int;
+
+  PreselectParticles *preselect = 0;
+
+  if (args_info.resubsample_given)
+    {
+      preselect = new PreselectParticles(args_info.resubsample_arg, args_info.resubsample_seed_arg);
+      preselect->reset();
+    }
+
+  if (preselect == 0)
+    return;
+
+  long pid_read = 0, pid_write = 0;
+  for (int s_id = 0; s_id < *num_snapshots; s_id++)
+    {
+      long previous_write = pid_write;
+      for (long q = 0; q < snapshot_split[s_id]; q++)
+        {
+          SingleParticle p;
+          p.ID = -1;
+          assert(pid_read < boxed->NumPart);
+          if (preselect->accept(p))
+           {
+             particle_id[pid_write] = particle_id[pid_read];
+             uniqueID[pid_write] = uniqueID[pid_read];
+             expansion_fac[pid_write] = expansion_fac[pid_read];
+             pid_write++; 
+           }
+          pid_read++;
+        }
+      snapshot_split[s_id] = pid_write - previous_write; 
+    }  
+  boxed->NumPart = pid_write; 
+  delete preselect;
+  cout << "Num accepted here = " << pid_write <<  " / input = " << pid_read << endl; 
+}
+
+void makeBoxFromSimulation(SimulationLoader *loader, SimuData* &boxed, MetricFunctor metric, prepSimulation_info& args_info)
+{
+  vector<long> targets, split;
+  long previous_target_num = 0;
+  PreselectParticles *preselect = 0;
+
+  if (args_info.resubsample_given)
+    {
+      preselect = new PreselectParticles(args_info.resubsample_arg, args_info.resubsample_seed_arg);
+      preselect->reset();
+    }
+
+  for (int nf = 0; nf < loader->num_files(); nf++)
+    {
+      SimuData *simu;
+      double *expfact;
+
+      cout << format("Analyzing and selecting targets in file number %d / %d") % (nf+1) % loader->num_files() << endl;
+      simu = loader->loadFile(nf);      
+
+      metric(simu, 0);
+
+      selectBox(simu, targets, args_info, preselect);
+      split.push_back(targets.size() - previous_target_num);
+      previous_target_num = targets.size();
+      
+      delete simu;
+    }
+
+  createBox(loader->getHeader(), targets, split, boxed, args_info);
+
+  if (preselect)
+    delete preselect;
+}
+
+int main(int argc, char **argv)
+{
+  prepSimulation_info args_info;
+  prepSimulation_conf_params args_params;
+  SimuData *simu, *simuOut;
+  SimulationLoader *loader;
+ 
+  prepSimulation_conf_init(&args_info);
+  prepSimulation_conf_params_init(&args_params);
+  
+  args_params.check_required = 0;
+  if (prepSimulation_conf_ext (argc, argv, &args_info, &args_params))
+    return 1;
+  
+  if (!args_info.configFile_given)
+    {
+      if (prepSimulation_conf_required (&args_info, PREPSIMULATION_CONF_PACKAGE))
+        return 1;
+    }
+  else
+    {
+      args_params.check_required = 1;
+      args_params.initialize = 0;
+      if (prepSimulation_conf_config_file (args_info.configFile_arg,
+					 &args_info,
+					 &args_params))
+	return 1;
+    }
+  
+  prepSimulation_conf_print_version();
+
+  SimulationPreprocessor *preselector = new PreselectParticles(args_info.subsample_arg, args_info.subsample_seed_arg);
+
+  if (args_info.ramsesBase_given || args_info.ramsesId_given)
+    {
+      if (args_info.ramsesBase_given && args_info.ramsesId_given)   {
+	loader = ramsesLoader(args_info.ramsesBase_arg, 
+			      args_info.ramsesId_arg, 
+			      true,  // double precision with ramses... set this to false if you are dealing with single precision
+			      NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID, preselector);
+      }
+      else
+	{
+	  cerr << "Both ramsesBase and ramsesId are required to be able to load snapshots" << endl;
+	  return 1;
+	}
+    }
+  else if (args_info.gadget_given)
+    {
+      loader = gadgetLoader(args_info.gadget_arg, 1/args_info.gadgetUnit_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID, 1, preselector);
+    }
+  else if (args_info.gadget2_given)
+    {
+      loader = gadgetLoader(args_info.gadget2_arg, 1/args_info.gadgetUnit_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID, 2, preselector);
+    }
+  else if (args_info.flash_given)
+    {
+      loader = flashLoader(args_info.flash_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID, preselector); 
+  }
+#ifdef SDF_SUPPORT
+  else if (args_info.multidark_given)
+    {
+      loader = multidarkLoader(args_info.multidark_arg, preselector);      
+    }
+  else if (args_info.sdf_given)
+    {
+      loader = sdfLoader(args_info.sdf_arg, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID, args_info.sdf_splitting_arg, preselector);
+    }
+#endif
+  else
+    {
+      cerr << "A simulation snapshot is required to generate a mock catalog." << endl;
+      return 1;
+    }
+     
+  if (loader == 0)
+    {
+      cerr << "Error while loading " << endl;
+      return 1;
+    }
+  simu = loader->getHeader();
+
+  {
+    SimuData *header = loader->getHeader();
+    cout << "Hubble = " << header->Hubble << endl;
+    cout << "Boxsize = " << header->BoxSize << endl;
+    cout << "Omega_M = " << header->Omega_M << endl;
+    cout << "Omega_Lambda = " << header->Omega_Lambda << endl;
+    cout << "Subsample fraction: " << (args_info.subsample_given ? args_info.subsample_arg : 1.0) << endl;
+  }
+  double *expfact;
+
+  boost::function2<void, SimuData*, double*> metricOperation=
+    boost::bind(metricTransform, _1, args_info.axis_arg, args_info.preReShift_flag,
+                      args_info.peculiarVelocities_flag, _2, 
+                      args_info.cosmo_flag);
+
+  if (args_info.inputParameter_given)
+    makeBoxFromParameter(loader->getHeader(), simuOut, args_info);
+  else
+    makeBoxFromSimulation(loader, simuOut, metricOperation, args_info);
+
+  // Reset the random number generator
+  preselector->reset();
+
+  long loaded = 0;
+  for (int nf = 0; nf < loader->num_files(); nf++)
+    {
+      long num_targets = simuOut->as<long>("snapshot_split")[nf];
+      cout << format("Building box from particles in %d / %d") % (nf+1) % loader->num_files() << endl;
+     
+      if (num_targets == 0)
+        {
+          cout << "No particles selected there. Skipping." << endl;
+          continue;
+        }
+      
+      SimuData *simu = loader->loadFile(nf);
+      double *efac = new double[simu->NumPart];
+      metricOperation(simu, efac);
+      buildBox(simu, num_targets, loaded, simuOut, efac);
+      
+      loaded += num_targets;
+      assert(loaded <= simuOut->NumPart);
+
+      delete simu;
+
+      delete[] efac;
+    }
+
+  if (args_info.joggleParticles_flag)
+    joggleParticles(simuOut);
+ 
+  saveBox(simuOut, args_info.outputParameter_arg, args_info);
+  generateOutput(simuOut, args_info.axis_arg, 
+                 args_info.output_arg);
+  delete preselector;
+ 
+  double subsample = 1.0;
+  if (args_info.subsample_given) subsample = args_info.subsample_arg;
+  if (args_info.resubsample_given) subsample = args_info.resubsample_arg;
+  printf("Done! %5.2e\n", subsample); 
+  return 0;
+}

c_source/prep/prepSimulation.ggo

@@ -0,0 +1,45 @@
+package "prepSimulation"
+version "0"
+
+option "configFile"	 -	"Configuration file path"	string optional
+
+# Ramses data
+option "ramsesBase" -	"Base directory for ramses" string optional
+option "ramsesId"   -   "Ramses snapshot id" int optional
+
+option "gadget"	    -	"Base name of gadget snapshot (without parallel writing extension)" string optional
+option "gadget2"	-	"Base name of gadget snapshot (version 2, without parallel writing extension)" string optional
+option "flash"	    -   "Base name for FLASH snapshot" string optional
+option "multidark"	    -   "Base name for multidark snapshot" string optional
+
+option "sdf"        -   "SDF snapshot name" string optional
+option "sdf_splitting" - "Number of artificial splitting of the SDF snapshot" int optional default="20"
+
+option "axis"	    -	"Redshift axis (X=0, Y=1, Z=2)" int optional default="2"
+
+option "output"	    -   "Output filename for particles" string required
+option "outputParameter" - "Output geometry parameter file for postprocessing" string required
+
+option "rangeX_min"   -   "Minimum range in X for making the box" double required
+option "rangeX_max"   -   "Maximum range in X for making the box" double required
+option "rangeY_min"   -   "Minimum range in Y for making the box" double required
+option "rangeY_max"   -   "Maximum range in Y for making the box" double required
+option "rangeZ_min"   -   "Minimum range in Z for making the box (after distortion)" double required
+option "rangeZ_max"   -   "Maximum range in Z for making the box (after distortion)" double required
+
+option "preReShift" - "Reshift the zero of the Z axis" flag off
+option "peculiarVelocities" - "Added peculiar velocities distortion" flag off
+
+option "cosmo"	- "Apply cosmological redshift" flag off
+
+option "subsample" - "Subsample the input simulation by the specified amount" double optional default="1.0"
+
+option "inputParameter" - "Input geometry (optional, warning!)" string optional
+option "gadgetUnit" - "Unit of length in gadget file in Mpc/h" double optional default="0.001"
+
+option "subsample_seed" - "Seed for random number generation to select the subsample" int optional default="190524"
+
+option "resubsample" - "Resubsampling factor compared to the subsampled simulation" double optional 
+option "resubsample_seed" - "Seed for resubsampling from a subsampled simulation" int optional default="20132011"
+
+option "joggleParticles"	- "Slightly joggle the input particle positions" flag off

c_source/pruning/CMakeLists.txt

@@ -0,0 +1,8 @@
+include_directories(${CMAKE_CURRENT_BINARY_DIR})
+
+SET(pruneVoids_SRCS pruneVoids.cpp)
+add_genopt(pruneVoids_SRCS pruneVoids.ggo pruneVoids_conf STRUCTNAME pruneVoids_info)
+add_executable(pruneVoids ${pruneVoids_SRCS})
+target_link_libraries(pruneVoids ${ZOB_LIBS})
+
+install(TARGETS pruneVoids DESTINATION ${VIDE_BIN})

c_source/pruning/pruneVoids.ggo

@@ -0,0 +1,44 @@
+package "removeEdgeVoids"
+version "alpha"
+
+
+option "configFile" - "Configuration filename" string optional
+
+option "partFile" - "Particle file from generateFromCatalog" string required
+
+option "extraInfo" - "Extra particle file from generateFromCatalog" string required
+
+option "voidDesc" - "Void list from ZOBOV" string required
+
+option "void2Zone" - "Zone file from ZOBOV" string required
+
+option "partVol" - "Particle volume file from ZOBOV" string required
+
+option "partAdj" - "Adjacency file from ZOBOV" string required
+
+option "zone2Part" - "Particle file from ZOBOV" string required
+option "mockIndex" - "Beginning index of mock particles" int required
+
+option "numVoids" - "Number of voids" int required
+
+option "isObservation" - "We are working with observational data" flag off
+
+option "useComoving" - "Void positions are in comoving coordinates" flag off
+
+option "omegaM" - "Omega_M for redshift convertion" double optional default="0.27"
+
+option "zMin" - "Minimum redshift of sample" double optional default="0.0"
+option "zMax" - "Maximum redshift of sample" double optional default="10.0"
+
+option "rMin" - "Minimum allowable void radius" double optional default="0.0"
+
+option "outputDir" - "Directory to place outputs" string required
+option "sampleName" - "unique string to assign to outputs" string required
+
+option "periodic" - "Set of edges which are periodic" string optional default="xy"
+
+option "tolerance" - "Fraction of void width to consider edge" double optional default="1.0"
+
+option "centralRadFrac" - "Fraction of void radii to consider central region" double optional default="4"
+
+option "maxCentralDen" - "Maximum central density to accept as a void" double optional default="0.2" 

c_source/util/CMakeLists.txt

@@ -0,0 +1,5 @@
+
+add_library(zobovTool loadZobov.cpp particleInfo.cpp contour_pixels.cpp)
+
+set_target_properties(zobovTool PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS})
+target_link_libraries(zobovTool ${COSMOTOOL_LIBRARY} ${GSL_LIBRARIES} ${NETCDF_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${HDF5_HL_LIBRARIES} ${HDF5_LIBRARIES} ${DL_LIBRARY})

c_source/util/contour_pixels.cpp

@@ -0,0 +1,107 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/contour_pixels.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <vector>
+#include <healpix_map.h>
+#include <healpix_map_fitsio.h>
+#include "contour_pixels.hpp"
+
+using namespace std;
+
+static const bool DEBUG = true;
+
+void computeFilledPixels(Healpix_Map<float>& m, vector<int>& filled)
+{
+  filled.clear();
+  for (int p = 0; p < m.Npix(); p++)
+    if (m[p] > 0)
+      filled.push_back(p);
+}
+
+void computeContourPixels(Healpix_Map<float>& m, vector<int>& contour)
+{
+  contour.clear();
+  for (int p = 0; p < m.Npix(); p++)
+    {
+      fix_arr<int, 8> result;
+
+      m.neighbors(p, result);
+      for (int q = 0; q < 8; q++)
+	{
+	  if (result[q] < 0)
+	    continue;
+
+	  float delta = (m[p]-0.5)*(m[result[q]]-0.5);
+	  if (delta < 0)
+	    {
+	      contour.push_back(p);
+	      // This is boundary go to next pixel
+	      break;
+	    }
+	}      
+    }
+
+  if (DEBUG)
+    {
+      Healpix_Map<int> contour_map;
+
+      contour_map.SetNside(m.Nside(), RING);
+      contour_map.fill(0);
+      for (int p = 0; p < contour.size(); p++)
+	{
+	  contour_map[contour[p]]=1;
+	}
+
+      fitshandle h;
+      h.create("!contour_map.fits");
+      write_Healpix_map_to_fits(h, contour_map, planckType<int>());
+    }
+}
+
+void computeMaskPixels(Healpix_Map<float>& m, vector<int>& contour)
+{
+  for (int p = 0; p < m.Npix(); p++)
+    {
+
+    if (m[p]>0)
+      {
+        contour.push_back(p);
+        // This is boundary go to next pixel
+      }
+    }
+
+  if (DEBUG)
+    {
+      Healpix_Map<int> contour_map;
+
+      contour_map.SetNside(m.Nside(), RING);
+      contour_map.fill(0);
+      for (int p = 0; p < contour.size(); p++)
+  {
+    contour_map[contour[p]]=1;
+  }
+
+      fitshandle h;
+      h.create("!mask_map.fits");
+      write_Healpix_map_to_fits(h, contour_map, planckType<int>());
+    }
+}
+

c_source/util/contour_pixels.hpp

@@ -0,0 +1,33 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/contour_pixels.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef  __CONTOUR_PIXELS_HPP
+#define  __CONTOUR_PIXELS_HPP
+
+#include <vector>
+#include <healpix_map.h>
+
+
+void computeContourPixels(Healpix_Map<float>& map, std::vector<int>& contour);
+void computeFilledPixels(Healpix_Map<float>& map, std::vector<int>& contour);
+void computeMaskPixels(Healpix_Map<float>& map, std::vector<int>& contour);
+
+#endif

c_source/util/gslIntegrate.hpp

@@ -0,0 +1,54 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/gslIntegrate.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef __MYGSL_INTEGRATE_HPP
+#define __MYGSL_INTEGRATE_HPP
+
+#include <gsl/gsl_integration.h>
+
+template<typename FunT>
+double gslSpecialFunction(double x, void *param)
+{
+  FunT *f = (FunT *)param;
+
+  return (*f)(x);
+}
+
+template<typename FunT>
+double gslIntegrate(FunT& v, double a, double b, double prec, int NPTS = 1024)
+{
+  gsl_integration_workspace *w = gsl_integration_workspace_alloc(NPTS);
+  gsl_function f;
+  double result;
+  double abserr;
+
+  f.function = &gslSpecialFunction<FunT>;
+  f.params = &v;
+
+  gsl_integration_qag(&f, a, b, prec, 0, NPTS, GSL_INTEG_GAUSS61,
+		      w, &result, &abserr);
+  
+  gsl_integration_workspace_free(w);
+
+  return result;
+}
+
+#endif

c_source/util/loadZobov.cpp

@@ -0,0 +1,245 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/loadZobov.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#include <cassert>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <cstdlib>
+#include <sstream>
+#include <algorithm>
+#include "loadZobov.hpp"
+#include <CosmoTool/fortran.hpp>
+
+using namespace std;
+using namespace CosmoTool;
+
+bool loadZobov(const char *descName, const char *adjName, const char *voidsName, 
+	       const char *volName, ZobovRep& z)
+{
+  ifstream descFile(descName);
+  ifstream adjFile(adjName);
+  ifstream volFile(voidsName);
+  int32_t numParticles, numZones, numPinZone;
+  int32_t totalParticles;
+  int32_t numVoids;
+  int32_t minParticlesInZone, maxParticlesInZone;
+
+  adjFile.read((char *)&numParticles, sizeof(numParticles));
+  adjFile.read((char *)&numZones, sizeof(numZones));
+  if (!adjFile)
+    return false;
+
+  cout << "Number of particles = " << numParticles << endl;
+  cout << "Number of zones = " << numZones << endl;
+
+  totalParticles = 0;
+
+  minParticlesInZone = -1;
+  maxParticlesInZone = -1;
+
+  z.allZones.resize(numZones);
+  for (int zone = 0; zone < numZones; zone++)
+    {
+      adjFile.read((char *)&numPinZone, sizeof(numPinZone));
+      if (!adjFile)
+	{
+	  cout << "Problem on the zone " << zone << " / " << numZones << endl;
+	  return false;
+	}
+      
+      z.allZones[zone].pId.resize(numPinZone);
+      adjFile.read((char *)&z.allZones[zone].pId[0], sizeof(int)*numPinZone);
+
+      if (maxParticlesInZone < 0 || numPinZone > maxParticlesInZone)
+	maxParticlesInZone = numPinZone;
+
+      if (minParticlesInZone < 0 || numPinZone < minParticlesInZone)
+	minParticlesInZone = numPinZone;
+
+      totalParticles += numPinZone;      
+    }
+  cout << "Zoned " << totalParticles << endl;
+
+  cout << "Minimum number of particles in zone = " << minParticlesInZone << endl;
+  cout << "Maximum number of particles in zone = " << maxParticlesInZone << endl;
+  
+  if (totalParticles != numParticles)
+    {
+      cerr << "The numbers of particles are inconsistent ! (" << totalParticles << " vs " << numParticles << ")"<< endl;
+      abort();
+    }
+
+  volFile.read((char *)&numVoids, sizeof(numVoids));
+  if (!volFile)
+    return false;
+
+  cout << "Number of voids = " << numVoids << endl;
+
+  z.allVoids.resize(numVoids);
+  for (int v = 0; v < numVoids; v++)
+    {
+      int32_t numZinV;
+
+      volFile.read((char *)&numZinV, sizeof(numZinV));
+      if (!volFile)
+	return false;
+
+      z.allVoids[v].zId.resize(numZinV);
+
+      int *zId = new int[numZinV];      
+
+      volFile.read((char *)zId, sizeof(int)*numZinV);
+      for (int k = 0; k < numZinV; k++)
+	z.allVoids[v].zId[k] = zId[k];
+      std::sort(&z.allVoids[v].zId[0], &z.allVoids[v].zId[numZinV]);
+
+      delete[] zId;
+    }
+
+  if (volName != 0)
+    {
+      cout << "Loading particle volumes (requested)" << endl;
+      ifstream f(volName);
+      int numParticles;
+
+      if (!f)
+	{
+	  cerr << "No such file " << volName << endl;
+	  abort();
+	}
+
+      f.read((char *)&numParticles, sizeof(int));
+      z.particleVolume.resize(numParticles);
+      f.read((char *)&z.particleVolume[0], sizeof(float)*numParticles);
+    }
+
+  cout << "Loading description" << endl;
+  int lineid = 1;
+
+  string line;
+  getline(descFile, line);
+  lineid++;
+  getline(descFile, line);
+  lineid++;
+  getline(descFile, line);
+  lineid++;
+  while (!descFile.eof())
+    {
+      istringstream lineStream(line.c_str());
+      int orderId, volId, coreParticle, numParticlesInZone, numZonesInVoid, numInVoid;
+      float coreDensity, volumeZone, volumeVoid, densityContrast;
+      float probability;
+
+      lineStream
+	>> orderId
+	   >> volId
+	   >> coreParticle
+	   >> coreDensity
+	   >> volumeZone
+	   >> numParticlesInZone
+	   >> numZonesInVoid
+	   >> volumeVoid
+	   >> numInVoid
+	   >> densityContrast
+	   >> probability;
+      if (!lineStream)
+	{
+	  cerr << "Error in text stream at line " << lineid << endl;
+	  abort();
+	}
+    
+/*  
+    sscanf(line.c_str(), "%d %d %d %f %f %d %d %f %d %f %f\n", &orderId, &volId,
+           &coreParticle, &coreDensity, &volumeZone, &numParticlesInZone, 
+           &numZonesInVoid,
+           &volumeVoid, &numInVoid, &densityContrast, &probability);
+*/
+
+      z.allVoids[volId].proba = probability;
+      z.allVoids[volId].volume = volumeVoid;
+      z.allVoids[volId].numParticles = numInVoid;
+      z.allVoids[volId].coreParticle = coreParticle;
+      
+      // Sanity check
+      int actualNumber = 0;
+      for (int j = 0; j < z.allVoids[volId].zId.size(); j++)
+        {
+          int zzid = z.allVoids[volId].zId[j];
+          assert(zzid < z.allZones.size());
+	  actualNumber += z.allZones[zzid].pId.size();
+        }
+      
+      if (actualNumber != numInVoid)
+	{
+	  cerr << "Sanity check failed."
+	       << " The number of particles in the description (" 
+	       << numInVoid
+	       << ") is different from the one in the file ("
+	       << actualNumber << ")" << endl;
+	}
+      getline(descFile, line);
+      lineid++;
+    }
+
+  cout << "Done loading" << endl;
+
+
+  
+
+  return true;
+}
+
+bool loadZobovParticles(const char *fname, std::vector<ZobovParticle>& particles)
+{
+  UnformattedRead f(fname);
+  int N;
+
+  f.beginCheckpoint();
+  N = f.readInt32();
+  f.endCheckpoint();
+
+  particles.resize(N);
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].x = f.readReal32();
+    }
+  f.endCheckpoint();
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].y = f.readReal32();
+    }
+  f.endCheckpoint();
+
+  f.beginCheckpoint();
+  for (int i = 0; i < N; i++)
+    {
+      particles[i].z = f.readReal32();
+    }
+  f.endCheckpoint();
+
+  return true;
+}

c_source/util/loadZobov.hpp

@@ -0,0 +1,60 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/loadZobov.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef __LOAD_ZOBOV_HPP
+#define __LOAD_ZOBOV_HPP
+
+#include <vector>
+
+struct ZobovZone
+{
+  std::vector<int> pId;
+};
+
+struct ZobovVoid
+{
+  std::vector<int> zId;
+  float proba;
+  int numParticles, coreParticle;
+  float volume;
+  float barycenter[3];
+  float nearestBoundary;
+};
+
+struct ZobovRep
+{
+  std::vector<ZobovZone> allZones;
+  std::vector<ZobovVoid> allVoids;
+  std::vector<float> particleVolume;
+};
+
+struct ZobovParticle
+{
+  float x, y, z;
+};
+
+bool loadZobov(const char *descName,
+	       const char *adjName, const char *voidName,
+	       const char *volName, ZobovRep& z);
+
+bool loadZobovParticles(const char *fname, std::vector<ZobovParticle>& particles);
+
+#endif

c_source/util/particleInfo.cpp

@@ -0,0 +1,137 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/particleInfo.cpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+#include "particleInfo.hpp"
+
+#include <CosmoTool/fortran.hpp>
+#include <cstdlib>
+#include <netcdf>
+
+using namespace std;
+using namespace CosmoTool;
+using namespace netCDF;
+
+template <bool failure>
+double retrieve_attr_safe_double(NcFile& f, const char* name, double defval) {
+  NcGroupAtt a = f.getAtt(name);
+  if (a.isNull()) {
+    if (failure) abort();
+    return defval;
+  }
+  if (a.getAttLength() != 1) {
+    abort();
+  }
+  double x;
+  a.getValues(&x);
+  return x;
+}
+
+template <bool failure>
+int retrieve_attr_safe_int(NcFile& f, const char* name, int defval) {
+  NcGroupAtt a = f.getAtt(name);
+  if (a.isNull()) {
+    if (failure) abort();
+    return defval;
+  }
+  if (a.getAttLength() != 1) {
+    abort();
+  }
+  int x;
+  a.getValues(&x);
+  return x;
+}
+
+bool loadParticleInfo(ParticleInfo& info, const std::string& particles,
+                      const std::string& extra_info) {
+  int numpart;
+  int isObservation;
+
+  try {
+    NcFile f_info(extra_info, NcFile::read);
+
+    info.ranges[0][0] =
+        retrieve_attr_safe_double<true>(f_info, "range_x_min", 0);
+    info.ranges[0][1] =
+        retrieve_attr_safe_double<true>(f_info, "range_x_max", 0);
+    info.ranges[1][0] =
+        retrieve_attr_safe_double<true>(f_info, "range_y_min", 0);
+    info.ranges[1][1] =
+        retrieve_attr_safe_double<true>(f_info, "range_y_max", 0);
+    info.ranges[2][0] =
+        retrieve_attr_safe_double<true>(f_info, "range_z_min", 0);
+    info.ranges[2][1] =
+        retrieve_attr_safe_double<true>(f_info, "range_z_max", 0);
+    info.mask_index = retrieve_attr_safe_int<true>(f_info, "mask_index", 0);
+    isObservation = retrieve_attr_safe_int<false>(f_info, "is_observation", 0);
+
+    for (int i = 0; i < 3; i++)
+      info.length[i] = info.ranges[i][1] - info.ranges[i][0];
+
+    try {
+      UnformattedRead f(particles);
+
+      float mul, offset;
+
+      f.beginCheckpoint();
+      numpart = f.readInt32();
+      f.endCheckpoint();
+
+      info.particles.resize(numpart);
+
+      offset = info.ranges[0][0];
+      // TEST PMS NON-COBIC BOXES
+      // mul = 1.0;
+      mul = info.ranges[0][1] - info.ranges[0][0];
+      f.beginCheckpoint();
+      for (int i = 0; i < numpart; i++)
+        info.particles[i].x = mul * f.readReal32();
+      f.endCheckpoint();
+
+      offset = info.ranges[1][0];
+      // mul = 1.0;
+      mul = info.ranges[1][1] - info.ranges[1][0];
+      f.beginCheckpoint();
+      for (int i = 0; i < numpart; i++)
+        info.particles[i].y = mul * f.readReal32();
+      f.endCheckpoint();
+
+      offset = info.ranges[2][0];
+      // mul = 1.0;
+      mul = info.ranges[2][1] - info.ranges[2][0];
+      f.beginCheckpoint();
+      for (int i = 0; i < numpart; i++)
+        info.particles[i].z = mul * f.readReal32();
+      f.endCheckpoint();
+
+      if (!isObservation) {
+        for (int i = 0; i < numpart; i++) {
+          info.particles[i].x += info.ranges[0][0];
+          info.particles[i].y += info.ranges[1][0];
+          info.particles[i].z += info.ranges[2][0];
+        }
+      }
+    } catch (NoSuchFileException const& e) {
+      return false;
+    }
+
+  } catch (exceptions::NcCantRead const&) {
+    return false;
+  }
+  return true;
+}

c_source/util/particleInfo.hpp

@@ -0,0 +1,46 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/particleInfo.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef _PARTICLE_INFO_HEADER_HPP
+#define _PARTICLE_INFO_HEADER_HPP
+
+#include <vector>
+#include <string>
+
+struct ParticleData {
+  float x, y, z;
+};
+
+typedef std::vector<ParticleData> ParticleVector;
+
+struct ParticleInfo
+{
+  ParticleVector particles;
+  float ranges[3][2];
+  float length[3];
+  int   mask_index; // PMS
+};
+
+bool loadParticleInfo(ParticleInfo& info,
+		      const std::string& particles, 
+		      const std::string& extra_info);
+
+#endif

c_source/util/voidTree.hpp

@@ -0,0 +1,367 @@
+/*+
+    VIDE -- Void IDentification and Examination -- ./c_tools/libzobov/voidTree.hpp
+    Copyright (C) 2010-2014 Guilhem Lavaux
+    Copyright (C) 2011-2014 P. M. Sutter
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; version 2 of the License.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
++*/
+
+
+
+#ifndef _VOID_TREE_HPP
+#define _VOID_TREE_HPP
+
+#include <iostream>
+#include <stdint.h>
+#include "loadZobov.hpp"
+#include <list>
+#include <set>
+#include <vector>
+
+struct VoidNode
+{
+  int vid;
+  VoidNode *parent;
+  std::list<VoidNode *> children;
+};
+
+struct VoidNodeOnDisk
+{
+  int vid;
+  int parent;
+};
+
+class VoidTree
+{
+protected:
+  uint32_t totalNumNodes, activeNodes;
+  VoidNode *nodes;
+  VoidNode *rootNode;  
+  ZobovRep& zobov;
+public:
+  typedef std::list<VoidNode *> VoidList;
+
+  void dumpTree(std::ostream& o)
+  {
+    VoidNodeOnDisk data;
+
+    o.write((char*)&activeNodes, sizeof(uint32_t));
+    for (uint32_t i = 0; i < activeNodes; i++)
+      {
+	data.vid = nodes[i].vid;
+	if (nodes[i].parent == 0)
+	  data.parent = -1;
+	else
+	  data.parent = nodes[i].parent - nodes;
+	o.write((char *)&data, sizeof(VoidNodeOnDisk));
+      }
+    
+  }  
+
+  int lookupParent(int voidId, const std::vector<std::list<int> >& voids_for_zones)
+  {
+    int lastSize = 0x7fffffff;
+    int goodParent = -1;
+    ZobovVoid &ref_void = zobov.allVoids[voidId];
+    const std::list<int>& candidateList = voids_for_zones[ref_void.zId.front()];
+    std::list<int>::const_iterator iter_candidate = candidateList.begin();
+
+    //    std::cout << "candidate list size is " << candidateList.size() << std::endl;
+
+    while (iter_candidate != candidateList.end())
+      {
+        int vid_candidate = *iter_candidate;
+
+        if (vid_candidate == voidId)
+          break;
+        ++iter_candidate;
+      }
+    if (iter_candidate == candidateList.end())
+      {
+//        std::cout << "Failure to lookup parent" << std::endl;
+        return -1;
+      }
+
+    // voidId must be in the list.
+    //    assert(iter_candidate != candidateList.end());
+    // Go back
+
+    iter_candidate = candidateList.end();
+
+    int vid_good_candidate = -1;
+    int old_good_candidate_size = zobov.allZones.size()+1;
+
+    do
+      {
+        int vid_candidate;
+
+        --iter_candidate;
+
+        vid_candidate = *iter_candidate;
+        std::vector<int>& candidate_zIds = zobov.allVoids[vid_candidate].zId;
+
+	if (voidId == vid_candidate)
+	  continue;
+
+        if (candidate_zIds.size() < ref_void.zId.size())
+          {
+            continue;
+          }
+
+	int counter = 0;
+        // All zones id are sorted in each void. So we just have parse the 
+        // vector of zones and check whether all the zones in ref_void.zId is
+        // in iter_candidate->zId, the list is analyzed only once.
+        // THOUGHT: candidateList may contain directly the information. It would suffice to have the void ids sorted according to volume. Then we just have to jump to the indice just smaller than voidId.
+
+        int k = 0;
+        for (int j = 0; j < candidate_zIds.size() && k < ref_void.zId.size(); j++)
+          {
+            if (candidate_zIds[j] == ref_void.zId[k])
+              k++;
+            else if (candidate_zIds[j] > ref_void.zId[k])
+              break;
+          }
+        if (k==ref_void.zId.size())
+	  {
+	    if (candidate_zIds.size() < old_good_candidate_size)
+	      {
+		vid_good_candidate = vid_candidate;
+		old_good_candidate_size = candidate_zIds.size();
+	      }
+	    //	    std::cout << "Found parent " << vid_candidate << std::endl;
+	    //	    return vid_candidate;
+	  }
+	
+        // Go bigger, though I would say we should not to.
+      }
+    while (iter_candidate != candidateList.begin()) ;
+    //if (vid_good_candidate < 0)
+    //  std::cout << "Failure to lookup parent (2)" << std::endl;
+    return vid_good_candidate;
+  }
+
+  VoidTree(ZobovRep& rep, std::istream& disk)
+    : zobov(rep)
+  {
+    totalNumNodes = rep.allVoids.size();
+
+    disk.read((char *)&activeNodes, sizeof(uint32_t));
+    nodes = new VoidNode[activeNodes];
+    rootNode = 0;
+    for (uint32_t i = 0; i < activeNodes; i++)
+      {
+	VoidNodeOnDisk data;
+
+	disk.read((char *)&data, sizeof(data));	
+	nodes[i].vid = data.vid;
+	if (data.parent < 0)
+	  {
+	    if (rootNode != 0)
+	      {
+		std::cerr << "Multiple root to the tree !!!" << std::endl;
+		abort();
+	      }
+	    nodes[i].parent = 0;
+	    rootNode = &nodes[i];
+	  }
+	else
+	  {
+	    nodes[i].parent = nodes + data.parent;
+	    nodes[i].parent->children.push_back(&nodes[i]);
+	  }
+      }
+    
+    computeMaxDepth();
+    computeChildrenByNode();
+  }
+  
+  VoidTree(ZobovRep& rep)
+    : zobov(rep)
+  {
+    totalNumNodes = rep.allVoids.size();
+
+    std::vector<std::list<int> > voids_for_zones;
+
+    voids_for_zones.resize(rep.allZones.size());
+    for (int i = 0; i < rep.allVoids.size(); i++) 
+       {
+         ZobovVoid& v = rep.allVoids[i];
+         for (int j = 0; j < v.zId.size(); j++)
+           voids_for_zones[v.zId[j]].push_back(i);
+       }
+
+    // One additional for the mega-root
+    nodes = new VoidNode[totalNumNodes+1];
+
+    for (int i = 0; i <= totalNumNodes; i++)
+      {
+	nodes[i].vid = i;
+	nodes[i].parent = 0;
+      }
+
+    std::cout << "Linking voids together..." << std::endl;
+    double volMin = 0;// 4*M_PI/3*pow(4.*512/500.,3);
+    int inserted = 0;
+    for (int i = 0; i < rep.allVoids.size(); i++)
+      {
+        if (rep.allVoids[i].volume < volMin) continue;
+
+	int p = lookupParent(i, voids_for_zones);
+        if ((i % 1000) == 0) std::cout << i << std::endl;
+
+	if (p >= 0)
+	  {
+	    nodes[p].children.push_back(&nodes[i]);
+	    nodes[i].parent = &nodes[p];
+	  }
+	  
+	inserted++;
+      }
+
+    assert(inserted <= totalNumNodes); 
+    rootNode = &nodes[inserted];
+    rootNode->vid = -1;
+    rootNode->parent = 0;
+
+    for (int i = 0; i < inserted; i++)
+      if (nodes[i].parent == 0)
+	{
+	  nodes[i].parent = rootNode;
+	  rootNode->children.push_back(&nodes[i]);
+	}
+    activeNodes = inserted+1;
+    
+    computeMaxDepth();
+    computeChildrenByNode();
+  }
+
+  ~VoidTree()
+  {
+    delete[] nodes;
+  }
+
+  int _depth_computer(VoidNode *node)
+  {
+    VoidList::iterator i = node->children.begin();
+    int d = 0;
+
+    while (i != node->children.end())
+      {
+	d = std::max(d,_depth_computer(*i));
+	++i;
+      }  
+
+    return d+1;
+  }
+
+  void computeMaxDepth()
+  {
+    std::cout << "maximum depth is " <<  _depth_computer(rootNode) << std::endl;
+  }
+
+  struct _children_stat {
+    int num, min_num, max_num, num_zero,num_one, num_multiple;
+  };
+
+  void _children_computer(VoidNode *node, _children_stat& s)
+  {
+    VoidList::iterator i = node->children.begin();
+    int d = 0, j = 0;
+
+    while (i != node->children.end())
+      {
+        _children_computer(*i, s);
+        ++i;
+        ++j;
+      }
+    s.num += j;
+    if (j!= 0)
+      s.min_num = std::min(s.min_num, j);
+    else s.num_zero ++;
+    if (j==1) s.num_one++;
+    if (j>1) s.num_multiple++;
+    s.max_num = std::max(s.max_num, j);
+  }
+
+  void computeChildrenByNode()
+  {
+    _children_stat s;
+   s.num = 0;
+   s.min_num = activeNodes+1;
+   s.max_num = s.num_zero = s.num_one =s.num_multiple= 0;
+   _children_computer(rootNode, s);
+    std::cout << "Average children by node " << s.num*1.0/activeNodes << " , " << s.min_num << " " << s.max_num << " " << s.num_zero << " " << s.num_one << " " << s.num_multiple << std::endl;
+  }
+
+  int getParent(int vid) const
+  {
+    assert(nodes[vid].parent != 0);
+    return nodes[vid].parent->vid;
+  }
+  
+  const VoidList& getChildren(int vid) const
+  {
+    return nodes[vid].children;
+  }
+
+
+  VoidNode *getRoot() { return rootNode; }
+
+  template<typename T>
+  void walkNode(VoidNode *node, T& traverse)
+  {
+    if (!traverse(node))
+      return;
+
+    VoidList::iterator i = node->children.begin();
+    
+    while (i != node->children.end())
+      {
+	walkNode(*i, traverse);
+	++i;
+      }
+  }
+
+  template<typename T>
+  void walk(T& traverse)
+  {
+    walkNode(rootNode, traverse);
+  }
+
+  template<typename T, typename T2>
+  void walkNodeWithMark(VoidNode *node, T& traverse, const T2& mark)
+  {
+    T2 new_mark = mark;
+
+    if (!traverse(node, new_mark))
+      return;
+
+    VoidList::iterator i = node->children.begin();
+    
+    while (i != node->children.end())
+      {
+	walkNodeWithMark(*i, traverse, new_mark);
+	++i;
+      }    
+  }
+
+  template<typename T,typename T2>
+  void walkWithMark(T& traverse, T2 mark)
+  {
+    walkNodeWithMark(rootNode, traverse, mark);
+  }
+};
+
+#endif

c_source/zobov/CMakeLists.txt

@@ -0,0 +1,16 @@
+
+add_executable(voz1b1 voz1b1.c readfiles.c vozutil.c voz.h)
+target_link_libraries(voz1b1 ${QHULL_LIBRARY} ${MATH_LIB})
+
+
+#add_executable(jozov jozov.c findrtop.c)
+#target_link_libraries(jozov ${MATH_LIB})
+
+add_executable(vozinit vozinit.c readfiles.c)
+target_link_libraries(vozinit ${MATH_LIB})
+
+
+add_executable(voztie voztie.c readfiles.c)
+target_link_libraries(voztie ${MATH_LIB})
+
+install(TARGETS voz1b1 vozinit voztie DESTINATION ${VIDE_BIN})

c_source/zobov/findrtop.c

@@ -0,0 +1,81 @@
+#include <math.h>
+#include <stdlib.h>
+
+/*------------------------------------------------------------------------------
+  Find nb elements of Real array a having the largest value.
+  Returns index iord of these elements, ordered so iord[0] corresponds
+  to element a[iord[0]] having the largest value.
+  If nb > na, then the last nb-na elements of iord are undefined.
+  Elements of a that are equal are left in their original order.
+
+  Courtesy Andrew Hamilton.
+*/
+void findrtop(double *a, int na, int *iord, int nb)
+{
+#undef	ORDER
+#define	ORDER(ia, ja)	a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja)
+
+  int i, ia, ib, it, n;
+  
+  n = (na <= nb)? na : nb;
+  if (n <= 0) return;
+  
+  /* heap first n elements, so smallest element is at top of heap */
+  for (ib = (n >> 1); ib < n; ib++) {
+    iord[ib] = ib;
+  }
+  for (ia = (n >> 1) - 1; ia >= 0; ia--) {
+    i = ia;
+    for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+      if (ib+1 < n) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(ia, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = ia;
+  }
+  
+  /* now compare rest of elements of array to heap */
+  for (ia = n; ia < na; ia++) {
+    /* if new element is greater than smallest, sift it into heap */
+    i = 0;
+    if (ORDER(ia, iord[i])) {
+      for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) {
+	if (ib+1 < n) {
+	  if (ORDER(iord[ib], iord[ib+1])) ib++;
+	}
+	if (ORDER(ia, iord[ib])) {
+	  iord[i] = iord[ib];
+	  i = ib;
+	} else {
+	  break;
+	}
+      }
+      iord[i] = ia;
+    }
+  }
+  
+  /* unheap iord so largest element is at top */
+  for (ia = n - 1; ia > 0; ia--) {
+    it = iord[ia];
+    i = 0;
+    iord[ia] = iord[i];
+    for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) {
+      if (ib+1 < ia) {
+	if (ORDER(iord[ib], iord[ib+1])) ib++;
+      }
+      if (ORDER(it, iord[ib])) {
+	iord[i] = iord[ib];
+	i = ib;
+      } else {
+	break;
+      }
+    }
+    iord[i] = it;
+  }
+}

c_source/zobov/jozov.c

@@ -0,0 +1,580 @@
+#include <assert.h>
+/* jovoz.c by Mark Neyrinck */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#define BIGFLT 1e30 /* Biggest possible floating-point number */
+#define NLINKS 1000 /* Number of possible links with the same rho_sl */
+#define FF fflush(stdout)
+
+typedef struct Particle {
+  float dens;
+  int nadj;
+  int ncnt;
+  int *adj;
+} PARTICLE;
+
+typedef struct Zone {
+  int core; /* Identity of peak particle */
+  int np; /* Number of particles in zone */
+  int npjoin; /* Number of particles in the joined void */
+  int nadj; /* Number of adjacent zones */
+  int nhl; /* Number of zones in final joined void */
+  float leak; /* Volume of last leak zone*/
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+  float denscontrast; /* density contrast */
+  double vol; /* Total volume of all particles in the zone */
+  double voljoin; /* Total volume of all particles in the joined void */
+} ZONE;
+
+typedef struct ZoneT {
+  int nadj; /* Number of zones on border */
+  int *adj; /* Each adjacent zone, with ... */
+  float *slv; /* Smallest Linking Volume */
+} ZONET;
+
+void findrtop(double *a, int na, int *iord, int nb);
+
+int main(int argc,char **argv) {
+
+  FILE *adj, *vol, *zon, *zon2, *txt;
+  PARTICLE *p;
+  ZONE *z;
+  ZONET *zt;
+  char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile;
+  int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2;
+  int *jumper, *jumped, *numinh;
+  int *zonenum, *zonelist, *zonelist2;
+  int link[NLINKS], link2, nl;
+  float lowvol, voltol, prob;
+
+  int q,q2;
+
+  int za, nin;
+  int testpart;
+  char already, interior, *inyet, *inyet2, added, beaten;
+  int *nm, **m;
+
+  float maxvol, minvol;
+  double *sorter, e1,maxdenscontrast;
+  int *iord;
+
+  int mockIndex;
+
+  e1 = exp(1.)-1.;
+
+  if (argc != 8) {
+    printf("Wrong number of arguments.\n");
+    printf("arg1: adjacency file\n");
+    printf("arg2: volume file\n");
+    printf("arg3: output file containing particles in each zone\n");
+    printf("arg4: output file containing zones in each void\n");
+    printf("arg5: output text file\n");
+    printf("arg6: Density threshold (0 for no threshold)\n");
+    printf("arg7: Beginning index of mock galaxies\n\n");
+    exit(0);
+  }
+  adjfile = argv[1];
+  volfile = argv[2];
+  zonfile = argv[3];
+  zonfile2 = argv[4];
+  txtfile = argv[5];
+  if (sscanf(argv[6],"%f",&voltol) == 0) {
+    printf("Bad density threshold.\n");
+    exit(0);
+  }
+  if (sscanf(argv[7],"%d",&mockIndex) == 0) {
+    printf("Bad mock galaxy index.\n");
+    exit(0);
+  }
+  printf("TOLERANCE: %f\n", voltol);
+  if (voltol <= 0.) {
+    printf("Proceeding without a density threshold.\n");
+    voltol = 1e30;
+  }
+
+  adj = fopen(adjfile, "r");
+  if (adj == NULL) {
+    printf("Unable to open %s\n",adjfile);
+    exit(0);
+  }
+  fread(&np,1, sizeof(int),adj);
+  if (mockIndex < 0)
+    mockIndex = np;
+  
+  printf("adj: %d particles\n", np);
+  FF;
+
+  p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
+  /* Adjacencies*/
+  for (i=0;i<np;i++) {
+    fread(&p[i].nadj,1,sizeof(pid_t),adj); 
+    /* The number of adjacencies per particle */
+    if (p[i].nadj > 0)
+      p[i].adj = (pid_t *)malloc(p[i].nadj*sizeof(pid_t));
+    else p[i].adj = 0;
+    p[i].ncnt = 0; /* Temporarily, it's an adj counter */
+  }
+  for (i=0;i<np;i++) {
+    fread(&nin,1,sizeof(pid_t),adj);
+    if (nin > 0)
+      for (k=0;k<nin;k++) {
+	fread(&j,1,sizeof(pid_t),adj);
+	if (j < np) {
+	  /* Set both halves of the pair */
+          assert(i < j);
+          if (p[i].ncnt == p[i].nadj)
+            {
+              int q;
+              printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", i, j);
+              for (q=0;q<p[i].nadj;q++)
+                printf("  %d\n", p[i].adj[q]);
+              //abort();
+            }
+         if (p[j].ncnt == p[j].nadj)
+            {
+              int q;
+              printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", j, i);
+              for (q=0;q<p[j].nadj;q++)
+                printf("  %d\n", p[j].adj[q]);
+              //abort();
+           }
+
+          p[i].adj[p[i].ncnt] = j;
+          p[j].adj[p[j].ncnt] = i;
+	  p[i].ncnt++; p[j].ncnt++;
+	} else {
+	  printf("%d: adj = %d\n",i,j);
+	}
+      }
+  }
+  fclose(adj);
+
+  /* Check that we got all the pairs */
+  /*  adj = fopen(adjfile, "r");
+      fread(&np,1, sizeof(int),adj);*/
+  for (i=0;i<np;i++) {
+    /*    fread(&nin,1,sizeof(int),adj); /* actually nadj */
+    // PMS
+    if (p[i].ncnt != p[i].nadj && i < mockIndex) {
+      /*if (p[i].ncnt != p[i].nadj) {*/
+    // END PMS
+      p[i].nadj = p[i].ncnt;
+      printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
+      /*exit(0);*/
+    }
+  }
+//  fclose(adj);
+
+  /* Volumes */
+  vol = fopen(volfile, "r");
+  if (vol == NULL) {
+    printf("Unable to open volume file %s.\n\n",volfile);
+    exit(0);
+  }
+  fread(&np2,1, sizeof(int),adj);
+  if (np != np2) {
+    printf("Number of particles doesn't match! %d != %d\n",np,np2);
+    exit(0);
+  }
+  printf("%d particles\n", np);
+  FF;
+  for (i=0;i<np;i++) {
+    fread(&p[i].dens,1,sizeof(float),vol);
+    // PMS
+    if ((p[i].dens < 1e-30) || (p[i].dens > 1e30) && i < mockIndex) {
+    //if ((p[i].dens < 1e-30) || (p[i].dens > 1e30)) {
+    // END PMS
+      printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens);
+      p[i].dens = 1.;
+    }
+    p[i].dens = 1./p[i].dens; /* Get density from volume */
+  }
+  fclose(vol);
+
+  jumped = (pid_t *)malloc(np*sizeof(pid_t));  
+  jumper = (pid_t *)malloc(np*sizeof(pid_t));
+  numinh = (pid_t *)malloc(np*sizeof(pid_t));
+
+  /* find jumper */
+  for (i = 0; i < np; i++) {
+    minvol = p[i].dens; jumper[i] = -1;
+    for (j=0; j< p[i].nadj; j++) {
+      if (p[p[i].adj[j]].dens < minvol) {
+	jumper[i] = p[i].adj[j];
+	minvol = p[jumper[i]].dens;
+      }
+    }
+    numinh[i] = 0;
+  }
+
+  printf("About to jump ...\n"); FF;
+  
+  /* Jump */
+  for (i = 0; i < np; i++) {
+    jumped[i] = i;
+    while (jumper[jumped[i]] > -1)
+      jumped[i] = jumper[jumped[i]];
+    numinh[jumped[i]]++;
+  }
+  printf("Post-jump ...\n"); FF;
+  
+  nzones = 0;
+  for (i = 0; i < np; i++)
+    if (numinh[i] > 0) nzones++;
+  printf("%d initial zones found\n", nzones);
+
+  z = (ZONE *)malloc(nzones*sizeof(ZONE));
+  if (z == NULL) {
+    printf("Unable to allocate z\n");
+    exit(0);
+  }
+  zt = (ZONET *)malloc(nzones*sizeof(ZONET));
+  if (zt == NULL) {
+    printf("Unable to allocate zt\n");
+    exit(0);
+  }
+  zonenum = (int *)malloc(np*sizeof(int));
+  if (zonenum == NULL) {
+    printf("Unable to allocate zonenum\n");
+    exit(0);
+  }
+
+  h = 0;
+  for (i = 0; i < np; i++)
+    if (numinh[i] > 0) {
+      z[h].core = i;
+      zonenum[i] = h;
+      h++;
+    } else {
+      zonenum[i] = -1;
+    }
+ 
+  /* Count border particles */
+  for (i = 0; i < np; i++)
+    for (j = 0; j < p[i].nadj; j++) {
+      testpart = p[i].adj[j];
+      if (jumped[i] != jumped[testpart])
+	zt[zonenum[jumped[i]]].nadj++;
+    }
+  
+  for (h=0;h<nzones;h++) {
+    zt[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
+    if (zt[h].adj == NULL) {
+      printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
+      exit(0);
+    }
+    zt[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
+    if (zt[h].slv == NULL) {
+      printf("Unable to allocate %d slv's of zone %d\n",zt[h].nadj,h);
+      exit(0);
+    }
+    zt[h].nadj = 0;
+  }
+
+  /* Find "weakest links" */
+  for (i = 0; i < np; i++) {
+    h = zonenum[jumped[i]];
+    for (j = 0; j < p[i].nadj; j++) {
+      testpart = p[i].adj[j];
+      if (h != zonenum[jumped[testpart]]) {
+	if (p[testpart].dens > p[i].dens) {
+	  /* there could be a weakest link through testpart */
+	  already = 0;
+	  for (za = 0; za < zt[h].nadj; za++)
+	    if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
+	      already = 1;
+	      if (p[testpart].dens < zt[h].slv[za]) {
+		zt[h].slv[za] = p[testpart].dens;
+	      }
+	    }
+	  if (already == 0) {
+	    zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
+	    zt[h].slv[zt[h].nadj] = p[testpart].dens;
+	    zt[h].nadj++;
+	  }
+	} else { /* There could be a weakest link through i */
+	  already = 0;
+	  for (za = 0; za < zt[h].nadj; za++)
+	    if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
+	      already = 1;
+	      if (p[i].dens < zt[h].slv[za]) {
+		zt[h].slv[za] = p[i].dens;
+	      }
+	    }
+	  if (already == 0) {
+	    zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
+	    zt[h].slv[zt[h].nadj] = p[i].dens;
+	    zt[h].nadj++;
+	  }
+	}
+      }
+    }
+  }
+  printf("Found zone adjacencies\n"); FF;
+
+  /* Free particle adjacencies */
+  for (i=0;i<np; i++) if (p[i].adj != 0) free(p[i].adj);
+
+  /* Use z instead of zt */
+  for (h=0;h<nzones;h++) {
+    /*printf("%d ",zt[h].nadj);*/
+    z[h].nadj = zt[h].nadj;
+    z[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
+    z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
+    for (za = 0; za<zt[h].nadj; za++) {
+      z[h].adj[za] = zt[h].adj[za];
+      z[h].slv[za] = zt[h].slv[za];
+    }
+    free(zt[h].adj);
+    free(zt[h].slv);
+    z[h].np = numinh[z[h].core];
+  }
+  free(zt);
+  free(numinh);
+
+  m = (pid_t **)malloc(nzones*sizeof(pid_t *));
+  /* Not in the zone struct since it'll be freed up (contiguously, we hope)
+     soon */
+  nm = (pid_t *)malloc(nzones*sizeof(pid_t));
+  for (h=0; h<nzones; h++) {
+    m[h] = (pid_t *)malloc(z[h].np*sizeof(pid_t));
+    nm[h] = 0;
+    z[h].vol = 0.;
+  }
+  for (i=0; i<np; i++) {
+    h = zonenum[jumped[i]];
+    if (i == z[h].core) {
+      m[h][nm[h]] = m[h][0];
+      m[h][0] = i; /* Put the core particle at the top of the list */
+    } else {
+      m[h][nm[h]] = i;
+    }
+    z[h].vol += 1.0/(double)p[i].dens;
+    nm[h] ++;
+  }
+  free(nm);
+
+
+  zon = fopen(zonfile,"w");
+  if (zon == NULL) {
+    printf("Problem opening zonefile %s.\n\n",zonfile);
+    exit(0);
+  }
+  fwrite(&np,1,sizeof(pid_t),zon);
+  fwrite(&nzones,1,sizeof(int),zon);
+  for (h=0; h<nzones; h++) {
+// PMS
+  //printf("%d %d %d", &(z[h].np), m[h], z[h].np);
+// END PMS
+    fwrite(&(z[h].np),1,sizeof(pid_t),zon);
+    fwrite(m[h],z[h].np,sizeof(pid_t),zon);
+    free(m[h]);
+  }
+  free(m);
+  close(zon);
+
+  inyet = (char *)malloc(nzones*sizeof(char));
+  inyet2 = (char *)malloc(nzones*sizeof(char));
+  zonelist = (int *)malloc(nzones*sizeof(int));
+  zonelist2 = (int *)malloc(nzones*sizeof(int));
+  sorter = (double *)malloc((nzones+1)*sizeof(double));
+
+  for (h = 0; h< nzones; h++) {
+    inyet[h] = 0;
+    inyet2[h] = 0;
+  }
+
+  nhl = 0; 
+
+  maxvol = 0.;
+  minvol = BIGFLT;
+  maxdenscontrast = 0.;
+  for(i=0;i<np; i++){
+    if (p[i].dens > maxvol) maxvol = p[i].dens;
+    if (p[i].dens < minvol) minvol = p[i].dens;
+  }
+  printf("Densities range from %e to %e.\n",minvol,maxvol);FF;
+
+  zon2 = fopen(zonfile2,"w");
+  if (zon2 == NULL) {
+    printf("Problem opening zonefile %s.\n\n",zonfile2);
+    exit(0);
+  }
+  fwrite(&nzones,1,sizeof(int),zon2);
+
+  for (h = 0; h<nzones; h++) {
+    nhlcount = 0;
+    for (hl = 0; hl < nhl; hl++)
+      inyet[zonelist[hl]] = 0;
+
+    zonelist[0] = h;
+    inyet[h] = 1;
+    nhl = 1;
+    z[h].npjoin = z[h].np;
+    do {
+      /* Find the lowest-volume (highest-density) adjacency */
+      lowvol = BIGFLT; nl = 0; beaten = 0;
+      for (hl = 0; hl < nhl; hl++) {
+	h2 = zonelist[hl];
+	if (inyet[h2] == 1) { /* If it's not already identified as 
+				 an interior zone, with inyet=2 */
+	  interior = 1;
+	  for (za = 0; za < z[h2].nadj; za ++) {
+	    if (inyet[z[h2].adj[za]] == 0) {
+	      interior = 0;
+	      if (z[h2].slv[za] == lowvol) {
+		link[nl] = z[h2].adj[za];
+		nl ++;
+		if (nl == NLINKS) {
+		  printf("Too many links with the same rho_sl!  Increase NLINKS from %d\n",nl);
+		  exit(0);
+		}
+	      }
+	      if (z[h2].slv[za] < lowvol) {
+		lowvol = z[h2].slv[za];
+		link[0] = z[h2].adj[za];
+		nl = 1;
+	      }
+	    }
+	  }
+	  if (interior == 1) inyet[h2] = 2; /* No bordering exter. zones */
+	}
+      }
+
+      if (nl == 0) {
+	beaten = 1;
+	z[h].leak = maxvol;
+	continue;
+      }
+	
+      if (lowvol > voltol) {
+	beaten = 1;
+	z[h].leak = lowvol;
+	continue;
+      }
+
+      for (l=0; l < nl; l++)
+	if (p[z[link[l]].core].dens < p[z[h].core].dens)
+	  beaten = 1;
+      if (beaten == 1) {
+	z[h].leak = lowvol;
+	continue;
+      }
+      /* Add everything linked to the link(s) */
+      nhl2 = 0;
+      for (l=0; l < nl; l++) {
+	if (inyet2[link[l]] == 0) {
+	  zonelist2[nhl2] = link[l];
+	  inyet2[link[l]] = 1;
+	  nhl2 ++;
+	  added = 1;
+	  while ((added == 1) && (beaten == 0)) {
+	    added = 0;
+	    for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) {
+	      h2 = zonelist2[hl];
+	      if (inyet2[h2] == 1) {
+		interior = 1; /* Guilty until proven innocent */
+		for (za = 0; za < z[h2].nadj; za ++) {
+		  link2 = z[h2].adj[za];
+		  if ((inyet[link2]+inyet2[link2]) == 0) {
+		    interior = 0;
+		    if (z[h2].slv[za] <= lowvol) {
+		      if (p[z[link2].core].dens < p[z[h].core].dens) {
+			beaten = 1;
+			break;
+		      }
+		      zonelist2[nhl2] = link2;
+		      inyet2[link2] = 1;
+		      nhl2++;
+		      added = 1;
+		    }
+		  }
+		}
+		if (interior == 1) inyet2[h2] = 2;
+	      }
+	    }
+	  }
+	}
+      }
+      for (hl = 0; hl < nhl2; hl++)
+	inyet2[zonelist2[hl]] = 0;
+      
+      /* See if there's a beater */
+      if (beaten == 1) {
+	z[h].leak = lowvol;
+      } else {
+	for (h2 = 0; h2 < nhl2; h2++) {
+	  zonelist[nhl] = zonelist2[h2];
+	  inyet[zonelist2[h2]] = 1;
+	  nhl++;
+	  z[h].npjoin += z[zonelist2[h2]].np;
+	}
+      }
+      if (nhl/10000 > nhlcount) {
+	nhlcount = nhl/10000;
+	printf(" %d",nhl); FF;
+      }
+    } while((lowvol < BIGFLT) && (beaten == 0));
+    
+    z[h].denscontrast = z[h].leak/p[z[h].core].dens;
+    if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
+    
+    /* find biggest denscontrast */
+    if (z[h].denscontrast > maxdenscontrast) {
+      maxdenscontrast = (double)z[h].denscontrast;
+    }
+
+    /* Don't sort; want the core zone to be first */
+
+    if (nhlcount > 0) { /* Outputs the number of zones in large voids */
+      printf(" h%d:%d\n",h,nhl);
+      FF;
+    }
+    /* Calculate volume */
+    z[h].voljoin = 0.;
+    for (q = 0; q<nhl; q++) {
+      z[h].voljoin += z[zonelist[q]].vol;
+    }
+
+    z[h].nhl = nhl;
+
+    fwrite(&nhl,1,sizeof(int),zon2);
+    fwrite(zonelist,nhl,sizeof(int),zon2);
+  }
+  fclose(zon2);
+
+  printf("Maxdenscontrast = %f.\n",maxdenscontrast);
+
+  /* Assign sorter by probability (could use volume instead) */
+  for (h=0; h< nzones; h++)
+    sorter[h] = (double)z[h].denscontrast;
+    
+  /* Text output file */
+
+  printf("about to sort ...\n");FF;
+
+  iord = (int *)malloc(nzones*sizeof(int));
+
+  findrtop(sorter, nzones, iord, nzones);
+
+  txt = fopen(txtfile,"w");
+  fprintf(txt,"%d particles, %d voids.\n", np, nzones);
+  fprintf(txt,"Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
+  for (h=0; h<nzones; h++) {
+    i = iord[h];
+    prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
+// PMS
+    if (z[i].np == 1) continue;
+// END PMS
+    fprintf(txt,"%d %d %d %e %e %d %d %e %d %f %6.2e\n",
+	    h+1, i, z[i].core, p[z[i].core].dens, z[i].vol, z[i].np, z[i].nhl, z[i].voljoin, z[i].npjoin, z[i].denscontrast, prob);
+
+  } /* h+1 to start from 1, not zero */
+  fclose(txt);
+
+  printf("Done!\n");FF;
+  return(0);
+}

c_source/zobov/readfiles.c

@@ -0,0 +1,121 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#define DL for (d=0;d<3;d++) /* Dimension loop */
+#define BF 1e30
+
+/* Positions */
+/* Returns number of particles read */
+int posread(char *posfile, float ***p, float fact) {
+
+  FILE *pos;
+  int np,dum,d,i;
+  float xmin,xmax,ymin,ymax,zmin,zmax;
+  float *ptemp;
+
+  pos = fopen(posfile, "r");
+  if (pos == NULL) {
+    printf("Unable to open position file %s\n\n",posfile);
+    exit(0);
+  }
+  /* Fortran77 4-byte headers and footers */
+  /* Delete "dum" statements if you don't need them */
+
+  /* Read number of particles */
+  fread(&dum,1,4,pos); fread(&np,1, sizeof(int),pos); fread(&dum,1,4,pos);
+
+  /* Allocate the arrays */
+  (*p) = (float **)malloc(np*sizeof(float *));
+  ptemp = (float *)malloc(np*sizeof(float));
+
+  printf("np = %d\n",np);
+
+  /* Fill the arrays */
+  fread(&dum,1,4,pos); 
+  fread(ptemp,np,4,pos);
+  for (i=0; i<np; i++) {
+    (*p)[i] = (float *)malloc(3*sizeof(float));
+    if ((*p)[i] == NULL) {
+      printf("Unable to allocate particle array in readfiles!\n");
+      fflush(stdout);
+      exit(0);
+    }
+    (*p)[i][0] = ptemp[i];
+  }
+  fread(&dum,1,4,pos); 
+  fread(&dum,1,4,pos); 
+  fread(ptemp,np,4,pos);
+  for (i=0; i<np; i++) (*p)[i][1] = ptemp[i];
+  fread(&dum,1,4,pos); 
+  fread(&dum,1,4,pos); 
+  fread(ptemp,np,4,pos);
+  for (i=0; i<np; i++) (*p)[i][2] = ptemp[i];
+  fread(&dum,1,4,pos); 
+
+  fclose(pos);
+  free(ptemp);
+
+  /* Get into physical units (Mpc/h) */
+  
+
+  printf("%f\n",fact);fflush(stdout);
+  for (i=0; i<np; i++) DL (*p)[i][d] *= fact;
+
+  /* Test range -- can comment out */
+  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+  for (i=0; i<np;i++) {
+    if ((*p)[i][0]<xmin) xmin = (*p)[i][0]; if ((*p)[i][0]>xmax) xmax = (*p)[i][0];
+    if ((*p)[i][1]<ymin) ymin = (*p)[i][1]; if ((*p)[i][1]>ymax) ymax = (*p)[i][1];
+    if ((*p)[i][2]<zmin) zmin = (*p)[i][2]; if ((*p)[i][2]>zmax) zmax = (*p)[i][2];
+  }
+  printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
+
+  return(np);
+}
+
+/* Velocities */
+/* Returns number of particles read */
+/* Used in voboz, but not zobov, which doesn't use velocities */
+int velread(char *velfile, float ***v, float fact) {
+
+  FILE *vel;
+  int np,dum,d,i;
+  float xmin,xmax,ymin,ymax,zmin,zmax;
+
+  vel = fopen(velfile, "r");
+  if (vel == NULL) {
+    printf("Unable to open velocity file %s\n\n",velfile);
+    exit(0);
+  }
+  /* Fortran77 4-byte headers and footers */
+  /* Delete "dum" statements if you don't need them */
+
+  /* Read number of particles */
+  fread(&dum,1,4,vel); fread(&np,1, sizeof(int),vel); fread(&dum,1,4,vel);
+
+  /* Allocate the arrays */
+  (*v) = (float **)malloc(3*sizeof(float*));
+  for (i=0;i<3;i++) (*v)[i] = (float *)malloc(np*sizeof(float));
+
+  /* Fill the arrays */
+  fread(&dum,1,4,vel); fread((*v)[0],np,4,vel); fread(&dum,1,4,vel);
+  fread(&dum,1,4,vel); fread((*v)[1],np,4,vel); fread(&dum,1,4,vel);
+  fread(&dum,1,4,vel); fread((*v)[2],np,4,vel); fread(&dum,1,4,vel);
+
+  fclose(vel);
+
+  /* Convert from code units into physical units (km/sec) */
+  
+  for (i=0; i<np; i++) DL (*v)[d][i] *= fact;
+
+  /* Test range -- can comment out */
+  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+  for (i=0; i<np;i++) {
+    if ((*v)[0][i] < xmin) xmin = (*v)[0][i]; if ((*v)[0][i] > xmax) xmax = (*v)[0][i];
+    if ((*v)[1][i] < ymin) ymin = (*v)[1][i]; if ((*v)[1][i] > ymax) ymax = (*v)[1][i];
+    if ((*v)[2][i] < zmin) zmin = (*v)[2][i]; if ((*v)[2][i] > zmax) zmax = (*v)[2][i];
+  }
+  printf("vx: %f,%f; vy: %f,%f; vz: %f,%f\n",xmin,xmax, ymin,ymax, zmin,zmax);
+  
+  return(np);
+}

c_source/zobov/voz.h

@@ -0,0 +1,11 @@
+#define MAXVERVER 100000
+#define NGUARD 84 /*Actually, the number of SPACES between guard points
+##define NGUARD 42 /*Actually, the number of SPACES between guard points
+		    in each dim */
+
+typedef int pid_t;
+
+typedef struct Partadj {
+  pid_t nadj;
+  pid_t *adj;
+} PARTADJ;

c_source/zobov/voz1b1.c

@@ -0,0 +1,352 @@
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+
+#define DL for (d=0;d<3;d++)
+#define BF 1e30
+
+#define MAX(a,b) ( ((a) < (b)) ? (a) : (b) )
+
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
+int posread(char *posfile, float ***p, float fact);
+
+int main(int argc, char *argv[]) {
+  int exitcode;
+  pid_t i, j, np;
+  float **r;
+  coordT rtemp[3], *parts;
+  coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
+  pointT intpoints[3*MAXVERVER];
+  FILE *pos, *out;
+  char *posfile, outfile[200], *suffix, *outDir;
+  PARTADJ *adjs;
+  float *vols;
+  realT predict, xmin,xmax,ymin,ymax,zmin,zmax;
+  pid_t *orig;
+  
+  int isitinbuf;
+  char isitinmain, d;
+  int numdiv;
+  pid_t nvp, nvpall, nvpbuf;
+  realT width, width2, totwidth, totwidth2, bf, s, g;
+  float border, boxsize;
+  realT c[3];
+  int b[3];
+  double totalvol;
+
+  if (argc != 10) {
+    printf("Wrong number of arguments.\n");
+    printf("arg1: position file\n");
+    printf("arg2: border size\n");
+    printf("arg3: boxsize\n");
+    printf("arg4: suffix\n");
+    printf("arg5: number of divisions\n");
+    printf("arg6-8: b[0-2]\n\n");
+    printf("arg9: output directory\n");
+    exit(0);
+  }
+  posfile = argv[1];
+  if (sscanf(argv[2],"%f",&border) != 1) {
+    printf("That's no border size; try again.\n");
+    exit(0);
+  }
+  if (sscanf(argv[3],"%f",&boxsize) != 1) {
+    printf("That's no boxsize; try again.\n");
+    exit(0);
+  }
+  suffix = argv[4];
+  if (sscanf(argv[5],"%d",&numdiv) != 1) {
+    printf("%s is no number of divisions; try again.\n",argv[5]);
+    exit(0);
+  }
+  if (numdiv == 1) {
+    printf("Only using one division; should only use for an isolated segment.\n");
+  }
+  if (numdiv < 1) {
+    printf("Cannot have a number of divisions less than 1.  Resetting to 1.\n");
+    numdiv = 1;
+  }
+  if (sscanf(argv[6],"%d",&b[0]) != 1) {
+    printf("That's no b index; try again.\n");
+    exit(0);
+  }
+  if (sscanf(argv[7],"%d",&b[1]) != 1) {
+    printf("That's no b index; try again.\n");
+    exit(0);
+  }
+  if (sscanf(argv[8],"%d",&b[2]) != 1) {
+    printf("That's no b index; try again.\n");
+    exit(0);
+  }
+  outDir = argv[9];
+  
+  /* Boxsize should be the range in r, yielding a range 0-1 */
+  np = posread(posfile,&r,1./boxsize);
+  printf("%d particles\n",np);fflush(stdout);
+
+  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+  for (i=0; i<np;i++) {
+    if (r[i][0]<xmin) xmin = r[i][0]; if (r[i][0]>xmax) xmax = r[i][0];
+    if (r[i][1]<ymin) ymin = r[i][1]; if (r[i][1]>ymax) ymax = r[i][1];
+    if (r[i][2]<zmin) zmin = r[i][2]; if (r[i][2]>zmax) zmax = r[i][2];
+  }
+  printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
+
+  width = 1./(float)numdiv;
+  width2 = 0.5*width;
+  if (border > 0.) bf = border;
+  else bf = 0.1;
+      /* In units of 0-1, the thickness of each subregion's buffer*/
+  totwidth = width+2.*bf;
+  totwidth2 = width2 + bf;
+  
+  s = width/(float)NGUARD;
+  if ((bf*bf - 2.*s*s) < 0.) {
+    printf("bf = %f, s = %f.\n",bf,s);
+    printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
+	   sqrt(2.)*width/bf);
+    exit(0);
+  }
+  g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
+  printf("s = %f, bf = %f, g = %f.\n",s,bf,g);
+
+  fflush(stdout);
+
+  adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
+  if (adjs == NULL) {
+    printf("Unable to allocate adjs\n");
+    exit(1);
+  }
+  
+  DL c[d] = ((float)b[d])*width;
+  printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
+  /* Assign temporary array*/
+  nvpbuf = 0; /* Number of particles to tesselate, including
+		 buffer */
+  nvp = 0; /* Without the buffer */
+  for (i=0; i<np; i++) {
+    isitinbuf = 1;
+    isitinmain = 1;
+    DL {
+      rtemp[d] = (double)r[i][d] - (double)c[d];
+      if (rtemp[d] > 0.5) rtemp[d] --;
+      if (rtemp[d] < -0.5) rtemp[d] ++;
+      isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
+      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
+    }
+  
+    if (isitinbuf) nvpbuf++;
+    if (isitinmain) nvp++;
+  }
+  
+  nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard
+					points */
+
+  parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
+  orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t));
+
+  if (parts == NULL) {
+    printf("Unable to allocate parts\n");
+    fflush(stdout);
+    exit(1);
+  }
+  if (orig == NULL) {
+    printf("Unable to allocate orig\n");
+    fflush(stdout);
+    exit(1);
+  }
+
+  nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */
+  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+  for (i=0; i<np; i++) {
+    isitinmain = 1;
+    DL {
+      rtemp[d] = (realT)r[i][d] - (realT)c[d];
+      if (rtemp[d] > 0.5) rtemp[d] --;
+      if (rtemp[d] < -0.5) rtemp[d] ++;
+      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
+    }
+    if (isitinmain) {
+      parts[3*nvp] = rtemp[0];
+      parts[3*nvp+1] = rtemp[1];
+      parts[3*nvp+2] = rtemp[2];
+      orig[nvp] = i;
+      nvp++;
+      if (rtemp[0] < xmin) xmin = rtemp[0];
+      if (rtemp[0] > xmax) xmax = rtemp[0];
+      if (rtemp[1] < ymin) ymin = rtemp[1];
+      if (rtemp[1] > ymax) ymax = rtemp[1];
+      if (rtemp[2] < zmin) zmin = rtemp[2];
+      if (rtemp[2] > zmax) zmax = rtemp[2];
+    }
+  }
+  printf("nvp = %d\n",nvp);
+  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
+  nvpbuf = nvp;
+  for (i=0; i<np; i++) {
+    isitinbuf = 1;
+    isitinmain = 1;
+
+    DL {
+      rtemp[d] = (realT)r[i][d] - (realT)c[d];
+      if (rtemp[d] > 0.5) rtemp[d] --;
+      if (rtemp[d] < -0.5) rtemp[d] ++;
+      isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
+      isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
+    }
+    if (isitinbuf && !isitinmain) {
+      /*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
+	rtemp[2]*rtemp[2]));
+	printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
+      parts[3*nvpbuf] = rtemp[0];
+      parts[3*nvpbuf+1] = rtemp[1];
+      parts[3*nvpbuf+2] = rtemp[2];
+      orig[nvpbuf] = i;
+
+      nvpbuf++;
+      if (rtemp[0] < xmin) xmin = rtemp[0];
+      if (rtemp[0] > xmax) xmax = rtemp[0];
+      if (rtemp[1] < ymin) ymin = rtemp[1];
+      if (rtemp[1] > ymax) ymax = rtemp[1];
+      if (rtemp[2] < zmin) zmin = rtemp[2];
+      if (rtemp[2] > zmax) zmax = rtemp[2];
+    }
+  }
+  printf("nvpbuf = %d\n",nvpbuf);
+  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
+  nvpall = nvpbuf;
+  predict = pow(width+2.*bf,3)*(float)np;
+  printf("There should be ~ %f points; there are %d\n",predict,nvpbuf);
+
+  for (i=0;i<np;i++) free(r[i]);
+  free(r);
+  
+  /* Add guard points */
+  for (i=0; i<NGUARD+1; i++) {
+    for (j=0; j<NGUARD+1; j++) {
+      /* Bottom */
+      parts[3*nvpall]   = -width2 + (realT)i * s;
+      parts[3*nvpall+1] = -width2 + (realT)j * s;
+      parts[3*nvpall+2] = -width2 - g;
+      nvpall++;
+      /* Top */
+      parts[3*nvpall]   = -width2 + (realT)i * s;
+      parts[3*nvpall+1] = -width2 + (realT)j * s;
+      parts[3*nvpall+2] = width2 + g;
+      nvpall++;
+    }
+  }
+  for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
+    for (j=0; j<NGUARD+1; j++) {
+      parts[3*nvpall]   = -width2 + (realT)i * s;
+      parts[3*nvpall+1] = -width2 - g;
+      parts[3*nvpall+2] = -width2 + (realT)j * s;
+      nvpall++;
+      
+      parts[3*nvpall]   = -width2 + (realT)i * s;
+      parts[3*nvpall+1] = width2 + g;
+      parts[3*nvpall+2] = -width2 + (realT)j * s;
+      nvpall++;
+    }
+  }
+  for (i=0; i<NGUARD+1; i++) {
+    for (j=0; j<NGUARD+1; j++) {
+      parts[3*nvpall]   = -width2 - g;
+      parts[3*nvpall+1] = -width2 + (realT)i * s;
+      parts[3*nvpall+2] = -width2 + (realT)j * s;
+      nvpall++;
+      
+      parts[3*nvpall]   = width2 + g;
+      parts[3*nvpall+1] = -width2 + (realT)i * s;
+      parts[3*nvpall+2] = -width2 + (realT)j * s;
+      nvpall++;
+    }
+  }
+  xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+  for (i=nvpbuf;i<nvpall;i++) {
+    if (parts[3*i] < xmin) xmin = parts[3*i];
+    if (parts[3*i] > xmax) xmax = parts[3*i];
+    if (parts[3*i+1] < ymin) ymin = parts[3*i+1];
+    if (parts[3*i+1] > ymax) ymax = parts[3*i+1];
+    if (parts[3*i+2] < zmin) zmin = parts[3*i+2];
+    if (parts[3*i+2] > zmax) zmax = parts[3*i+2];
+  }
+  
+  printf("Added guard points to total %d points (should be %d)\n",nvpall,
+	 nvpbuf + 6*(NGUARD+1)*(NGUARD+1));
+  printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
+  
+  /* Do tesselation*/
+  printf("File read.  Tessellating ...\n"); fflush(stdout);
+  exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs);
+  if (exitcode != 0)
+   {
+     printf("Error while tesselating. Stopping here."); fflush(stdout);
+     exit(1);
+   }
+  
+  /* Now calculate volumes*/
+  printf("Now finding volumes ...\n"); fflush(stdout);
+  vols = (float *)malloc(nvp*sizeof(float));
+  if (vols == NULL) {
+    printf("Unable to allocate vols\n");
+    exit(1);
+  }
+  
+  for (i=0; i<nvp; i++) { /* Just the original particles
+			     Assign adjacency coordinate array*/
+    /* Volumes */
+    for (j = 0; j < adjs[i].nadj; j++)
+      DL {
+	deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d];
+	if (deladjs[3*j+d] < -0.5) deladjs[3*j+d]++;
+	if (deladjs[3*j+d] > 0.5) deladjs[3*j+d]--;
+      }
+    
+    exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
+    vols[i] *= (float)np;
+    if (i % 1000 == 0)
+      printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
+  }
+
+  /* Get the adjacencies back to their original values */
+
+  for (i=0; i<nvp; i++)
+    for (j = 0; j < adjs[i].nadj; j++)
+      adjs[i].adj[j] = orig[adjs[i].adj[j]];
+  
+  totalvol = 0.;
+  for (i=0;i<nvp; i++) {
+    totalvol += (double)vols[i];
+  }
+  printf("Average volume = %g\n",totalvol/(float)nvp);
+ 
+  /* Now the output!
+     First number of particles */
+  sprintf(outfile,"%s/part.%s.%02d.%02d.%02d",outDir,suffix,b[0],b[1],b[2]);
+
+  printf("Output to %s\n\n",outfile);
+  out = fopen(outfile,"w");
+  if (out == NULL) {
+    printf("Unable to open %s\n",outfile);
+    exit(1);
+  }
+  fwrite(&np,1, sizeof(int),out);
+  fwrite(&nvp,1, sizeof(int),out);
+  printf("nvp = %d\n",nvp);
+
+  /* Tell us where the original particles were */
+  fwrite(orig,sizeof(pid_t),nvp,out);
+  /* Volumes*/
+  fwrite(vols,sizeof(float),nvp,out);
+  /* Adjacencies */
+  for (i=0;i<nvp;i++) {
+    fwrite(&(adjs[i].nadj),1,sizeof(pid_t),out);
+    if (adjs[i].nadj > 0)
+      fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out);
+    else printf("0");
+  }
+  fclose(out);
+  
+  return(0);
+}

c_source/zobov/vozinit.c

@@ -0,0 +1,197 @@
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+
+#define NUMCPU 2
+#define DL for (d=0;d<3;d++)
+#define BF 1e30
+
+#define QHULL_MAX_PARTICLES ((1L<<24)-1)
+
+int posread(char *posfile, float ***p, float fact);
+
+int main(int argc, char *argv[]) {
+  int i, np;
+  float **rfloat, rtemp[3];
+  FILE *pos, *scr;
+  char *posfile, scrfile[200], systemstr[90], *suffix, *outDir, *vobozPath;
+  float xmin,xmax,ymin,ymax,zmin,zmax;
+  
+  int isitinbuf;
+  char isitinmain, d;
+  int numdiv;
+  int p; 
+  int nvp, nvpall, nvpbuf, nvpmin, nvpmax, nvpbufmin, nvpbufmax; /* yes, the insurance */
+  float width, width2, totwidth, totwidth2, bf, s, g;
+  float widthX, widthY, widthZ;
+  float border, boxsize, boxsizeX, boxsizeY, boxsizeZ;
+  float c[3];
+  int numGuards;
+  int b[3];
+  int numThreads;
+  int mockIndex;
+
+  if (argc != 10) {
+    printf("Wrong number of arguments.\n");
+    printf("arg1: position file\n");
+    printf("arg2: buffer size (default 0.1)\n");
+    printf("arg3: box size\n");
+    printf("arg6: number of divisions (default 2)\n");
+    printf("arg7: suffix describing this run\n");
+    printf("arg8: number of parallel threads\n");
+    printf("arg9: location of voboz executables\n");
+    printf("arg10: output directory\n");
+    printf("arg11: index of mock galaxies\n\n");
+    exit(0);
+  }
+  posfile = argv[1];
+  suffix = argv[2];
+  if (sscanf(suffix,"%f",&border) != 1) {
+    printf("That's no border size; try again.\n");
+    exit(0);
+  }
+  suffix = argv[3];
+  if (sscanf(suffix,"%f",&boxsize) != 1) {
+    printf("That's no boxsize; try again.\n");
+    exit(0);
+  }
+  suffix = argv[4];
+  if (sscanf(suffix,"%d",&numdiv) != 1) {
+    printf("That's no number of divisions; try again.\n");
+    exit(0);
+  }
+  if (numdiv < 1) {
+    printf("Cannot have a number of divisions less than 1.  Resetting to 1:\n");
+    numdiv = 1;
+  }
+
+  suffix = argv[5];
+  vobozPath = argv[6];
+  if (sscanf(vobozPath,"%d",&numThreads) != 1) {
+    printf("That's no number of threads; try again.\n");
+    exit(0);
+  }
+  vobozPath = argv[7];
+  outDir = argv[8];
+  if (sscanf(argv[9],"%d",&mockIndex) != 1) {
+    printf("That's no mock galaxy index; try again.\n");
+    exit(0);
+  }
+
+
+  /* Read the position file */
+  np = posread(posfile,&rfloat,1./boxsize);
+  /* Boxsize should be the range in r, yielding a range 0-1 */
+
+  width = boxsize/(float)numdiv;
+  //width = 1./(float)numdiv;
+  width2 = 0.5*width;
+  if (border > 0.) bf = border;
+  else bf = 0.1;
+
+  /* In units of 0-1, the thickness of each subregion's buffer*/
+  totwidth = width+2.*bf;
+  totwidth2 = width2 + bf;
+  
+  s = width/(float)NGUARD;
+  if ((bf*bf - 2.*s*s) < 0.) {
+    printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
+	   totwidth/sqrt(0.5*bf*bf));
+    printf("bf = %f\n",bf);
+    exit(0);
+  }
+  g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
+  printf("s = %f, bf = %f, g = %f.\n",s,bf,g);
+  
+  nvpmax = 0; nvpbufmax = 0; nvpmin = np; nvpbufmin = np;
+  
+  for (b[0] = 0; b[0] < numdiv; b[0]++) {
+   c[0] = ((float)b[0]+0.5)*width;
+   for (b[1] = 0; b[1] < numdiv; b[1]++) {
+    c[1] = ((float)b[1]+0.5)*width;
+    for (b[2] = 0; b[2] < numdiv; b[2]++) {
+      c[2] = ((float)b[2]+0.5)*width;
+
+      nvp = 0; /* Number of particles excluding buffer */
+      nvpbuf = 0; /* Number of particles to tesselate, including
+		     buffer */
+      xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
+      for (i=0; i<np; i++) {
+	isitinbuf = 1; isitinmain = 1;
+	for (d=0; d<3; d++) {
+	  rtemp[d] = rfloat[i][d] - c[d];
+	  if (rtemp[d] > 0.5) rtemp[d] --;
+	  if (rtemp[d] < -0.5) rtemp[d] ++;
+	  isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
+	  isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
+	}
+	if (isitinbuf) {
+	  nvpbuf++;
+	}
+	if (isitinmain) {
+	  nvp++;
+	  if (rtemp[0] < xmin) xmin = rtemp[0];
+	  if (rtemp[0] > xmax) xmax = rtemp[0];
+	  if (rtemp[1] < ymin) ymin = rtemp[1];
+	  if (rtemp[1] > ymax) ymax = rtemp[1];
+	  if (rtemp[2] < zmin) zmin = rtemp[2];
+	  if (rtemp[2] > zmax) zmax = rtemp[2];
+	}
+      }
+      if (nvp > nvpmax) nvpmax = nvp;
+      if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf;
+      if (nvp < nvpmin) nvpmin = nvp;
+      if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf;
+
+      printf("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n",
+	     b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf);
+    }
+   }
+  }
+  printf("Nvp range: %d,%d\n",nvpmin,nvpmax);
+  printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax);
+
+  numGuards = 6*(NGUARD+1)*(NGUARD+1);
+  printf("Total max particles: %d\n" , nvpbufmax+numGuards);
+  if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES)
+    {
+      printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n");
+      fflush(stdout);
+      exit(1);
+    }
+
+  /* Output script file */
+  sprintf(scrfile,"scr%s",suffix);
+  printf("Writing script file to %s.\n",scrfile);fflush(stdout);
+  scr = fopen(scrfile,"w");
+  if (scr == NULL) {
+    printf("Problem opening script file.\n");
+    fflush(stdout);
+    exit(1);
+  }
+  fprintf(scr,"#!/bin/bash -f\n");
+  p = 0;
+  for (b[0]=0;b[0]<numdiv; b[0]++) {
+   for (b[1] = 0; b[1] < numdiv; b[1]++) {
+    for (b[2] = 0; b[2] < numdiv; b[2]++) {
+      //fprintf(scr,"%s/../c_tools/zobov2/voz1b1/voz1b1_2 %s %f %f %f %f %s %d %d %d %d %d %d %s&\n",
+     //   vobozPath,
+     //   posfile,border,boxsize,boxsize,boxsize,suffix,numdiv,numdiv, numdiv,b[0],b[1],b[2],
+     //   outDir);
+      fprintf(scr,"%s/voz1b1 %s %f %f %s %d %d %d %d %s&\n",
+              vobozPath,
+              posfile,border,boxsize,suffix,numdiv,b[0],b[1],b[2], 
+              outDir);
+      p++;
+      if ((p == numThreads)) { fprintf(scr, "wait\n"); p = 0; }
+    }
+   }
+  }
+  fprintf(scr,"wait\n");
+  fprintf(scr,"%s/voztie %d %s %s %d\n", vobozPath, numdiv,suffix, outDir, mockIndex);
+  fclose(scr);
+
+  sprintf(systemstr,"chmod u+x %s",scrfile);
+  system(systemstr);
+
+  return(0);
+}

c_source/zobov/voztie.c

@@ -0,0 +1,303 @@
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "voz.h"
+
+int compare_int(const void *a, const void *b)
+{
+  const int *ia = (const int *)a;
+  const int *ib = (const int *)b;
+  if ((*ia) < (*ib))
+    return -1;
+  else if ((*ia) > (*ib))
+    return 1;
+  else
+    return 0;
+}
+
+int main(int argc, char *argv[]) {
+
+  FILE *part, *adj, *vol;
+  char partfile[200], *suffix, adjfile[200], volfile[200], *outDir;
+  float *vols, volstemp;
+  int *cnt_adj;
+  
+  PARTADJ *adjs;
+
+  int numdiv,np,np2,na;
+
+  pid_t i,j,k,p,nout;
+  pid_t nvp,npnotdone,nvpmax, nvpsum, *orig;
+  double avgnadj, avgvol;
+
+  int numRemoved = 0;
+  
+  // PMS
+  int mockIndex;
+  // END PMS
+
+  if (argc != 5) {
+    printf("Wrong number of arguments.\n");
+    printf("arg1: number of divisions (default 2)\n");
+    printf("arg2: suffix describing this run\n");
+    printf("arg3: file directory\n");
+    printf("arg4: Beginning index of mock particles\n\n");
+    exit(0);
+  }
+  if (sscanf(argv[1],"%d",&numdiv) != 1) {
+    printf("That's no number of divisions; try again.\n");
+    exit(0);
+  }
+  if (numdiv < 1) {
+    printf("Cannot have a number of divisions less than 1.  Resetting to 1:\n");
+    numdiv = 1;
+  }
+
+  suffix = argv[2];
+  outDir = argv[3];
+  if (sscanf(argv[4],"%d",&mockIndex) != 1) {
+    printf("That's no mock galaxy index; try again.\n");
+    exit(0);
+  }
+  
+  np = -1; nvpmax = -1; nvpsum = 0;
+
+  for (i = 0; i < numdiv; i++) {
+   for (j = 0; j < numdiv; j++) {
+    for (k = 0; k < numdiv; k++) {
+      sprintf(partfile,"%s/part.%s.%02d.%02d.%02d",outDir,suffix,i,j,k);
+      part = fopen(partfile,"r");
+      if (part == NULL) {
+	printf("Unable to open file %s.\n\n",partfile);
+	exit(0);
+      }
+      fread(&np2,1,sizeof(int),part);
+      fread(&nvp,1,sizeof(int),part);
+      nvpsum += nvp;
+      if (np == -1)
+	np = np2;
+      else 
+	if (np2 != np) {
+	  printf("Incompatible total particle numbers: %d,%d\n\n",np,np2);
+	  exit(0);
+	}
+      if (nvp > nvpmax) nvpmax = nvp;
+      fclose(part);
+    }
+   }
+  }
+
+  printf("We have %d particles to tie together.\n",np); fflush(stdout);
+  printf("The maximum number of particles in a file is %d.\n",nvpmax);
+
+  // PMS
+  if (mockIndex == -1) mockIndex = np;
+  // END PMS
+
+  cnt_adj = (int *)malloc(np*sizeof(int));
+  if (cnt_adj == NULL)
+    printf("Could not allocate cnt_adj.\n");
+
+  adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ));
+  if (adjs == NULL) printf("Couldn't allocate adjs.\n");
+  vols = (float *)malloc(np*sizeof(float));
+  if (vols == NULL) printf("Couldn't allocate vols.\n");
+  orig = (pid_t *)malloc(nvpmax*sizeof(pid_t));
+  if (orig == NULL) printf("Couldn't allocate orig.\n");
+  if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) {
+    printf("Not enough memory to allocate. Exiting.\n");
+    exit(0);
+  }
+  for (p=0;p<np;p++)
+    {
+      vols[p] = -1.;
+      adjs[p].nadj = 0;
+      adjs[p].adj = 0;
+    }
+
+  nvpsum = 0;
+  for (i = 0; i < numdiv; i++) {
+   for (j = 0; j < numdiv; j++) {
+    for (k = 0; k < numdiv; k++) {
+      sprintf(partfile,"%s/part.%s.%02d.%02d.%02d",outDir,suffix,i,j,k);
+      part = fopen(partfile,"r");
+      if (part == NULL) {
+	printf("Unable to open file %s.\n\n",partfile);
+	exit(0);
+      }
+      fread(&np2,1,sizeof(pid_t),part);
+      fread(&nvp,1,sizeof(pid_t),part);
+      /*printf("nvp = %d\n",nvp);fflush(stdout);*/
+
+      nvpsum += nvp;
+
+      fread(orig,nvp,sizeof(pid_t),part);
+      for (p=0;p<nvp;p++) {
+	fread(&volstemp,1,sizeof(float),part);
+	if (vols[orig[p]] > -1.)
+	  if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) {
+	    printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
+		   orig[p],vols[orig[p]],volstemp);
+	    //exit(0);
+	  }
+	vols[orig[p]] = volstemp;
+      }
+      
+      for (p=0;p<nvp;p++) {
+	fread(&na,1,sizeof(int),part);
+        pid_t pid = orig[p];
+	if (na > 0) {
+          pid_t *previous_adj = adjs[pid].adj;
+          assert(adjs[pid].nadj == 0 || adjs[pid].nadj == na);
+	  adjs[pid].nadj = na;
+          adjs[pid].adj = (pid_t *)malloc(na*sizeof(pid_t));
+	  if (adjs[pid].adj == 0) {
+	    printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
+	    exit(0);
+	  }
+	  fread(adjs[pid].adj,na,sizeof(pid_t),part);
+          if (previous_adj != 0)
+           {
+             qsort(previous_adj, na, sizeof(pid_t), compare_int);
+             qsort(adjs[pid].adj, na, sizeof(pid_t), compare_int);
+             if (memcmp(previous_adj, adjs[pid].adj, sizeof(pid_t)*na) != 0)
+               {
+                 printf("Inconsistent adjacencies between two divisions.\n");
+                 abort();
+               }
+           }
+	} else {
+	  printf("0"); fflush(stdout);
+	}
+      }
+      fclose(part);
+      printf("%d ",k);
+    }
+   }
+  }
+  printf("\n");
+
+  // PMS : remove mock galaxies and anything adjacent to a mock galaxy
+  printf("\nRemoving mock galaxies...\n"); 
+ 
+  // completely unlink mock particles
+  for (i = mockIndex; i < np; i++) {
+    vols[i] = 1.e-29;
+    adjs[i].nadj = 0;
+  }
+  
+  // unlink particles adjacent to mock galaxies
+  for (i = 0; i < mockIndex; i++) {
+    for (j = 0; j < adjs[i].nadj; j++) {
+      if (adjs[i].adj[j] > mockIndex) {
+//printf("KILLING %d\n", i);
+        vols[i] = 1.e-29;
+        adjs[i].nadj = 0;
+        numRemoved++;
+        break;
+      }
+    }
+  }
+
+  // update all other adjacencies
+  for (i = 0; i < mockIndex; i++) {
+    
+    int numAdjSaved = 0;
+    for (j = 0; j < adjs[i].nadj; j++) {
+
+      if ( adjs[adjs[i].adj[j]].nadj != 0) {
+        adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
+        numAdjSaved++;
+      }
+
+    }
+    adjs[i].nadj = numAdjSaved; 
+  }
+
+  printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex, 
+                                                                 numRemoved); 
+  printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
+
+  // END PMS
+
+  npnotdone = 0; avgnadj = 0.; avgvol = 0.;
+  for (p=0;p<np;p++) {
+    // PMS
+    if (vols[p] == 1.e-29) continue;
+    // END PMS
+    if (vols[p] == -1.) npnotdone++;
+    avgnadj += (double)(adjs[p].nadj);
+    avgvol += (double)(vols[p]);
+  }
+  if (npnotdone > 0)
+    printf("%d particles not done!\n", npnotdone);
+  printf("%d particles done more than once.\n",nvpsum-np);
+  avgnadj /= (double)np;
+  avgvol /= (double)np;
+  printf("Average # adjacencies = %lf (%f for Poisson)\n",avgnadj,
+	 48.*3.141593*3.141593/35.+2.);
+  printf("Average volume = %lf\n",avgvol);
+    
+  /* Now the output! */
+
+  sprintf(adjfile,"%s/adj%s.dat",outDir,suffix);
+  sprintf(volfile,"%s/vol%s.dat",outDir,suffix);
+
+  printf("Outputting to%s, %s\n\n", adjfile, volfile);
+
+  adj = fopen(adjfile,"w");
+  if (adj == NULL) {
+    printf("Unable to open %s\n",adjfile);
+    exit(0);
+  }
+  fwrite(&mockIndex,1, sizeof(int),adj);
+  /* Adjacencies: first the numbers of adjacencies, 
+     and the number we're actually going to write per particle */
+
+  // Recount the number of adjacencies after merge
+  for(i=0;i<mockIndex;i++)
+    cnt_adj[i] = adjs[i].nadj;
+
+  for (i=0;i<mockIndex;i++)
+    fwrite(&cnt_adj[i],1,sizeof(int),adj);
+    
+  /* Now the lists of adjacencies (without double counting) */
+  for (i=0;i<mockIndex;i++) {
+      nout = 0;
+      for (j=0;j<adjs[i].nadj; j++) 
+         if (adjs[i].adj[j] > i)
+           nout++;
+      fwrite(&nout,1,sizeof(int),adj);      
+      for (j=0;j<adjs[i].nadj; j++) {
+        pid_t id = adjs[i].adj[j];
+	if (id > i) 
+	  fwrite(&id,1,sizeof(pid_t),adj);
+      }
+    }
+
+  fclose(adj);
+  
+  /* Volumes */
+  vol = fopen(volfile,"w");
+  if (vol == NULL) {
+    printf("Unable to open %s\n",volfile);
+    exit(0);
+  }
+// PMS
+  fwrite(&mockIndex,1, sizeof(int),vol);
+  fwrite(vols,sizeof(float),mockIndex,vol);
+  //fwrite(&np,1, sizeof(int),vol);
+  //fwrite(vols,sizeof(float),np,vol);
+// END PMS
+
+  fclose(vol);
+
+  free(vols);
+  free(cnt_adj);
+  free(adjs);
+  free(orig);
+
+  return(0);
+}

c_source/zobov/vozutil.c

@@ -0,0 +1,245 @@
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+
+#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2)
+
+/*char qh_version[] = "user_eg 3.1 2001/10/04";  */
+
+int compar(const void * n1, const void * n2) {
+  int i1,i2;
+
+  i1 = *(int *)n1;
+  i2 = *(int *)n2;
+  return 2*(i1 > i2) - 1 + (i1 == i2);
+}
+
+/* Finds Delaunay adjacencies of a set of points */
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) {
+
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= stdout;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+  int i, ver, count;
+  int numfacets, numsimplicial, numridges, totneighbors, numneighbors,
+    numcoplanars, numtricoplanars;
+  setT *vertices, *vertices2, *vertex_points, *coplanar_points;
+  vertexT *vertex, **vertexp;
+  vertexT *vertex2, **vertex2p;
+  int vertex_i, vertex_n;
+  facetT *facet, **facetp, *neighbor, **neighborp;
+  pointT *point, **pointp;
+  int numdiv;
+
+  PARTADJ adjst;
+
+  int errorReported = 0;
+
+  adjst.adj = (int *)malloc(MAXVERVER*sizeof(int));
+  if (adjst.adj == NULL) {
+    printf("Unable to allocate adjst.adj\n");
+    exit(0);
+  }
+
+  /* Delaunay triangulation*/
+  sprintf (flags, "qhull s d Qt");
+  exitcode= qh_new_qhull (dim, nvpall, points, ismalloc,
+                      flags, outfile, errfile);
+
+
+  if (!exitcode) {                  /* if no error */
+    /* 'qh facet_list' contains the convex hull */
+
+    /* From qh_printvneighbors */
+    qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial,
+		   &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
+    qh_vertexneighbors();
+    vertices= qh_facetvertices (qh facet_list, NULL, 0);
+    vertex_points= qh_settemp (nvpall);
+    coplanar_points= qh_settemp (nvpall);
+    qh_setzero (vertex_points, 0, nvpall);
+    qh_setzero (coplanar_points, 0, nvpall);
+    FOREACHvertex_(vertices)
+      qh_point_add (vertex_points, vertex->point, vertex);
+    FORALLfacet_(qh facet_list) {
+      FOREACHpoint_(facet->coplanarset)
+	qh_point_add (coplanar_points, point, facet);
+    }
+    ver = 0;
+    FOREACHvertex_i_(vertex_points) {
+      (*adjs)[ver].nadj = 0;
+      if (vertex) {
+	/* Count the neighboring vertices, check that all are real
+	   neighbors */
+	adjst.nadj = 0;
+	FOREACHneighbor_(vertex) {
+	  if ((*adjs)[ver].nadj > -1) {
+	    if (neighbor->visitid) {
+	      vertices2 = neighbor->vertices;
+	      FOREACHvertex2_(vertices2) {
+		if (ver != qh_pointid(vertex2->point)) {
+		  adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point);
+		  adjst.nadj ++;
+		  if (adjst.nadj > MAXVERVER) {
+		    printf("Increase MAXVERVER to at least %d!\n",adjst.nadj);
+		    exit(0);
+		  }
+		}
+	      }
+	    } else {
+	      printf(" %d",ver);
+	      (*adjs)[ver].nadj = -1; /* There are unreal vertices here */
+	    }
+	  }
+	}
+      } else (*adjs)[ver].nadj = -2;
+
+      /* Enumerate the unique adjacencies*/
+      if (adjst.nadj >= 4) {
+	qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar);
+	count = 1;
+
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    if (adjst.adj[i] >= nvpbuf && !errorReported) {
+        errorReported = 1;
+	      printf("Guard point encountered.  Increase border and/or nguard.\n");
+	      printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2],
+		     points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]);
+	    }
+	    count++;
+	  }
+	(*adjs)[ver].adj = (int *)malloc(count*sizeof(int));
+	if ((*adjs)[ver].adj == NULL) {
+	  printf("Unable to allocate (*adjs)[ver].adj\n");
+	  exit(0);
+	}
+	(*adjs)[ver].adj[0] = adjst.adj[0];
+	count = 1;
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    (*adjs)[ver].adj[count] = adjst.adj[i];
+	    count++;
+	  }
+	(*adjs)[ver].nadj = count;
+      } else {
+	printf("Number of adjacencies %d < 4, particle %d -> %d\n",adjst.nadj,ver,ver);
+	exit(0);
+      }
+      ver++;
+      if (ver == nvp) break;
+    }
+    qh_settempfree (&coplanar_points);
+    qh_settempfree (&vertex_points);
+    qh_settempfree (&vertices);
+  }
+  qh_freeqhull(!qh_ALL);                 /* free long memory */
+  qh_memfreeshort (&curlong, &totlong);  /* free short memory and memory allocator */
+  if (curlong || totlong)
+    fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  free(adjst.adj);
+  return exitcode;
+}
+
+/* Calculates the Voronoi volume from a set of Delaunay adjacencies */
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) {
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= NULL;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  facetT *facet;	    /* set by FORALLfacets */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+
+  coordT *point, *normp, *coordp, *feasiblep, *deladj;
+  int i, j, k;
+  boolT zerodiv;
+  float runsum;
+  char region;
+  /*coordT *points;
+    pointT *intpoints;*/
+
+  /* make point array from adjacency coordinates (add offset)*/
+  /*points = (coordT *)malloc(4*numpoints*sizeof(coordT));
+  if (points == NULL) {
+    printf("Unable to allocate points\n");
+    exit(0);
+    }*/
+  for (i=0; i<numpoints; i++) {
+    runsum = 0.;
+    deladj = deladjs + 3*i;
+    point = points + 4*i;
+    for (j=0;j<3;j++) {
+      runsum += deladj[j]*deladj[j];
+      point[j] = deladj[j];
+    }
+    point[3] = -0.5*runsum;
+  }
+  sprintf (flags, "qhull H0");
+
+  exitcode= qh_new_qhull (4, numpoints, points, ismalloc,
+			  flags, outfile, errfile);
+
+  numpoints = 0;
+  if (!exitcode) {                  /* if no error */
+    FORALLfacets {
+      numpoints++;
+    }
+    /* Now we know how many points */
+    /*intpoints = (pointT *)malloc(dim*numpoints*sizeof(pointT));
+    if (intpoints == NULL) {
+    printf("Unable to allocate intpoints\n");
+    exit(0);
+    }*/
+
+    j = 0;
+    FORALLfacets {
+      if (!qh feasible_point) {
+	fprintf (stdout, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
+	qh_errexit( qh_ERRinput, NULL, NULL);
+      }
+      point= coordp= intpoints + j*3;
+      j++;
+      normp= facet->normal;
+      feasiblep= qh feasible_point;
+      if (facet->offset < -qh MINdenom) {
+	for (k= qh hull_dim; k--; )
+	  *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
+      }else {
+	for (k= qh hull_dim; k--; ) {
+	  *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
+				   &zerodiv) + *(feasiblep++);
+	  if (zerodiv) {
+	    qh_memfree (point, qh normal_size);
+	    printf("LABELprintinfinite\n");
+	    exit(0);
+	  }
+	}
+      }
+    }
+  }
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+
+  /* Now we calculate the volume */
+  sprintf (flags, "qhull FA");
+  exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc,
+			  flags, outfile, errfile);
+
+  qh_getarea(qh facet_list);
+  *vol = qh totvol;
+
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+  if (curlong || totlong)
+    fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  /*free(points); free(intpoints);*/
+
+  return exitcode;
+}

c_source/zobov/zobov_readme.txt

@@ -0,0 +1,25 @@
+This is the readme file for ZOBOV (ZOnes Bordering On Voidness),
+version 1.0 (Jan 5, 2008) a void-finding algorithm for sets of
+particles (e.g. a cosmological N-body simulation) by Mark Neyrinck.
+It is an inversion of the halo-finding algorithm VOBOZ (VOronoi BOund
+Zones), which was also written by Mark Neyrinck, with the guidance of
+Nick Gnedin and Andrew Hamilton.
+
+The gz file which once contained this file should include the
+following files:
+
+Makefile    jozov.c      readme.txt  voz1b1.c   voztie.c
+findrtop.c  readfiles.c  voz.h       vozinit.c  vozutil.c
+
+For instructions on installing and running ZOBOV, please see
+http://www.ifa.hawaii.edu/~neyrinck/zobov/.  To use ZOBOV, you will
+need to download the Qhull package, at http://www.qhull.org.  For help
+in understanding ZOBOV, please see the paper describing it
+(arXiv:0712.3049), If these resources are insufficient, feel free to
+email Mark.Neyrinck@colorado.edu with questions.
+
+This is free software.  It may be freely copied, modified, and
+redistributed, as long as ZOBOV and its author are acknowledged.
+There is no warranty or other guarantee of fitness for ZOBOV; it is
+provided "as is."  I welcome bug reports (but do not guarantee their
+fixing), which should be sent to Mark.Neyrinck@colorado.edu.