added support for hyperuniform lattices using existing hard-sphere jamming routine

21 files changed, 2551 insertions, 125 deletions

diff --git a/.gitignore b/.gitignore
index 2eae92a..bebedfe 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,5 @@
 bin/*
 obj/*
 data/raw/*
+data/raw
+makefile_*
diff --git a/makefile b/makefile
index d788743..ed1142d 100644
--- a/makefile
+++ b/makefile
@@ -1,12 +1,18 @@
 
 CC = clang
+CCP = clang++
 CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp -march=native 
+CPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp -march=native 
+CCPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -shared -fopenmp -march=native 
 LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc
 
 OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand
-BIN = corr_test fracture anal_process big_anal_process cen_anal_process
+OBJP = spheres_poly/box spheres_poly/event spheres_poly/heap spheres_poly/neighbor spheres_poly/sphere
+OBJCCP = spheres_poly/spheres
+BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process
 
-all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+
+all: opt ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
 
 opt:
 	@echo build options:
@@ -14,18 +20,26 @@ opt:
 	@echo "CFLAGS   = ${CFLAGS}"
 	@echo "LDFLAGS  = ${LDFLAGS}"
 
-obj/%.o: src/%.c
+${OBJ:%=obj/%.o} ${BIN:%=obj/%.o}: obj/%.o: src/%.c
 	@echo CC -c -o $@
-	@${CC} -c -o $@ $< ${CFLAGS}
+	@${CC} -c ${CFLAGS} $< -o $@
+
+${OBJP:%=obj/%.o}: obj/%.o: src/%.cpp
+	@echo CCP -c -o $@
+	@${CCP} -c ${CPFLAGS} $< -o $@
+
+${OBJCCP:%=obj/%.o}: obj/%.o: src/%.cpp ${OBJP:%=obj/%.o} 
+	@echo ${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} $< -o $@
+	@${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} -fuse-ld=gold $< -o $@
 
-bin/%: obj/%.o ${OBJ:%=obj/%.o}
-	@echo CC -o $@
-	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
+bin/%: obj/%.o ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o}
+	@echo ${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@
+	@${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@
 
 clean:
 	@echo cleaning:
-	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+	@echo rm -f ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${OBJCCP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
+	@rm -f ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${OBJCCP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
 
 .PHONY: all clean
 
diff --git a/makefile_ept b/makefile_ept
deleted file mode 100644
index c38af25..0000000
--- a/makefile_ept
+++ /dev/null
@@ -1,31 +0,0 @@
-
-CC = clang
-CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp=libiomp5 -march=native -I/usr/include/suitesparse -I/usr/lib/gcc/x86_64-linux-gnu/5/include
-LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc
-
-OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand
-BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process
-
-all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-
-opt:
-	@echo build options:
-	@echo "CC       = ${CC}"
-	@echo "CFLAGS   = ${CFLAGS}"
-	@echo "LDFLAGS  = ${LDFLAGS}"
-
-obj/%.o: src/%.c
-	@echo CC -c -o $@
-	@${CC} -c -o $@ $< ${CFLAGS}
-
-bin/%: obj/%.o ${OBJ:%=obj/%.o}
-	@echo CC -o $@
-	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
-
-clean:
-	@echo cleaning:
-	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-
-.PHONY: all clean
-
diff --git a/makefile_hal b/makefile_hal
deleted file mode 100644
index e12d1d7..0000000
--- a/makefile_hal
+++ /dev/null
@@ -1,31 +0,0 @@
-
-CC = clang
-CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp
-LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler #-ltcmalloc
-
-OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand
-BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process
-
-all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-
-opt:
-	@echo build options:
-	@echo "CC       = ${CC}"
-	@echo "CFLAGS   = ${CFLAGS}"
-	@echo "LDFLAGS  = ${LDFLAGS}"
-
-obj/%.o: src/%.c
-	@echo CC -c -o $@
-	@${CC} -c -o $@ $< ${CFLAGS}
-
-bin/%: obj/%.o ${OBJ:%=obj/%.o}
-	@echo CC -o $@
-	@${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS}
-
-clean:
-	@echo cleaning:
-	@echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-	@rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-
-.PHONY: all clean
-
diff --git a/src/fracture.c b/src/fracture.c
index d155b02..bb7701b 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -105,8 +105,12 @@ int main(int argc, char *argv[]) {
 						lattice = SQUARE_LATTICE;
 						lattice_c = 's';
 						break;
+					case 2:
+						lattice = VORONOI_HYPERUNIFORM_LATTICE;
+						lattice_c = 'h';
+						break;
 					default:
-						printf("lattice specifier must be 0 (VORONOI_LATTICE) or 1 (SQUARE_LATTICE).\n");
+						printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (SQUARE_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
 						exit(EXIT_FAILURE);
 				}
 				break;
diff --git a/src/fracture.h b/src/fracture.h
index f4f7334..8f883bc 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -33,7 +33,8 @@
 
 typedef enum lattice_t {
 	VORONOI_LATTICE,
-	SQUARE_LATTICE
+	SQUARE_LATTICE,
+	VORONOI_HYPERUNIFORM_LATTICE
 } lattice_t;
 
 typedef enum bound_t {
@@ -185,3 +186,5 @@ bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, i
 unsigned long int rand_seed();
 
 long double rand_dist_pow(const gsl_rng *r, double beta);
+
+double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate);
diff --git a/src/graph_create.c b/src/graph_create.c
index 635b12b..b071c49 100644
--- a/src/graph_create.c
+++ b/src/graph_create.c
@@ -672,6 +672,8 @@ graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cho
 			if (bound == EMBEDDED_BOUND) side_bounds = true;
 			else side_bounds = false;
 			return ini_square_network(L, bound, side_bounds, c);
+		case (VORONOI_HYPERUNIFORM_LATTICE):
+			return ini_voro_graph(L, bound, dual, genfunc_hyperuniform, c);
 		default:
 			printf("lattice type unsupported\n");
 			exit(EXIT_FAILURE);
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
index 54d9daa..8dcd0df 100644
--- a/src/graph_genfunc.c
+++ b/src/graph_genfunc.c
@@ -20,58 +20,7 @@ double g(double rho, double dist) {
 double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) {
 	*num = 2 * pow(L / 2, 2);
 
-	// necessary to prevent crashing when underflow occurs
-	gsl_set_error_handler_off();
-
-	double *lattice = (double *)malloc(2 * (*num) * sizeof(double));
-	double rho = *num;
-	unsigned int to_gen = *num;
-	unsigned int start = 0;
-
-	if (boundary == EMBEDDED_BOUND) {
-		for (unsigned int i = 0; i < L / 2; i++) {
-			lattice[2 * i + 1] = 0;
-			lattice[2 * i] = (2. * i + 1.) / L;
-
-			lattice[L + 2 * i + 1] = 1;
-			lattice[L + 2 * i] = (2. * i + 1.) / L;
-
-			lattice[2 * L + 2 * i + 1] = (2. * i + 1.) / L;
-			lattice[2 * L + 2 * i] = 0;
-
-			lattice[3 * L + 2 * i + 1] = (2. * i + 1.) / L;
-			lattice[3 * L + 2 * i] = 1;
-		}
-
-		to_gen -= 2 * L;
-		start = 2 * L;
-	}
-
-	for (unsigned int i = 0; i < to_gen; i++) {
-		bool reject = true;
-		double x, y;
-		while(reject) {
-			x = gsl_ran_flat(r, 0, 1);
-			y = gsl_ran_flat(r, 0, 1);
-			double pp = 1;
-			for (unsigned int j = 0; j < i; j++) {
-				double ds0, ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8;
-				ds0 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2);
-				ds1 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2);
-				ds2 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2);
-				ds3 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2);
-				ds4 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2);
-				ds5 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] + 1,2);
-				ds6 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] - 1,2);
-				ds7 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] + 1,2);
-				ds8 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] - 1,2);
-				pp *= g(rho, ds0) * g(rho, ds1) * g(rho, ds2) * g(rho, ds3) * g(rho, ds4) * g(rho, ds5) * g(rho, ds6) * g(rho, ds7) * g(rho, ds8);
-			} 
-			if (pp > gsl_ran_flat(r, 0, 1)) reject = false;
-		}
-		lattice[2*start + 2 * i] = x;
-		lattice[2*start + 2 * i + 1] = y;
-	}
+	double *lattice = spheres(*num, 0.8, 0.5, 0.9, 100, 100000);
 
 	return lattice;
 }
diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp
new file mode 100644
index 0000000..716090d
--- /dev/null
+++ b/src/spheres_poly/box.cpp
@@ -0,0 +1,1157 @@
+/* 
+   Updated July 24, 2009 to include hardwall boundary condition option, as 
+   well as polydisperse spheres.  
+
+   Packing of hard spheres via molecular dynamics
+   Developed by Monica Skoge, 2006, Princeton University
+   Contact: Monica Skoge (mskoge.princeton.edu) with questions
+   This code may be used, modified and distributed freely.
+   Please cite:
+   
+   "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+   	M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+	
+   if you use these codes.	
+*/
+
+#include "box.h"
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#include <iomanip>
+
+//==============================================================
+//==============================================================
+//  Class Box: Fills box with hardspheres to given packing fraction
+//  and evolves spheres using molecular dynamics!
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+box::box(int N_i, double r_i, double growthrate_i, double maxpf_i, 
+	 double bidispersityratio_i, double bidispersityfraction_i, 
+	 double massratio_i, int hardwallBC_i):
+  r(r_i),          
+  N(N_i),
+  growthrate(growthrate_i),
+  h(N_i+1),
+  maxpf(maxpf_i),
+  bidispersityratio(bidispersityratio_i),
+  bidispersityfraction(bidispersityfraction_i),
+  massratio(massratio_i),
+  hardwallBC(hardwallBC_i)
+{
+  ngrids = Optimalngrids2(r);
+  cells.set_size(ngrids);
+ 
+  s = new sphere[N];
+  binlist = new int[N];
+  x = new vector<DIM>[N];        
+  h.s = s;
+
+  gtime = 0.;
+  rtime = 0.;
+  ncollisions = 0;
+  ntransfers = 0;
+  nchecks = 0;
+  neventstot = 0;
+  ncycles = 0;
+  xmomentum = 0.; 
+  pressure = 0.;
+  collisionrate = 0.;
+
+  cells.set_size(ngrids);
+  cells.initialize(-1);      // initialize cells to -1
+  srandom(::time(0));        // initialize the random number generator
+  for (int i=0; i<N; i++)    // initialize binlist to -1
+    binlist[i] = -1;
+  
+  time(&start);
+}
+
+
+//==============================================================
+// Destructor
+//==============================================================
+box::~box() 
+{
+  delete[] s;
+  delete[] binlist;
+  delete[] x;
+}
+
+
+//==============================================================
+// ReadFile
+//==============================================================
+void box::ReadPositions(const char* filename)
+{  
+  // open file to read in arrays
+  std::ifstream infile(filename);
+  
+  infile.ignore(256, '\n');  // ignore the dim line
+  infile.ignore(256, '\n');  // ignore the #sphere 1 line
+  infile.ignore(256, '\n');  // ignore the #sphere line
+  infile.ignore(256, '\n');  // ignore the diameter line
+  infile.ignore(1000, '\n'); // ignore the 100 010 001 line
+  infile.ignore(256, '\n');  // ignore the T T T line
+
+  for (int i=0; i<N; i++)
+    {
+      infile >> s[i].r;      // read in radius    
+      infile >> s[i].gr;     // read in growth rate
+      infile >> s[i].m;      // read in mass
+      for (int k=0; k<DIM; k++)  
+	infile >> s[i].x[k]; // read in position 
+    }
+  infile.close();
+}
+
+
+//==============================================================
+// Recreates all N spheres at random positions
+//==============================================================
+void box::RecreateSpheres(const char* filename, double temp)
+{
+  ReadPositions(filename);  // reads in positions of spheres
+  VelocityGiver(temp);      // gives spheres initial velocities
+  AssignCells();            // assigns spheres to cells
+  SetInitialEvents();
+}
+
+
+//==============================================================
+// Creates all N spheres at random positions
+//==============================================================
+void box::CreateSpheres(double temp)
+{
+  int Ncurrent = 0;
+  for(int i=0; i<N; i++)
+    {
+      CreateSphere(Ncurrent);
+      Ncurrent++;
+    }
+  if (Ncurrent != N)
+    std::cout << "problem! only made " << Ncurrent << " out of " << N << " desired spheres" << std::endl;
+  
+  VelocityGiver(temp);
+  SetInitialEvents();
+}
+
+   
+//==============================================================
+// Creates a sphere of random radius r at a random unoccupied position 
+//==============================================================
+void box::CreateSphere(int Ncurrent)
+{
+  int keeper;    // boolean variable: 1 means ok, 0 means sphere already there
+  int counter = 0;   // counts how many times sphere already exists
+  vector<DIM> xrand;  // random new position vector
+  double d = 0.;
+  double radius;
+  double growth_rate;
+  double mass;
+  int species;
+  
+  if (Ncurrent < bidispersityfraction*N)
+    {
+      radius = r;
+      growth_rate = growthrate;
+      mass = 1.;
+      species = 1;
+    }
+  else
+    {
+      radius = r*bidispersityratio;
+      growth_rate = growthrate*bidispersityratio;
+      mass = massratio;
+      species = 2;
+    }
+  
+  while (counter<1000)
+    {
+      keeper = 1;
+      
+      for(int k=0; k<DIM; k++) 
+	xrand[k] = ((double)random()/(double)RAND_MAX)*SIZE;
+      
+      for (int i=0; i<Ncurrent; i++)  // check if overlapping other spheres
+	{
+	  d=0.;
+	  if (hardwallBC)
+	    {
+	      for (int k=0; k<DIM; k++) 
+		{     
+		  d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
+		}
+	    }
+	  else
+	    {
+	      for (int k=0; k<DIM; k++)
+		{
+		  if ((xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]) > SIZE*SIZE/4.)
+		    {
+		      if (xrand[k] > SIZE/2.)  // look at right image
+			d += (xrand[k] - (s[i].x[k]+SIZE))*
+			  (xrand[k] - (s[i].x[k]+SIZE));
+		      else                     // look at left image
+			d += (xrand[k] - (s[i].x[k]-SIZE))*
+			  (xrand[k] - (s[i].x[k]-SIZE));
+		    }
+		  else
+		    d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
+		}
+	    }
+	  if (d <= (radius + s[i].r)*(radius + s[i].r)) // overlapping!
+	    {
+	      keeper = 0;
+	      counter++;
+	      break;
+	    }
+	}
+      
+      if (hardwallBC)
+	{
+	  for (int k=0; k<DIM; k++)    // check if overlapping wall
+	    {     
+	      if ((xrand[k] <= radius) || (SIZE - xrand[k] <= radius)) // touching wall
+		{
+		  keeper = 0;
+		  counter++; 
+		  break;
+		}
+	    }
+	}
+      
+      if (keeper == 1)
+	break;   
+      
+      if (counter >= 1000)
+	{
+	  std::cout << "counter >= 1000" << std::endl;
+	  exit(-1);
+	}
+    }
+
+  // now convert xrand into index vector for cells
+  vector<DIM,int> cell;
+  cell = vector<DIM>::integer(xrand*((double)(ngrids))/SIZE);
+  
+  s[Ncurrent] = sphere(Ncurrent, xrand, cell, gtime, radius, growth_rate, 
+		       mass, species);
+  
+  //first check to see if entry at cell
+  if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
+    cells.get(cell) = Ncurrent;
+
+  else  // if no, add i to right place in binlist
+    {  
+      int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
+      int pointer = iterater;
+      while (iterater != -1)
+	{
+	  pointer = iterater;
+	  iterater = binlist[iterater];
+	}
+      binlist[pointer] = Ncurrent;
+    }
+
+}
+
+
+//==============================================================
+// Assign cells to spheres read in from existing configuration
+//==============================================================
+void box::AssignCells()
+{
+  for (int i=0; i<N; i++)
+    {
+      // now convert x into index vector for cells
+      vector<DIM,int> cell;
+      cell = vector<DIM>::integer(s[i].x*((double)(ngrids))/SIZE);
+      s[i].cell = cell;
+      
+      //first check to see if entry at cell
+      if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield
+	cells.get(cell) = i;
+      
+      else  // if no, add i to right place in binlist
+	{  
+	  int iterater = cells.get(cell); // now iterate through to end and add Ncurrent
+	  int pointer = iterater;
+	  while (iterater != -1)
+	    {
+	      pointer = iterater;
+	      iterater = binlist[iterater];
+	    }
+	  binlist[pointer] = i;
+	}
+    }
+}
+
+	
+//==============================================================
+// Velocity Giver, assigns initial velocities from Max/Boltz dist.
+//==============================================================
+void box::VelocityGiver(double T)
+{
+  for (int i=0; i<N; i++)
+    {
+      for (int k=0; k<DIM; k++)
+	{
+	  if (T==0.)
+	    s[i].v[k] = 0.;
+	  else
+	    s[i].v[k] = Velocity(T);
+	}
+    }
+}
+
+
+//==============================================================
+// Velocity, gives a single velocity from Max/Boltz dist.
+//==============================================================
+double box::Velocity(double T)
+{
+  double rand;                       // random number between -0.5 and 0.5
+  double sigmasquared = T;    // Assumes M = mass of sphere = 1
+  double sigma = sqrt(sigmasquared); // variance of Gaussian
+  double stepsize = 1000.;           // stepsize for discretization of integral
+  double vel = 0.0;                  // velocity
+  double dv=sigma/stepsize;
+  double p=0.0;
+  
+  rand = (double)random() / (double)RAND_MAX - 0.5;
+  if(rand < 0) 
+    {
+      rand = -rand;
+      dv = -dv;
+    }
+  
+  while(fabs(p) < rand) // integrate until the integral equals rand
+    {
+      p += dv * 0.39894228 * exp(-vel*vel/(2.*sigmasquared))/sigma;
+      vel += dv;
+    }
+  return vel;
+}
+
+
+//==============================================================
+// Finds next events for all spheres..do this once at beginning
+//==============================================================
+void box::SetInitialEvents()
+{
+  for (int i=0; i<N; i++)  // set all events to checks
+    {
+      event e(gtime, i, INF); 
+      s[i].nextevent = e;
+      h.insert(i);
+    }
+}
+
+
+//==============================================================
+// Finds next event for sphere i 
+//==============================================================
+event box::FindNextEvent(int i)
+{
+  double outgrowtime;
+  outgrowtime = (.5/ngrids - (s[i].r+s[i].gr*gtime))/s[i].gr + gtime;
+
+  event t = FindNextTransfer(i);
+  event c = FindNextCollision(i);
+
+  if ((outgrowtime < c.time)&&(outgrowtime < t.time)&&(ngrids>1))
+    {
+      event o = event(outgrowtime,i,INF-1);
+      return o;
+    }
+
+  if ((c.time < t.time)&&(c.j == INF)) // next event is check at DBL infinity
+    return c;
+  else if (c.time < t.time) // next event is collision!
+    {
+      CollisionChecker(c); 
+      return c; 
+    }
+  else // next event is transfer!
+    return t;  
+} 
+
+
+//==============================================================
+// Checks events of predicted collision partner to keep collisions
+// symmetric
+//==============================================================
+void box::CollisionChecker(event c)
+{
+  int i = c.i;
+  int j = c.j;
+  event cj(c.time,j,i,c.v*(-1));
+
+  // j should have NO event before collision with i!
+  if (!(c.time  < s[j].nextevent.time))
+    std::cout << i << " " <<  j << " error collchecker, s[j].nextevent.time= " << s[j].nextevent.time << " " << s[j].nextevent.j << ", c.time= " << c.time << std::endl;
+  
+  int k = s[j].nextevent.j; 
+  if ((k < N) && (k!=i)) // j's next event was collision so give k a check
+    s[k].nextevent.j = INF;
+  
+  // give collision cj to j
+  s[j].nextevent = cj;
+  h.upheap(h.index[j]);
+}
+
+
+//==============================================================
+// Find next transfer for sphere i 
+//==============================================================
+event box::FindNextTransfer(int i)
+{
+  double ttime = dblINF;  
+  int wallindex = INF;   // -(k+1) left wall, (k+1) right wall
+
+  vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
+  vector<DIM> vi = s[i].v;
+
+  for (int k=0; k<DIM; k++)
+    {
+      double newtime;
+      if (vi[k]==0.) 
+	newtime= dblINF;
+      else if (vi[k]>0)  // will hit right wall, need to check if last wall
+	{
+	  if ((hardwallBC)&&(s[i].cell[k] == ngrids - 1))
+	    newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
+		       - (xi[k]+s[i].r+s[i].gr*gtime))/(vi[k]+s[i].gr);
+	  else
+	    newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids))
+		       - xi[k])/(vi[k]);
+	  
+	  if (newtime<ttime)
+	    {
+	      wallindex = k+1;
+	      ttime = newtime;
+	    }
+	}
+      else if (vi[k]<0)  // will hit left wall
+	{
+	  if ((hardwallBC)&&(s[i].cell[k] == 0))
+	    newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) 
+		       - (xi[k]-(s[i].r+s[i].gr*gtime)))/(vi[k]-s[i].gr);
+	  else
+	    newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids)) 
+		       - xi[k])/(vi[k]);
+	  
+	  if (newtime<ttime)
+	    {
+	      wallindex = -(k+1);
+	      ttime = newtime;
+	    }
+	}
+    }
+  
+
+  if (ttime < 0)
+    ttime = 0;
+  // make the event and return it
+  event e = event(ttime+gtime,i,wallindex+DIM+N+1);
+  return e;
+}
+
+
+//==============================================================
+// Check all nearest neighbor cells for collision partners
+//==============================================================
+void box::ForAllNeighbors(int i, vector<DIM,int> vl, vector<DIM,int> vr,
+			  neighbor& operation)
+{
+  vector<DIM,int> cell = s[i].cell;
+
+  // now iterate through nearest neighbors
+  vector<DIM, int> offset;          // nonnegative neighbor offset
+  vector<DIM, int> pboffset;        // nearest image offset
+
+   vector<DIM,int> grid;
+
+   int ii ;
+
+   grid=vl;
+   while(1)
+   {
+     //if (vr[0] > 1)
+     //std::cout << grid << "..." << cell+grid << "\n";
+     for(int k=0; k<DIM; k++)
+     {
+        offset[k]=grid[k]+ngrids;  // do this so no negatives 	
+        if (cell[k]+grid[k]<0) //out of bounds to left
+          pboffset[k] = -1;
+	else if (cell[k]+grid[k]>=ngrids) // out of bounds to right
+	  pboffset[k] = 1;
+        else
+          pboffset[k] = 0;
+     }     
+     int j = cells.get((cell+offset)%ngrids);
+     while(j!=-1)
+       {
+	 operation.Operation(j,pboffset);
+	 j = binlist[j];
+       }
+
+     // A. Donev:     
+     // This code makes this loop dimension-independent
+     // It is basically a flattened-out loop nest of depth DIM
+     for(ii=0;ii<DIM;ii++)
+     {
+       grid[ii] += 1;
+       if(grid[ii]<=vr[ii]) break;
+       grid[ii]=vl[ii];
+     }
+     if(ii>=DIM) break;
+   }  
+}
+
+
+//==============================================================
+// PredictCollision
+//==============================================================
+void box::PredictCollision(int i, int j, vector<DIM, int> pboffset, 
+			     double& ctime, int& cpartner, 
+			     vector<DIM, int>& cpartnerpboffset)
+{
+  double ctimej;
+  
+  if (i!=j)
+    {	 
+      ctimej = CalculateCollision(i,j,pboffset.Double())+gtime;
+      
+      if (ctimej < gtime)
+	std::cout << "error in find collision ctimej < 0" << std::endl;
+      
+      if ((ctimej < ctime)&&(ctimej < s[j].nextevent.time))
+	{
+	  ctime = ctimej;
+	  cpartner = j;
+	  cpartnerpboffset = pboffset;
+	}	
+    }
+}
+
+
+//==============================================================
+// Find next collision
+//==============================================================
+event box::FindNextCollision(int i)
+{
+  collision cc(i, this);
+  
+  vector<DIM, int> vl, vr;
+
+  for (int k=0; k<DIM; k++)  // check all nearest neighbors
+    {
+      vl[k] = -1;
+      vr[k] = 1;
+    }
+  
+  ForAllNeighbors(i,vl,vr,cc);
+
+  event e;
+  if (cc.cpartner == i)  // found no collisions in neighboring cells
+    {
+      if (cc.ctime != dblINF)
+	std::cout << "ctime != dblINF" << std::endl;
+      e = event(dblINF,i,INF);  // give check at double INF
+    }
+  else
+    e = event(cc.ctime,i,cc.cpartner,cc.cpartnerpboffset);
+
+  return e;
+}
+
+
+//==============================================================
+// Calculates collision time between i and image of j using quadratic formula
+//==============================================================
+double box::CalculateCollision(int i, int j, vector<DIM> pboffset)
+{
+  if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12))
+    {
+      return dblINF;
+    }
+  else
+    {      
+      // calculate updated position and velocity of i and j
+      vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
+      vector<DIM> vi = s[i].v;
+      vector<DIM> xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime);
+      vector<DIM> vj = s[j].v;
+      
+      double r_now_i = s[i].r + gtime*s[i].gr;
+      double r_now_j = s[j].r + gtime*s[j].gr;
+      
+      double A,B,C;
+      A = vector<DIM>::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr);
+      B = vector<DIM>::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr);
+      C = vector<DIM>::norm_squared(xi - xj) - (r_now_i+r_now_j)*(r_now_i+r_now_j);
+
+      if (C < -1E-12*(r_now_i+r_now_j))
+	{
+	  std::cout << "error, " << i << " and " << j << 
+	    " are overlapping at time "<< gtime << std::cout;
+	  std::cout << "A, B, C = "  << A << " " << " " << B << 
+	    " " << " " << C <<  std::endl;
+	  if (CheckSphereDiameters()>0)
+	    std::cout << "a sphere has grown greater than unit cell" << 
+	      std::endl;
+	  else
+	    std::cout << "unknown error" << std::cout;
+	  exit(-1);
+	}
+      
+      return QuadraticFormula(A, B, C);
+    }  
+}
+
+
+//==============================================================
+// Quadratic Formula ax^2 + bx + c = 0
+//==============================================================
+ double box::QuadraticFormula(double a, double b, double c)
+{
+  double x = dblINF;
+  double xpos;
+  double xneg;
+  double det = b*b - a*c;
+
+  if (c <= 0.)
+    {
+      if(b < 0.) // spheres already overlapping and approaching
+	{
+	  //std::cout << "spheres overlapping and approaching" << std::endl;
+	  //std::cout << "# events= " << neventstot << std::endl;
+	  x = 0.;	
+	}
+    }
+  else if (det > -10.*DBL_EPSILON)
+    {
+      if (det < 0.)  // determinant can be very small for double roots
+	det = 0.;    
+      if (b < 0.)
+	x = c/(-b + sqrt(det));
+      else if ((a < 0.)&&(b > 0.))
+	x = -(b + sqrt(det))/a;
+      else
+	x = dblINF;
+    }
+  return x;
+}
+
+
+//==============================================================
+// Returns first event
+//==============================================================
+void box::ProcessEvent()
+{  
+  neventstot++;
+  // Extract first event from heap
+  int i = h.extractmax();   
+  event e = s[i].nextevent; // current event
+  event f;                  // replacement event
+
+  if ((e.j>=0)&&(e.j<N))  // collision!
+    {
+      ncollisions++;
+      //std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
+      Collision(e);
+      f = FindNextEvent(i);
+      s[i].nextevent = f;
+      h.downheap(1);
+      if (f.time < e.time)
+	{
+	  std::cout << "error, replacing event with < time" << std::endl;
+	  exit(-1);
+	}
+      
+      // make sure collision was symmetric and give j a check
+      if ((s[e.j].nextevent.j != i)||(s[e.j].nextevent.time != gtime))
+	{
+	  std::cout << "error collisions not symmetric" << std::endl;
+	  std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
+	  std::cout << "but " << e.j << " thinks it has " << s[e.j].nextevent.j<< " "  << s[e.j].nextevent.time << std::endl;
+	  exit(-1);
+	}
+      else  // give j a check
+	s[e.j].nextevent.j = INF;
+    }
+  else if (e.j==INF)      // check!  
+    {
+      nchecks++;
+      //std::cout << "check for " << e.i << " at time " << e.time << std::endl;
+      f = FindNextEvent(i);
+      s[i].nextevent = f;
+      h.downheap(1);
+    }
+  else if (e.j==INF-1)      // sphere outgrowing unit cell, decrease ngrids!
+    {
+      gtime = e.time;
+      Synchronize(false);
+      ngrids = ngrids - 1;
+//      std::cout << "need to reduce ngrids to " << ngrids << std::endl;   
+      ChangeNgrids(ngrids);
+      h.downheap(1);
+    }
+  else                    // transfer!
+    {
+      ntransfers++;
+      //std::cout << "transfer for " << e.i << " at time " << e.time << std::endl;
+      Transfer(e);
+      f = FindNextEvent(i);
+      s[i].nextevent = f;
+      h.downheap(1);
+      //r = FindNextEvent(i, e.j-N-DIM-1);
+      //if (f.time <= e.time)
+      if (f.time < e.time)
+	{
+	  std::cout << "error after transfer, replacing new event with < time" << " " << std::endl;
+	  std::cout << std::setprecision(16) << "e.time= " << e.time << ", f.time= " << f.time << ", f.i= " << f.i << ", f.j= " << f.j << "e.i= " << e.i << ", e.j= " << e.j << std::endl;
+	  std::cout << std::setprecision(16) << "difference= " << e.time - f.time << std::endl;
+	  exit(-1);
+	}
+    }
+}
+
+
+//==============================================================
+// Processes a collision
+//=============================================================
+void box::Collision(event e)
+{
+  double ctime = e.time;
+  int i = e.i;
+  int j = e.j;
+  vector<DIM,int> v = e.v;  // virtual image
+  gtime = ctime;
+
+  // Update positions and cells of i and j to ctime
+  s[i].x += s[i].v*(gtime-s[i].lutime);
+  s[j].x += s[j].v*(gtime-s[j].lutime);
+
+  // Check to see if a diameter apart
+  double r_sum = s[i].r + s[j].r + (s[i].gr+s[j].gr)*gtime;
+  double distance = vector<DIM>::norm_squared(s[i].x - s[j].x- v.Double()*SIZE) - r_sum*r_sum;
+  if (distance*distance > 10.*DBL_EPSILON)
+    std::cout << "overlap " << distance << std::endl;
+
+  s[i].lutime = gtime;
+  s[j].lutime = gtime;
+  
+  vector<DIM,double> vipar;          // parallel comp. vi
+  vector<DIM,double> vjpar;          // parallel comp. vj
+  vector<DIM,double> viperp;         // perpendicular comp. vi
+  vector<DIM,double> vjperp;         // perpendicular comp. vj
+
+  // make unit vector out of displacement vector
+  vector<DIM,double> dhat;
+  dhat = s[i].x - s[j].x - v.Double()*SIZE;  // using image of j!!
+  double dhatmagnitude = sqrt(dhat.norm_squared());
+  dhat /= dhatmagnitude;
+
+  vipar = dhat*vector<DIM>::dot(s[i].v, dhat);
+  vjpar = dhat*vector<DIM>::dot(s[j].v, dhat);
+  viperp = s[i].v - vipar;
+  vjperp = s[j].v - vjpar;
+
+  s[i].v = vjpar + dhat*(s[i].gr+s[j].gr)*2 + viperp;
+  s[j].v = vipar - dhat*(s[i].gr+s[j].gr)*2 + vjperp;
+
+  // momentum exchange
+  double xvelocity;   // exchanged velocity
+  xvelocity = vector<DIM>::dot(s[i].v - s[j].v, dhat) - (s[i].gr+s[j].gr);
+  xmomentum += xvelocity*dhatmagnitude*s[i].m*s[j].m*2/(s[i].m+s[j].m);
+}
+
+
+//==============================================================
+// Transfer, takes care of boundary events too
+//=============================================================
+void box::Transfer(event e)
+{
+  gtime = e.time;
+  int i = e.i;
+  int j = e.j;
+  int k=0;           // dimension perpendicular to wall it crosses
+ 
+  // update position and lutime (velocity doesn't change)
+  s[i].x += s[i].v*(gtime-s[i].lutime);
+  s[i].lutime = gtime;
+
+  vector<DIM,int> celli;  // new cell for i
+  celli = s[i].cell;  // this is not redundant
+  
+  // update cell
+  if (j>N+DIM+1)  // right wall
+    {
+      k = j-N-DIM-2;
+      celli[k] = s[i].cell[k] + 1;
+
+      if (hardwallBC)
+	{
+	  // if in right-most cell, reflect back
+	  if (s[i].cell[k] == ngrids - 1)
+	    s[i].v[k] = -s[i].v[k] - 4.*s[i].gr;
+	  else
+	    {
+	      if (ngrids>1)
+		UpdateCell(i, celli); 
+	    }
+	}
+      else
+	{
+	  // if in right-most cell, translate x and cell
+	  if (s[i].cell[k] == ngrids - 1)
+	    {
+	      s[i].x[k] -= SIZE;
+	      celli[k] -= ngrids;
+	    }
+	  if (ngrids>1)
+	    UpdateCell(i, celli); 
+	}
+    }
+  else if (j<N+DIM+1)  // left wall
+    {
+      k = -j+N+DIM;
+      celli[k] = s[i].cell[k] - 1;
+
+      if (hardwallBC)
+	{
+	  // if in left-most cell, reflect back
+	  if (s[i].cell[k] == 0)
+	    s[i].v[k] = -s[i].v[k] + 4.*s[i].gr;
+	  else
+	    {
+	      if (ngrids>1)
+		UpdateCell(i, celli); 
+	    }
+	}
+      else
+	{
+	  // if in left-most cell, translate x and cell
+	  if (s[i].cell[k] == 0)
+	    {
+	      s[i].x[k] += SIZE;
+	      celli[k] += ngrids;
+	    }
+	  if (ngrids>1)
+	    UpdateCell(i, celli); 
+	}
+    }
+  else
+    std::cout << "error in Transfer" << std::endl;    
+  
+}
+
+
+//==============================================================
+// Updates cell of a sphere to time
+//=============================================================
+void box::UpdateCell(int i, vector<DIM,int>& celli)
+{
+  if (celli == s[i].cell)
+    std::cout << "error in update cell..shouldn't be the same" << std::endl;
+  
+  // delete i from cell array at cell
+
+  if (cells.get(s[i].cell) == i) 
+    {
+      if (binlist[i] == -1)
+	cells.get(s[i].cell) = -1;
+      else
+	{
+	  cells.get(s[i].cell) = binlist[i];
+	  binlist[i] = -1;
+	}
+    }
+
+  else if (cells.get(s[i].cell) == -1)
+    {
+      std::cout << "error " << i << " not in claimed cell UpdateCell" << std::endl;
+      OutputCells();
+    }
+
+  else  // if no, find i in binlist
+    {  
+      int iterater = cells.get(s[i].cell);
+      int pointer = iterater;
+      while ((iterater != i)&&(iterater != -1))
+	{
+	  pointer = iterater;
+	  iterater = binlist[iterater];
+	}
+      if (iterater == -1)  // got to end of list without finding i
+	{
+	  std::cout << "problem " << i << " wasn't in claimed, cell iterater = -1" << std::endl;
+	  OutputCells();
+	}
+      else  // we found i!
+	{
+	  binlist[pointer] = binlist[i]; 
+	  binlist[i] = -1;
+	}	  
+    } 
+
+  // now add i to cell array at celli
+  s[i].cell = celli;
+  
+  //first check to see if entry at celli
+  if (cells.get(celli) == -1) //if yes, add i to cells gridfield
+    cells.get(celli) = i;
+  else  // if no, add i to right place in binlist
+    {
+      int iterater = cells.get(celli);  // now iterate through to end and add i
+      int pointer = iterater;
+      while (iterater != -1)  // find the end of the list
+	{
+	  pointer = iterater;
+	  iterater = binlist[iterater];
+	}
+      binlist[pointer] = i;
+      binlist[i] = -1; // redundant
+    }
+}
+
+
+//==============================================================
+// Output event heap...purely used for debugging
+//==============================================================
+void box::OutputEvents()
+{
+  h.print();
+}
+
+
+//==============================================================
+// Output positions of spheres and their cells...purely used for debugging
+//==============================================================
+void box::OutputCells()
+{
+  for (int i=0; i<N; i++)
+    std::cout << i << " " << s[i].x << " " << s[i].v << " " << s[i].cell << std::endl;
+}
+
+
+//==============================================================
+// Update positions...purely for graphical display
+//==============================================================
+void box::TrackPositions()
+{
+  for (int i=0; i<N; i++)
+    x[i] = s[i].x + s[i].v*(gtime-s[i].lutime);
+}
+
+
+//==============================================================
+// Computes the total energy
+//==============================================================
+double box::Energy()
+{
+  double E=0;
+  for (int i=0; i<N; i++)
+    E += 0.5*s[i].m*s[i].v.norm_squared();
+
+  return E/N;
+}
+
+
+//==============================================================
+// Calculates the packing fraction
+//==============================================================
+double box::PackingFraction()
+{
+  double rfactor = 0.;
+  for (int i=0; i<N; i++)
+    {
+      rfactor += pow(s[i].r + gtime*s[i].gr, DIM);
+    }
+  double v = (rfactor*pow(sqrt(PI), DIM))/(exp(lgamma(1.+((double)(DIM))/2.)));
+  return v/(pow(SIZE, DIM));
+}
+
+
+//==============================================================
+// Checks to make sure all sphere diameters are less than dimension
+// of unit cell
+//==============================================================
+int box::CheckSphereDiameters()
+{
+  int offender = 0;
+  for (int i=0; i<N; i++)
+    {
+      if (s[i].r*2 > 1./ngrids){
+	offender = i;
+	break;
+      }
+    }
+  return offender;
+}
+
+
+//==============================================================
+// Change ngrids
+//==============================================================
+void box::ChangeNgrids(int newngrids)
+{
+  cells.set_size(newngrids);
+  cells.initialize(-1);      // initialize cells to -1
+  for (int i=0; i<N; i++)    // initialize binlist to -1
+    binlist[i] = -1;
+  AssignCells();
+  for (int i=0; i<N; i++)
+    s[i].nextevent = event(0., i, INF); 
+  Process(N); 
+}	 
+
+
+//==============================================================
+// Calculates the optimal ngrids for the initial configuration
+// and assumes that ngrids gets updated (reduced) as the packing
+// proceeds
+//==============================================================
+int box::Optimalngrids2(double currentradius)
+{
+  return (int)(1./(2.*currentradius));
+}
+
+
+//==============================================================
+// Calculates the optimal ngrids, assuming ngrids is not updated
+// automatically and is very conservative
+//==============================================================
+int box::Optimalngrids()
+{
+  double maxr;
+
+  maxr = pow(exp(lgamma(1.+((double)(DIM))/2.))*maxpf/
+	     (N*(bidispersityfraction + (1.-bidispersityfraction)*
+		 pow(bidispersityratio, DIM))), 
+	     1./DIM)/sqrt(PI);
+
+  return (int)(1./(2.*maxr));
+}
+
+
+//==============================================================
+// Processes n events
+//==============================================================
+void box::Process(int n)
+{
+  double deltat = gtime;
+  for (int i=0; i<n; i++)
+    {
+      ProcessEvent();
+    }
+  pf = PackingFraction();   // packing fraction
+  deltat = gtime - deltat;
+  double oldenergy = energy;
+  energy = Energy();        // kinetic energy
+
+  energychange = ((oldenergy - energy)/oldenergy)*100; // percent change in energy
+
+  if (deltat != 0.)
+    {
+      pressure = 1+xmomentum/(2.*energy*N*deltat);
+      collisionrate = ((double)(ncollisions))/deltat;
+    }  
+
+  // reset to 0
+  ncollisions = 0;
+  ntransfers = 0;
+  nchecks = 0;
+  xmomentum = 0.;
+  ncycles++;
+}
+
+
+//==============================================================
+// Prints statistics for n events
+//==============================================================
+void box::PrintStatistics()
+{
+  std::cout << "packing fraction = " << pf << std::endl;
+  std::cout << "gtime = " << gtime << std::endl; 
+  std::cout << "total time = " << rtime+gtime << std::endl;
+  std::cout << "kinetic energy = " << energy << std::endl;
+  std::cout << "total # events = " << neventstot << std::endl;
+  std::cout << "# events = " << ncollisions+ntransfers+nchecks << ", # collisions = " << ncollisions << ", # transfers = " << ntransfers << ", # checks =" << nchecks << std::endl;
+  std::cout << "growthrate = " << growthrate << std::endl;
+  std::cout << "collisionrate = " << collisionrate << std::endl;
+  std::cout << "reduced pressure = " << pressure << std::endl;
+  std::cout << "-----------------" << std::endl;
+}
+
+
+//==============================================================
+// Updates spheres to gtime, synchronizes, and can change growth rate
+//==============================================================
+void box::Synchronize(bool rescale)
+{
+  double vavg = sqrt(2.*M*energy);
+
+  for (int i=0; i<N; i++)
+    {
+      s[i].x = s[i].x + s[i].v*(gtime-s[i].lutime);
+      s[i].nextevent.time -= gtime;
+
+      if (s[i].nextevent.time < 0.)
+	std::cout << "error, event times negative after synchronization" << std::endl;
+      if (rescale == true)   // give everyone checks
+	{
+	  s[i].nextevent = event(0., i, INF); 
+	  s[i].v /= vavg;
+	}
+	  
+      s[i].lutime = 0.;
+      s[i].r += gtime*s[i].gr;   
+    }
+
+  //r += gtime*growthrate;       // r defined at gtime = 0
+  rtime += gtime;
+  gtime = 0.;
+  
+  if (rescale == true)
+    Process(N);
+}
+
+
+//==============================================================
+// Run time
+//==============================================================
+void box::RunTime()
+{
+  time(&end);
+  std::cout << "run time = " << difftime(end, start) << std::endl;
+}
+
+
+//==============================================================
+// Write configuration
+//==============================================================
+void box::WriteConfiguration(double *data)
+{
+  int count1;
+
+  if (gtime != 0.)   // synchronize spheres if not currently synchronized
+    Synchronize(false);
+      
+  
+ for (int i=0; i<N; i++)  // output radius and positions
+   {
+     for (int k=0; k<DIM; k++)
+			 data[DIM * i + k] = s[i].x[k];
+   }
+ 
+}
diff --git a/src/spheres_poly/box.h b/src/spheres_poly/box.h
new file mode 100644
index 0000000..69418f4
--- /dev/null
+++ b/src/spheres_poly/box.h
@@ -0,0 +1,169 @@
+/* 
+   Packing of hard spheres via molecular dynamics
+   Developed by Monica Skoge, 2006, Princeton University
+   Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions
+   This code may be used, modified and distributed freely.
+   Please cite:
+   
+   "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+   	M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+	
+   if you use these codes.	
+*/
+
+//-----------------------------------------------------------------------------
+// Box maker
+//---------------------------------------------------------------------------
+
+#ifndef  BOX_H
+#define  BOX_H
+
+
+#include <vector>
+#include <math.h>
+
+#include "vector.h"
+#include "grid_field.h"
+#include "event.h"
+#include "sphere.h"
+#include "heap.h"
+
+
+#define PI     3.141592653589793238462643
+#define SIZE   1.0            // size of box
+#define VOLUMESPHERE pow(PI,((double)(DIM))/2.)/exp(lgamma(1+((double)(DIM))/2.)) // volume prefactor for sphere
+#define DBL_EPSILON  2.2204460492503131e-016 // smallest # such that 1.0+DBL_EPSILON!=1.0
+#define M 1.0
+
+//---------------------------------------------------------------------------
+// Class neighbor
+//---------------------------------------------------------------------------
+class neighbor
+{
+ public:
+  int i;
+  
+  neighbor(int i_i);
+
+ public:
+  virtual void Operation(int j, vector<DIM, int>& pboffset) = 0;
+};
+
+
+class box {
+  
+ public:
+  
+  // constructor and destructor
+  box(int N_i, double r_i, double growthrate_i, double maxpf_i, 
+      double bidispersityratio, double bidispersityfraction, 
+      double massratio, int hardwallBC);
+  ~box();
+
+  // Creating configurations
+  int Optimalngrids();
+  int Optimalngrids2(double maxr);
+  void CreateSpheres(double temp);
+  void CreateSphere(int Ncurrent);   
+  double Velocity(double temp);
+  void VelocityGiver(double temp);  
+  void SetInitialEvents();
+  void RecreateSpheres(const char* filename, double temp);
+  void ReadPositions(const char* filename);
+  void AssignCells();
+ 
+  // Predicting next event
+  event FindNextEvent(int i);
+  void CollisionChecker(event c);
+  event FindNextTransfer(int i);
+  event FindNextCollision(int i);
+  void ForAllNeighbors(int, vector<DIM, int>, vector<DIM,int>, neighbor&);
+  void PredictCollision(int i, int j, vector<DIM, int> pboffset, 
+			double& ctime, int& cpartner, 
+			vector<DIM, int>& cpartnerpboffset);
+  double CalculateCollision(int i, int j, vector<DIM>  pboffset);
+  double QuadraticFormula(double a, double b, double c);
+  
+  // Processing an event
+  void Process(int n);
+  void ProcessEvent();
+  void Collision(event e);
+  void Transfer(event e);
+  void UpdateCell(int i, vector<DIM,int>& celli);
+  void Synchronize(bool rescale);
+  void ChangeNgrids(int newngrids);
+
+  // Debugging
+  void TrackPositions();
+  void OutputEvents();
+  void OutputCells();
+  void GetInfo();  
+  int CheckSphereDiameters();
+
+  // Statistics
+  double Energy();
+  double PackingFraction();  
+  void PrintStatistics();
+  void RunTime();
+  void WriteConfiguration(double* data);
+  
+  
+  //variables
+
+  const int N;                   // number of spheres
+  
+  int ngrids;                    // number of cells in one direction
+  double maxpf;
+  double growthrate;             // growth rate of the spheres
+  double r;                      // radius, defined at gtime = 0
+  double gtime;                  // this is global clock
+  double rtime;                  // reset time, total time = rtime + gtime
+  double collisionrate;          // average rate of collision between spheres
+  double bidispersityratio;      // ratio of sphere radii
+  double bidispersityfraction;   // fraction of smaller spheres
+  double massratio;              // ratio of sphere masses
+  int hardwallBC;                // =0 for periodic BC, =1 for hard wall
+  
+
+  // statistics
+  double pressure;               // pressure
+  double xmomentum;              // exchanged momentum
+  double pf;                     // packing fraction
+  double energy;                 // kinetic energy
+  double energychange;
+  int ncollisions;               // number of collisions
+  int ntransfers;                // number of transfers
+  int nchecks;                   // number of checks
+  int neventstot;                // total number of events 
+  int ncycles;                   // counts # cycles for output
+
+  time_t start, error, end;      // run time of program
+
+  // arrays
+  sphere *s;                      // array of spheres
+  grid_field<DIM, int> cells; // array that keeps track of spheres in each cell
+  int *binlist;                   // linked-list for cells array
+  heap h;                         // event heap
+  vector<DIM> *x;                 // positions of spheres.used for graphics
+};
+
+
+//---------------------------------------------------------------------------
+// Predicts collisions, inherits neighbor operation
+//---------------------------------------------------------------------------
+class collision : public neighbor 
+{
+ public:
+  
+  box *b; 
+  double ctime;
+  int cpartner;
+  vector<DIM,int> cpartnerpboffset;
+
+ public:
+  collision(int i_i, box *b);
+
+  virtual void Operation(int j, vector<DIM, int>& pboffset);
+};
+
+#endif 
diff --git a/src/spheres_poly/event.cpp b/src/spheres_poly/event.cpp
new file mode 100644
index 0000000..c8c104d
--- /dev/null
+++ b/src/spheres_poly/event.cpp
@@ -0,0 +1,65 @@
+#include "event.h" 
+
+//==============================================================
+//==============================================================
+//  Class Event 
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+event::event(double time_i, int i_i, int j_i, vector<DIM,int> v_i):
+  time(time_i),
+  i(i_i),
+  j(j_i),
+  v(v_i)
+{
+}
+
+event::event(double time_i, int i_i, int j_i):
+  time(time_i),
+  i(i_i),
+  j(j_i)
+{
+}
+
+event::event(const event& e)
+{
+  time = e.time;
+  i = e.i;
+  j = e.j;
+  v = e.v;
+}
+
+event::event()
+{
+}
+
+  
+//==============================================================
+// Destructor
+//==============================================================
+event::~event() 
+{
+}
+
+void event::erase()
+{
+  time = dblINF;
+  i = 0;
+  j = 0;
+}
+
+bool event::operator<(const event& e) const
+{
+  return e.time < time;
+}
+
+bool event::operator>(const event& e) const
+{
+  return e.time > time;
+}
+
+ 
diff --git a/src/spheres_poly/event.h b/src/spheres_poly/event.h
new file mode 100644
index 0000000..55dd1fc
--- /dev/null
+++ b/src/spheres_poly/event.h
@@ -0,0 +1,58 @@
+#ifndef  EVENT_H
+#define  EVENT_H
+
+#include "vector.h"
+#define INF    100000000
+#define dblINF 100000000.
+
+class event {
+
+ public:
+
+  // constructor and destructor
+  event(double time_i, int i_i, int j_i, vector<DIM,int> v_i);
+  event(double time_i, int i_i, int j_i);
+  event(const event& e);
+  event();
+  
+  ~event();
+
+  bool operator<(const event&) const;
+  bool operator>(const event&) const;
+  void erase();
+
+ //variables
+
+
+  double time;             // time of next collision
+  int i;        // collision partner with lower number
+  int j;        // collision partner with higher number
+  vector<DIM,int> v;        // virtual image
+
+  /* 0<=j<=N                binary collision between i and j
+     j=N+DIM+1+x            transfer where x=-(k+1) for left wall
+                            and x=k+1 for right wall
+     j=INF                 both check after event that did not altered motion of                           i and check after event that altered motion of i, i.e                           rescaling of velocities. I currently don't see need t                           o separate the two
+
+     j=-1                  check after collision
+
+     Virtual identifiers as scalars...I think bad idea, but here's my work
+     there will be easily fixed problems if DIM>=10
+     -x<=v<=x               x=k+1, image in k direction
+     v=xy                   x,y
+      =-xy                  -x,y
+      =-yx                  x,-y
+      =yx                   -x,-y
+     v=xyz                  x,y,z
+      =-xyz                 -x,y,z
+      =-yxz                 x,-y,z
+      =-zxy                 x,y,-z
+      =zyx                  -x,-y,z
+      =xzy                  x,-y,-z
+      =yzx                  -x,y,-z
+      =-zyx                 -x,-y,-z
+  */      
+
+};
+
+#endif 
diff --git a/src/spheres_poly/grid_field.h b/src/spheres_poly/grid_field.h
new file mode 100644
index 0000000..10430ed
--- /dev/null
+++ b/src/spheres_poly/grid_field.h
@@ -0,0 +1,148 @@
+/* 
+   Packing of hard spheres via molecular dynamics
+   Developed by Monica Skoge, 2006, Princeton University
+   Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions
+   This code may be used, modified and distributed freely.
+   Please cite:
+   
+   "Packing Hyperspheres in High-Dimensional Euclidean Spaces"
+   	M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006
+	
+   if you use these codes.	
+*/
+
+
+#ifndef GRID_FIELD_H
+#define GRID_FIELD_H
+
+#include "vector.h"
+
+// ======================================================================
+// grid_field 
+// ======================================================================
+
+// A field of V-vectors on a D dimensional manifold
+
+template<int D, class T>
+class grid_field {
+  
+ public:
+  int elements;
+
+ private:
+  T* f;
+  vector<D, int> size;           // number of grid points for each dimension
+  vector<D, int> offset;
+ 
+ public:
+
+  grid_field();
+  grid_field(const vector<D, int>&);
+  ~grid_field();
+
+  T& get(const vector<D, int>&);
+
+  vector<D, int> get_size() const;
+  void set_size(const vector<D, int>&);
+  void set_size(const int);
+  void initialize(const int i);
+};
+
+
+// grid_field
+// ~~~~~~~~~~~~
+template<int D, class T>
+grid_field<D, T>::grid_field()
+  : f(0), elements(0)
+{
+}
+
+
+// grid_field
+// ~~~~~~~~~~~~
+template<int D, class T>
+grid_field<D, T>::grid_field(const vector<D, int>& s)
+  : f(0)
+{
+  set_size(s);
+}
+
+
+// ~grid_field
+// ~~~~~~~~~~~~~
+template <int D, class T>
+grid_field<D, T>::~grid_field()
+{
+  if(f != 0)
+    delete[] f;
+}
+
+
+// get_size
+// ~~~~~~~~
+template<int D, class T>
+inline vector<D, int> grid_field<D, T>::get_size() const
+{
+  return size;
+}
+
+
+// set_size
+// ~~~~~~~~
+template<int D, class T>
+void grid_field<D, T>::set_size(const vector<D, int>& s)
+{
+  if(f != 0)
+    delete[] f;
+
+  size = s;
+
+  elements = 1;
+  for(int i=0; i<D; i++) {
+    offset[i] = elements;
+    elements *= size.x[i];
+  }
+
+  f = new T[elements];
+}
+
+
+// set_size
+// ~~~~~~~~
+template<int D, class T>
+void grid_field<D, T>::set_size(const int s)
+{
+  vector<D, int> square;
+
+  for(int k=0; k<D; k++)
+    square[k] = s;
+
+  set_size(square);
+}
+
+
+// get
+// ~~~
+template<int D, class T>
+inline T& grid_field<D, T>::get(const vector<D, int>& pos)
+{
+  int p=0;
+  for(int i=0; i<D; i++)
+    p += pos.x[i]*offset[i];
+
+  return f[p];
+}
+
+
+// initialize
+// ~~~
+template<int D, class T>
+void grid_field<D, T>::initialize(const int value)
+{
+  for(int i=0; i<elements; i++)
+    f[i] = value;
+}
+
+
+
+#endif
diff --git a/src/spheres_poly/heap.cpp b/src/spheres_poly/heap.cpp
new file mode 100644
index 0000000..c83865c
--- /dev/null
+++ b/src/spheres_poly/heap.cpp
@@ -0,0 +1,186 @@
+#include "heap.h"
+#include "event.h"
+#include <iostream>
+#include "box.h"
+
+
+//==============================================================
+//==============================================================
+//  Class heap: Event heap used in hard spheres calculation
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+heap::heap(int maxsize_i):
+  maxsize(maxsize_i)
+{
+  a = new int[maxsize];
+  index = new int[maxsize];
+  
+  N = 0;   // current number of events in heap
+}
+
+
+//==============================================================
+// Constructor
+//==============================================================
+heap::heap(const heap &h)
+{
+  maxsize = h.maxsize;
+  a = h.a;
+  index = h.index;
+  N = h.N;                // current number of events in heap
+  s = h.s;
+}
+
+
+//==============================================================
+// Destructor
+//==============================================================
+heap::~heap()
+{
+  delete[] a;
+  delete[] index;
+}
+
+//==============================================================
+// Upheap
+//==============================================================
+void heap::upheap(int k)
+{
+  int i = a[k];
+
+  while ((k>1)&&(s[a[k/2]].nextevent.time > s[i].nextevent.time))
+    {
+      a[k] = a[k/2];
+      index[a[k/2]] = k;
+      k = k/2;
+    }
+  a[k] = i;
+  index[i] = k;
+}
+
+//==============================================================
+// Downheap
+//==============================================================
+void heap::downheap(int k)
+{
+  int j;
+  int i = a[k];
+  
+  while(k <= N/2)
+    {
+      j = k+k;
+      if ((j < N)&&(s[a[j]].nextevent.time > s[a[j+1]].nextevent.time))
+	j++;
+      if (s[i].nextevent.time <= s[a[j]].nextevent.time)
+	break;
+      a[k] = a[j];
+      index[a[j]] = k;
+      k = j;
+    }
+  a[k] = i;
+  index[i] = k;
+}
+
+//==============================================================
+// Insert
+//==============================================================
+void heap::insert(int i)
+{
+  if (N >= maxsize)
+    std::cout << "error, N >= maxsize, cannot insert another event" << std::endl;
+  else
+    {
+      N++;
+      a[N] = i;
+      index[i] = N;
+      upheap(N);
+    }  
+}
+
+
+//==============================================================
+// Extract max
+//==============================================================
+int heap::extractmax()
+{
+  return a[1];
+}
+
+/*
+//==============================================================
+// Replace
+//==============================================================
+void heap::replace(int i)
+{
+  a[1] = i;
+  s[i].nextevent = t;
+
+  if (!(e.time > s[a[1]].nextevent.time))
+    {
+      if (!(s[a[1]].nextevent.j == INF))// must be check i'm changing to coll. at same time
+	std::cout << "error replaced non check with earlier time" << std::endl;
+      a[1] = i;
+      index[i] = 1;
+    }
+  else
+    {    
+      a[0] = i;
+      downheap(0);
+    }
+}
+*/
+
+//==============================================================
+// Search
+//==============================================================
+int heap::search(int j)
+{
+  for (int k=1; k<=N; k++)
+    {
+      if (a[k] == j)
+	return k;
+    }
+  return -1;
+}
+
+/*
+//==============================================================
+// Change
+//==============================================================
+void heap::change(int i)
+{
+  int iindex = index[i];
+  
+  if (s[i].nextevent.time == s[a[iindex]].nextevent.time)  
+    std::cout << "error changing an event to an equal time" << std::endl;
+  else if (s[i].nextevent.time > s[a[iindex]].nextevent.time)
+    std::cout << "error changing an event to a greater time" << std::endl;
+
+  a[iindex] = i;
+  upheap(iindex);
+}
+*/
+
+//==============================================================
+// Print
+//==============================================================
+void heap::print()
+{
+  for (int k=1; k<=N; k++)
+    std::cout << k << " " << a[k] << " " << s[a[k]].nextevent.j << " " << s[a[k]].nextevent.time << std::endl;
+}
+
+
+//==============================================================
+// Check index array
+//==============================================================
+void heap::checkindex()
+{
+  for (int k=1; k<=N; k++)
+    std::cout << k << " " << a[k] << " " << index[a[k]] << std::endl;
+}
diff --git a/src/spheres_poly/heap.h b/src/spheres_poly/heap.h
new file mode 100644
index 0000000..91180e3
--- /dev/null
+++ b/src/spheres_poly/heap.h
@@ -0,0 +1,42 @@
+//---------------------------------------------------------------------------
+// Event heap maker
+//---------------------------------------------------------------------------
+
+#ifndef  HEAP_H
+#define  HEAP_H
+
+#include "event.h"
+#include "sphere.h"
+
+class heap {
+
+ public:
+
+  // constructor and destructor
+  heap(int maxsize);
+  heap(const heap &h);
+  ~heap();
+
+  // variables
+  int maxsize;   // max allowed number of events
+  int N;         // current number of events
+  int *a;
+  sphere *s;
+  int *index;     // array of indices for each sphere
+  //event minevent;
+
+
+  // functions which operate on a binary heap
+  
+  void upheap(int k);
+  void downheap(int k);
+  void insert(int i);
+  void replace(int i);
+  int search(int j);
+  void change(int i); 
+  int extractmax();
+  void print();
+  void checkindex();
+  
+};
+#endif 
diff --git a/src/spheres_poly/neighbor.cpp b/src/spheres_poly/neighbor.cpp
new file mode 100644
index 0000000..920c099
--- /dev/null
+++ b/src/spheres_poly/neighbor.cpp
@@ -0,0 +1,40 @@
+#include "vector.h"
+#include <iostream>
+#include "box.h"
+
+
+//==============================================================
+//==============================================================
+//  Class neighbor
+//==============================================================
+//==============================================================
+
+neighbor::neighbor(int i_i)
+ : i(i_i)
+{
+}
+
+
+//==============================================================
+//==============================================================
+//  Class collision
+//==============================================================
+//==============================================================
+
+collision::collision(int i_i, box *b_i)
+  : neighbor(i_i),
+    b(b_i)
+{
+  ctime = dblINF;
+  cpartner = i;
+}
+
+
+//==============================================================
+// Operation is finding the next collision from a given cell
+//==============================================================
+void collision::Operation(int j, vector<DIM, int>& pboffset)
+{
+  b->PredictCollision(i, j, pboffset, ctime, cpartner, cpartnerpboffset);
+}
+
diff --git a/src/spheres_poly/sphere.cpp b/src/spheres_poly/sphere.cpp
new file mode 100644
index 0000000..0fbbd68
--- /dev/null
+++ b/src/spheres_poly/sphere.cpp
@@ -0,0 +1,66 @@
+#include "box.h"
+#include <iostream>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+
+#include "vector.h"
+
+//==============================================================
+//==============================================================
+//  Class Sphere: 
+//==============================================================
+//==============================================================
+
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere()
+{
+}
+
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere(const sphere& s)
+{
+  i = s.i;
+  x = s.x;
+  v = s.v;
+  cell = s.cell;
+  lutime = s.lutime;
+  nextevent = s.nextevent;
+  nextcollision = s.nextcollision;
+  r = s.r;
+  gr = s.gr;
+  m = s.m;
+  species=s.species;
+}
+
+//==============================================================
+// Constructor
+//==============================================================
+sphere::sphere(int i_i, vector<DIM> x_i, vector<DIM, int> cell_i, 
+	       double lutime_i, double r_i, double gr_i, double m_i, int species_i):
+  i(i_i),
+  x(x_i),
+  cell(cell_i),
+  lutime(lutime_i),
+  r(r_i),
+  gr(gr_i),
+  m(m_i),
+  species(species_i)
+{
+}
+
+//==============================================================
+// Destructor
+//==============================================================
+sphere::~sphere() 
+{
+
+}
+
+ 
diff --git a/src/spheres_poly/sphere.h b/src/spheres_poly/sphere.h
new file mode 100644
index 0000000..28c64b4
--- /dev/null
+++ b/src/spheres_poly/sphere.h
@@ -0,0 +1,44 @@
+#ifndef  SPHERE_H
+#define  SPHERE_H
+
+
+#include "vector.h"
+
+
+class sphere {
+
+ public:
+
+  // constructor and destructor
+
+  sphere();
+  sphere(const sphere& s);
+  sphere(int i_i, vector<DIM> x, vector<DIM, int> cell_i, double lutime_i, 
+	 double r_i, double gr_i, double m_i, int species_i);
+  ~sphere();
+
+ //variables
+  
+  int i;                          // sphere ID
+
+  // impending event
+  event nextevent;                // next event...can be collision or transfer
+  event nextcollision;            // next collision if next event is transfer
+  // maybe nextnext event
+  
+  // past information
+  double lutime;                  // last update time
+  vector<DIM, int> cell;          // cell that it belongs to
+  vector<DIM, double> x;          // position
+  vector<DIM, double> v;          // velocity
+  double r;                       // sphere radius
+  double gr;                      // sphere growth rate
+  double m;                       // sphere mass
+  int species;                    // species number (not used during the MD)
+  // make sure efficent in memory
+
+ 
+
+};
+
+#endif 
diff --git a/src/spheres_poly/spheres.cpp b/src/spheres_poly/spheres.cpp
new file mode 100644
index 0000000..9b14b01
--- /dev/null
+++ b/src/spheres_poly/spheres.cpp
@@ -0,0 +1,64 @@
+//===========================================================
+//===========================================================
+//===========================================================
+//
+//  Molecular dynamics simulation of hardspheres
+//
+//===========================================================
+//===========================================================
+//===========================================================
+
+#include <iostream>
+#include <math.h>
+#include <fstream>
+#include <vector>
+#include <time.h>
+#include <string.h>
+#include <stdlib.h>
+#include <cstdlib>
+
+#include "box.h"
+#include "sphere.h"
+#include "event.h"
+#include "heap.h"
+
+
+double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate)
+{
+	int eventspercycle = 20;              // # events per sphere per cycle
+	double initialpf = 0.01;              // initial packing fraction
+	double temp = 0.2;                    // initial temp., use 0 for zero velocities
+	double growthrate = 0.001;            // growth rate
+	double massratio = 1.;               // ratio of sphere masses
+	int hardwallBC = 0;                   // 0 for periodic, 1 for hard wall BC
+
+  double d, r;   // initial diameter and radius of spheres
+
+  r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM)));
+
+  box b(N, r, growthrate, maxpf, bidispersityratio, 
+	bidispersityfraction, massratio, hardwallBC);
+
+	b.CreateSpheres(temp);
+  
+  
+  while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure)) 
+    {
+      b.Process(eventspercycle*N);
+      b.Synchronize(true);
+    }
+  
+	double *data = (double *)malloc(DIM * N * sizeof(double));
+
+  b.WriteConfiguration(data);
+  
+  return data;
+}
+
+extern "C" {
+	double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) {
+		return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate);
+	}
+}
+
+
diff --git a/src/spheres_poly/spheres.h b/src/spheres_poly/spheres.h
new file mode 100644
index 0000000..44d8052
--- /dev/null
+++ b/src/spheres_poly/spheres.h
@@ -0,0 +1,4 @@
+
+#pragma once
+
+double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate);
diff --git a/src/spheres_poly/vector.h b/src/spheres_poly/vector.h
new file mode 100644
index 0000000..9ee481f
--- /dev/null
+++ b/src/spheres_poly/vector.h
@@ -0,0 +1,471 @@
+#ifndef VECTOR_H
+#define VECTOR_H
+
+#include <iostream>
+#include <fstream>
+
+#define DIM 2
+
+template <int D, typename T=double>
+class vector {
+  
+  public:
+  T x[D];
+  
+  public:
+  
+  vector();
+  vector(const T[D]);           
+  vector(const vector&);
+  ~vector();
+
+  vector<D, T>& operator+=(const vector<D, T>&);
+  vector<D, T>& operator-=(const vector<D, T>&);
+  vector<D, T>& operator*=(const T);
+  vector<D, T>& operator/=(const T);
+  vector<D, T> operator+(const vector<D, T>&) const;
+  vector<D, T> operator-(const vector<D, T>&) const;
+  vector<D, T> operator*(const T) const;
+  vector<D, T> operator/(const T) const;
+  vector<D, T> operator%(const T) const;
+  bool operator==(const vector<D, T> &a) const;
+
+  vector<D, int> integer() const;
+  vector<D, double> Double() const;
+  static vector<D, int> integer(const vector<D, T>&); 
+  static vector<D, double> Double(const vector<D, T>&); 
+  
+  T& operator[](const unsigned int);
+ 
+  double dot(const vector<D, T>&) const;
+  static double dot(const vector<D, T>&, const vector<D, T>&);
+  
+  double norm_squared() const;
+  static double norm_squared(const vector<D, T>&);
+
+  void read(std::ifstream&);
+  void write(std::ofstream&) const;
+};
+
+
+template <int D, typename T>
+std::ostream& operator<<(std::ostream&, const vector<D, T>&);
+
+
+// constructor
+// ~~~~~~~~~~~
+template <int D, typename T>
+vector<D, T>::vector()
+{
+  for(int k=0; k<D; k++)
+    x[k] = 0;
+}
+
+template <int D, typename T>
+vector<D, T>::vector(const T x_i[D])
+{
+  for(int k=0; k<D; k++)
+    x[k] = x_i[k];
+}
+
+template <int D, typename T>
+vector<D, T>::vector(const vector<D, T> &v)
+{
+  for(int k=0; k<D; k++)
+    x[k] = v.x[k];
+}
+
+
+// destructor
+// ~~~~~~~~~~
+template <int D, typename T>
+vector<D, T>::~vector()
+{
+}
+
+
+// +=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator+=(const vector<D, T> &v)
+{
+  for(int k=0; k<D; k++)
+    x[k] += v.x[k];
+
+  return *this;
+}
+
+
+// -=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator-=(const vector<D, T> &v)
+{
+  for(int k=0; k<D; k++)
+    x[k] -= v.x[k];
+
+  return *this;
+}
+
+
+// *=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator*=(const T s)
+{
+  for(int k=0; k<D; k++)
+    x[k] *= s;
+
+  return *this;
+}
+
+// /=
+// ~~
+template <int D, typename T>
+inline vector<D, T>& vector<D, T>::operator/=(const T s)
+{
+  for(int k=0; k<D; k++)
+    x[k] /= s;
+
+  return *this;
+}
+
+
+// +
+// ~ 
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator+(const vector<D, T> &a) const
+{
+  vector<D, T> c;
+  
+  for(int k=0; k<D; k++)
+    c.x[k] = x[k] + a.x[k];
+
+  return c;
+}
+
+
+// -
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator-(const vector<D, T> &a) const
+{
+  vector<D, T> c;
+
+  for(int k=0; k<D; k++)
+    c[k] = x[k] - a.x[k];
+
+  return c;
+}
+
+
+// *
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator*(const T s) const
+{
+  vector<D, T> c;
+  
+  for(int k=0; k<D; k++)
+    c[k] = x[k] * s;
+
+  return c;
+}
+
+
+// /
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator/(const T s) const
+{
+  vector<D, T> c;
+
+  for(int k=0; k<D; k++)
+    c[k] = x[k] / s;
+
+  return c;
+}
+
+
+// ==
+// ~
+template <int D, typename T>
+inline bool vector<D, T>::operator==(const vector<D, T> &a) const
+{
+  for(int k=0; k<D; k++)
+    {
+      if (!(x[k]==a.x[k]))
+	return false;
+    }
+  return true;
+}
+
+
+// %
+// ~
+template <int D, typename T>
+inline vector<D, T> vector<D, T>::operator%(const T s) const
+{
+  vector<D, T> c;
+
+  for(int k=0; k<D; k++)
+    c[k] = x[k] % s;
+
+  return c;
+}
+
+
+// integer
+// ~~~~~~~
+template <int D, typename T>
+inline vector<D, int> vector<D, T>::integer() const
+{
+  vector<D, int> c;
+
+  for(int k=0; k<D; k++)
+    c[k] = (int)x[k];
+
+  return c;
+}
+
+template <int D, typename T>
+inline vector<D, int> vector<D, T>::integer(const vector<D, T>& v)
+{
+  return v.integer();
+}
+
+
+// double
+// ~~~~~~~
+template <int D, typename T>
+inline vector<D, double> vector<D, T>::Double() const
+{
+  vector<D, double> c;
+
+  for(int k=0; k<D; k++)
+    c[k] = (double)x[k];
+
+  return c;
+}
+
+template <int D, typename T>
+inline vector<D, double> vector<D, T>::Double(const vector<D, T>& v)
+{
+  return v.Double();
+}
+
+
+
+// []
+// ~~
+template <int D, typename T>
+inline T& vector<D, T>::operator[](const unsigned int i)
+{
+  return x[i];
+}
+
+
+// Dot
+// ~~~
+template <int D, typename T>
+inline double vector<D, T>::dot(const vector<D, T> &a) const
+{
+  double d=0;
+
+  for(int k=0; k<D; k++)
+    d += x[k] * a.x[k];
+
+  return d;
+}
+
+template <int D, typename T>
+inline double vector<D, T>::dot(const vector<D, T> &a, const vector<D, T> &b)
+{
+  return a.dot(b);
+}
+
+
+// NormSquared
+// ~~~~~~~~~~~
+template <int D, typename T>
+inline double vector<D, T>::norm_squared() const
+{
+  return dot(*this, *this);
+}
+
+template <int D, typename T>
+inline double vector<D, T>::norm_squared(const vector<D, T>& v)
+{
+  return v.norm_squared();
+}
+
+
+// read
+// ~~~~
+template <int D, typename T>
+void vector<D, T>::read(std::ifstream& in)
+{
+  in.read((char*)x, sizeof(T)*D);
+}
+
+// write
+// ~~~~~
+template <int D, typename T>
+void vector<D, T>::write(std::ofstream& out) const
+{
+  out.write((const char*)x, sizeof(T)*D);
+}
+
+
+
+// Insertion
+// ~~~~~~~~~
+template <int D, typename T>
+std::ostream& operator<<(std::ostream& os, const vector<D, T>& v)
+{
+  os << "(";
+
+  for(int k=0; k<D-1; k++)
+    os << v.x[k] << ", ";
+
+  os << v.x[D-1] << ")";
+
+  return os;
+}
+
+
+// ======================================================================
+// vector_field 
+// ======================================================================
+
+// A field of V-vectors on a D dimensional manifold
+
+template<int V, int D, typename T=double>
+class vector_field {
+  
+ public:
+  int elements;
+
+ private:
+  vector<V, T>* f;
+  vector<D, int> size;           // number of grid points for each dimension
+  vector<D, int> offset;
+ 
+ public:
+
+  vector_field();
+  vector_field(const vector<D, int>&);
+  ~vector_field();
+
+  vector<D, int> get_size() const;
+  void set_size(const vector<D, int>&);
+
+  vector<V, T>& get(const vector<D, int>&);
+
+  void read(std::ifstream&);
+  void write(std::ofstream&) const;
+
+  static void swap(vector_field<V, D, T>&, vector_field<V, D, T>&);
+};
+
+
+// vector_field
+// ~~~~~~~~~~~~
+template<int V, int D, typename T>
+vector_field<V, D, T>::vector_field()
+  : f(0), elements(0)
+{
+}
+
+
+// vector_field
+// ~~~~~~~~~~~~
+template<int V, int D, typename T>
+vector_field<V, D, T>::vector_field(const vector<D, int>& s)
+  : f(0)
+{
+  set_size(s);
+}
+
+// ~vector_field
+// ~~~~~~~~~~~~~
+template <int V, int D, typename T>
+vector_field<V, D, T>::~vector_field()
+{
+  if(f != 0)
+    delete[] f;
+}
+
+// get_size
+// ~~~~~~~~
+template<int V, int D, typename T>
+inline vector<D, int> vector_field<V, D, T>::get_size() const
+{
+  return size;
+}
+
+// set_size
+// ~~~~~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::set_size(const vector<D, int>& s)
+{
+  if(f != 0)
+    delete[] f;
+
+  size = s;
+
+  elements = 1;
+  for(int i=0; i<D; i++) {
+    offset[i] = elements;
+    elements *= size.x[i];
+  }
+
+  f = new vector<V, T>[elements];
+}
+
+// get
+// ~~~
+template<int V, int D, typename T>
+inline vector<V, T>& vector_field<V, D, T>::get(const vector<D, int>& pos)
+{
+  int p=0;
+  for(int i=0; i<D; i++)
+    p += pos.x[i]*offset[i];
+
+  return f[p];
+}
+
+
+// read
+// ~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::read(std::ifstream& in)
+{
+  in.read((char*)f, elements*sizeof(T)*V);
+}
+
+
+// write
+// ~~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::write(std::ofstream& out) const
+{
+  out.write((const char*)f, elements*sizeof(T)*V);
+}
+
+// swap
+// ~~~~
+template<int V, int D, typename T>
+void vector_field<V, D, T>::swap(vector_field<V, D, T>& v1,
+				 vector_field<V, D, T>& v2)
+{
+  vector<V, T>* f;
+
+  f = v1.f;
+  v1.f = v2.f;
+  v2.f = f;
+}
+
+
+
+#endif