diff --git a/.gitignore b/.gitignore
index 487c453..c0ba1bd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,8 +1,4 @@
 *~
 *.o
 *.prog
-voz1b1
-jozov
-vozinit
-voztie
 .mydepend
diff --git a/c_tools/zobov2/voz1b1/voz.h b/c_tools/zobov2/voz1b1/voz.h
new file mode 100644
index 0000000..1d48e92
--- /dev/null
+++ b/c_tools/zobov2/voz1b1/voz.h
@@ -0,0 +1,29 @@
+#ifndef __VOZ_H
+#define __VOZ_H
+
+#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;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+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);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp
new file mode 100644
index 0000000..7912082
--- /dev/null
+++ b/c_tools/zobov2/voz1b1/voz1b1.cpp
@@ -0,0 +1,337 @@
+#include <iostream>
+#include <CosmoTool/miniargs.hpp>
+#include <cmath>
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+#include "voz_io.hpp"
+
+using namespace std;
+
+#define DL for (d=0;d<3;d++)
+
+bool checkParameters(int *numdiv, int *b)
+{
+  for (int i = 0; i < 3; i++)
+    {
+      if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i])
+        return false;
+    }
+  return true;
+}
+
+
+int main(int argc, char *argv[]) {
+  const float BF = std::numeric_limits<float>::max();
+
+  int exitcode;
+  pid_t i, j, np;
+  PositionData pdata;
+  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;
+  pid_t nvp, nvpall, nvpbuf;
+  realT width, width2, totwidth, totwidth2, bf, s, g;
+  double border, boxsize[3];
+  realT c[3];
+  int b[3];
+  int numdiv[3];
+  double totalvol;
+  
+  CosmoTool::MiniArgDesc args[] = {
+    { "POSITION FILE", MINIARG_STRING, &posfile },
+    { "BORDER SIZE", MINIARG_DOUBLE, &border },
+    { "BOX_X", MINIARG_DOUBLE, &boxsize[0] },
+    { "BOX_Y", MINIARG_DOUBLE, &boxsize[1] },
+    { "BOX_Z", MINIARG_DOUBLE, &boxsize[2] },
+    { "SUFFIX", MINIARG_STRING, &suffix },
+    { "NUM_DIVISION_X", MINIARG_INT, &numdiv[0] },
+    { "NUM_DIVISION_Y", MINIARG_INT, &numdiv[1] },
+    { "NUM_DIVISION_Z", MINIARG_INT, &numdiv[2] },
+    { "B0", MINIARG_INT, &b[0] },
+    { "B1", MINIARG_INT, &b[1] },
+    { "B2", MINIARG_INT, &b[2] },
+    { "OUTPUT DIRECTORY", MINIARG_STRING, &outDir },
+    { 0, MINIARG_NULL, 0 }
+  };
+
+  if (!CosmoTool::parseMiniArgs(argc, argv, args))
+    return 1;
+
+  if (!checkParameters(numdiv, b))
+    return 2;
+  
+  /* Boxsize should be the range in r, yielding a range 0-1 */
+  if (!pdata.readFrom(posfile))
+    return 3;
+
+  (cout << pdata.np << " particles" << endl).flush();
+
+  pdata.findExtrema();
+
+  (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") 
+    % pdata.np
+    % xyz_min[0] % xyz_max[0]
+    % xyz_min[1] % xyz_max[1]
+    % xyz_min[2] % xyz_max[2]).flush();
+
+
+  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(0);
+  }
+  
+  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);
+  }
+  if (orig == NULL) {
+    printf("Unable to allocate orig\n");
+    fflush(stdout);
+  }
+
+  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);
+
+  pdata.destroy();
+  
+  /* 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));
+  
+  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(0);
+  }
+  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);
+}
diff --git a/c_tools/zobov2/voz1b1/voz_io.cpp b/c_tools/zobov2/voz1b1/voz_io.cpp
new file mode 100644
index 0000000..6bf033f
--- /dev/null
+++ b/c_tools/zobov2/voz1b1/voz_io.cpp
@@ -0,0 +1,64 @@
+#include <iostream>
+#include <fstream>
+#include "voz_io.hpp"
+#include "CosmoTool/fortran.hpp"
+#include <cmath>
+
+using namespace CosmoTool;
+using namespace std;
+
+bool PositionData::readFrom(const string& fname)
+{
+  try
+    {
+      UnformattedRead f(fname.c_str());
+      
+      f.beginCheckpoint();
+      np = f.readInt32();
+      f.endCheckPoint();
+
+      for (int j = 0; j < 3; j++)
+        {
+          xyz[j] = new float[np];
+        
+          f.beginCheckpoint();
+          for (int p = 0; p < np; p++)
+            xyz[j][p] = f.readReal32();
+          f.endCheckpoint();
+        }
+    }
+  catch (const NoSuchFileException& e)
+    {
+      return false;
+    }
+  catch (const InvalidUnformattedAccess& e)
+    {
+      return false;
+    }
+  catch (const EndOfFileException& e)
+    {
+      return false;
+    }
+  
+  return true;
+}
+
+void PositionData::findExtrema()
+{
+  const float BF = std::numeric_limits<float>::max();
+
+  for (int i = 0; i < 3; i++)
+    {
+      xyz_min[i] = BF; 
+      xyz_max[i] = -BF;
+    }
+
+  for (int i = 0; i < 3; i++)
+    {
+      for (pid_t p = 0; p < np; p++) 
+        {
+          xyz_min[p] = min(xyz_min[p], xyz[p][i]);
+          xyz_max[p] = max(xyz_max[p], xyz[p][i]);
+        }
+    }
+}
diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp
new file mode 100644
index 0000000..dbeb632
--- /dev/null
+++ b/c_tools/zobov2/voz1b1/voz_io.hpp
@@ -0,0 +1,22 @@
+#ifndef __VOZ_IO_HPP
+#define __VOZ_IO_HPP
+
+#include <string>
+
+struct PositionData
+{
+  float *xyz[3];
+  pid_t np;
+  
+  void destroy()
+  {
+    for (int j = 0; j < 3; j++)
+      delete[] xyz[j];
+  }
+
+  void findExtrema();
+
+  bool readFrom(const std::string& fname);
+};
+
+#endif
diff --git a/c_tools/zobov2/voz1b1/vozutil.c b/c_tools/zobov2/voz1b1/vozutil.c
new file mode 100644
index 0000000..c56da6f
--- /dev/null
+++ b/c_tools/zobov2/voz1b1/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;
+}