changed code to rely on jst

51 files changed, 242 insertions, 5003 deletions

diff --git a/.gitignore b/.gitignore
index bebedfe..fb544d9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,3 @@
-bin/*
-obj/*
 data/raw/*
 data/raw
-makefile_*
+build/*
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..84b8f0f
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,22 @@
+
+cmake_minimum_required(VERSION 3.7)
+project(fracture)
+
+include_directories(src ~/.local/include)
+link_directories(~/.local/lib)
+
+file(GLOB SOURCES lib/*.c)
+add_executable(fracture src/fracture.c ${SOURCES})
+add_executable(big_anal_process src/big_anal_process.c ${SOURCES})
+add_executable(anal_process src/anal_process.c ${SOURCES})
+add_executable(cen_anal_process src/cen_anal_process.c ${SOURCES})
+add_executable(long_anal_process src/long_anal_process.c ${SOURCES})
+add_executable(corr_test src/corr_test.c ${SOURCES})
+
+target_link_libraries(fracture gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst profiler)
+target_link_libraries(big_anal_process gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler)
+target_link_libraries(cen_anal_process gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler)
+target_link_libraries(long_anal_process gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler)
+target_link_libraries(corr_test gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler)
+target_link_libraries(anal_process gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis m jst tcmalloc profiler)
+
diff --git a/src/bound_set.c b/lib/bound_set.c
index 65f1e6f..65f1e6f 100644
--- a/src/bound_set.c
+++ b/lib/bound_set.c
diff --git a/src/break_edge.c b/lib/break_edge.c
index 2f112c2..2f112c2 100644
--- a/src/break_edge.c
+++ b/lib/break_edge.c
diff --git a/lib/correlations.c b/lib/correlations.c
new file mode 100644
index 0000000..63532ec
--- /dev/null
+++ b/lib/correlations.c
@@ -0,0 +1,47 @@
+
+#include "fracture.h"
+
+double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) {
+	uint_t nv = instance->graph->dnv;
+	uint_t ne = instance->graph->ne;
+	uint_t *ev = instance->graph->dev;
+	bool nulldists = false;
+	if (dists == NULL) {
+		dists = graph_dists(instance->graph, false);
+		nulldists = true;
+	}
+	double *corr = calloc(nv, sizeof(double));
+	uint_t *numat = calloc(nv, sizeof(uint_t));
+
+	for (uint_t i = 0; i < ne; i++) {
+		uint_t v1 = ev[2 * i];
+		uint_t v2 = ev[2 * i + 1];
+		for (uint_t j = 0; j < ne; j++) {
+			uint_t v3 = ev[2 * j];
+			uint_t v4 = ev[2 * j + 1];
+			uint_t dist1 = dists[v1][v3];
+			uint_t dist2 = dists[v1][v4];
+			uint_t dist3 = dists[v2][v3];
+			uint_t dist4 = dists[v2][v4];
+			uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4;
+			corr[dist] += instance->fuses[i] && instance->fuses[j];
+			numat[dist]++;
+		}
+	}
+
+	for (uint_t i = 0; i < nv; i++) {
+		if (numat[i] > 0) {
+			corr[i] /= numat[i];
+		}
+	}
+
+	if (nulldists) {
+		for (int i = 0; i < nv; i++) {
+			free(dists[i]);
+		}
+		free(dists);
+	}
+
+	return corr;
+}
+
diff --git a/src/data.c b/lib/data.c
index 2047c44..2047c44 100644
--- a/src/data.c
+++ b/lib/data.c
diff --git a/src/factor_update.c b/lib/factor_update.c
index a51705a..a51705a 100644
--- a/src/factor_update.c
+++ b/lib/factor_update.c
diff --git a/lib/fracture.h b/lib/fracture.h
new file mode 100644
index 0000000..b1114fb
--- /dev/null
+++ b/lib/fracture.h
@@ -0,0 +1,123 @@
+
+#pragma once
+
+#include <assert.h>
+#include <cholmod.h>
+#include <float.h>
+#include <getopt.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_sf_erf.h>
+#include <gsl/gsl_sf_exp.h>
+#include <gsl/gsl_sf_log.h>
+#include <inttypes.h>
+#include <math.h>
+#include <omp.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <unistd.h>
+
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+
+// these defs allow me to switch to long int cholmod in a sitch
+#define int_t int
+#define uint_t unsigned int
+#define CINT_MAX INT_MAX
+#define CHOL_F(x) cholmod_##x
+
+#define GSL_RAND_GEN gsl_rng_mt19937
+
+typedef struct {
+	const graph_t *graph;
+	bool *fuses;
+	long double *thres;
+	double inf;
+	cholmod_dense *boundary_cond;
+	cholmod_factor *factor;
+	bool voltage_bound;
+	uint_t num_broken;
+	uint_t dim;
+	uint_t nep;
+	uint_t *evp;
+	cholmod_sparse *voltcurmat;
+} net_t;
+
+typedef struct {
+	uint_t num_broken;
+	uint_t *break_list;
+	double *conductivity;
+	long double *extern_field;
+} data_t;
+
+intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
+
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
+
+cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
+
+int edge_to_verts(uint_t width, bool periodic, uint_t edge,
+									bool index);
+
+int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge,
+											 bool index);
+
+double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
+													bool index);
+
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
+void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c);
+
+void net_notch(net_t *net, double notch_len, cholmod_common *c);
+data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
+double *net_voltages(const net_t *net, cholmod_common *c);
+double *net_currents(const net_t *net, const double *voltages, cholmod_common *c);
+double net_conductivity(const net_t *net, const double *voltages);
+
+void update_boundary(net_t *instance, const double *avg_field);
+
+FILE *get_file(const char *prefix, uint_t width, uint_t crack,
+							 double beta, uint_t iter, uint_t num_iter,
+							 uint_t num, bool read);
+
+double update_beta(double beta, uint_t width, const double *stress,
+									 const double *damage, double bound_total);
+
+cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts,
+															 uint_t *edges_to_verts, cholmod_common *c);
+
+net_t *net_copy(const net_t *net, cholmod_common *c);
+
+void net_free(net_t *instance, cholmod_common *c);
+
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
+
+bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
+
+components_t *get_clusters(net_t *instance);
+
+uint_t *get_cluster_dist(net_t *instance);
+
+void randfunc_flat(gsl_rng *r, double *x, double *y);
+void randfunc_gaus(gsl_rng *r, double *x, double *y);
+
+double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c);
+
+double *bin_values(graph_t *network, uint_t width, double *values);
+
+cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c);
+
+data_t *data_create(uint_t num_edges);
+void data_free(data_t *data);
+void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
+
+long double rand_dist_pow(const gsl_rng *r, double beta);
+
+bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/src/gen_laplacian.c b/lib/gen_laplacian.c
index 1cc93a4..1cc93a4 100644
--- a/src/gen_laplacian.c
+++ b/lib/gen_laplacian.c
diff --git a/src/gen_voltcurmat.c b/lib/gen_voltcurmat.c
index e870140..e870140 100644
--- a/src/gen_voltcurmat.c
+++ b/lib/gen_voltcurmat.c
diff --git a/src/geometry.c b/lib/geometry.c
index ec788f1..ec788f1 100644
--- a/src/geometry.c
+++ b/lib/geometry.c
diff --git a/lib/get_dual_clusters.c b/lib/get_dual_clusters.c
new file mode 100644
index 0000000..eaf7562
--- /dev/null
+++ b/lib/get_dual_clusters.c
@@ -0,0 +1,30 @@
+
+#include "fracture.h"
+
+components_t *get_clusters(net_t *instance) {
+	components_t *c = graph_components_get(instance->graph, instance->fuses, true);
+	return c;
+}
+
+unsigned int *get_cluster_dist(net_t *instance) {
+	components_t *c = get_clusters(instance);
+	unsigned int *cluster_dist = (unsigned int *)calloc(
+			instance->graph->dnv, sizeof(unsigned int));
+
+	for (uint32_t i = 1; i <= c->n; i++) {
+		unsigned int num_in_cluster = 0;
+		for (unsigned int j = 0; j < instance->graph->dnv; j++) {
+			if (c->labels[j] == i)
+				num_in_cluster++;
+		}
+
+		if (num_in_cluster == 0)
+			break;
+
+		cluster_dist[num_in_cluster - 1]++;
+	}
+
+	graph_components_free(c);
+
+	return cluster_dist;
+}
diff --git a/src/net.c b/lib/net.c
index a1047e2..eb44721 100644
--- a/src/net.c
+++ b/lib/net.c
@@ -6,7 +6,7 @@ long double *get_thres(uint_t ne, double beta) {
 	assert(thres != NULL);
 
 	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
-	gsl_rng_set(r, rand_seed());
+	gsl_rng_set(r, jst_rand_seed());
 
 	for (uint_t i = 0; i < ne; i++) {
 		thres[i] = rand_dist_pow(r, beta);
@@ -102,6 +102,8 @@ net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, b
 		CHOL_F(free_sparse)(&laplacian, c);
 	}
 
+	net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
+
 	return net;
 }
 
@@ -136,6 +138,7 @@ void net_free(net_t *net, cholmod_common *c) {
 	free(net->thres);
 	CHOL_F(free_dense)(&(net->boundary_cond), c);
 	CHOL_F(free_factor)(&(net->factor), c);
+	CHOL_F(free_sparse)(&(net->voltcurmat), c);
 	free(net->evp);
 	free(net);
 }
diff --git a/src/net_conductivity.c b/lib/net_conductivity.c
index e9325bb..e9325bb 100644
--- a/src/net_conductivity.c
+++ b/lib/net_conductivity.c
diff --git a/src/net_currents.c b/lib/net_currents.c
index 9bb2874..a078336 100644
--- a/src/net_currents.c
+++ b/lib/net_currents.c
@@ -4,7 +4,7 @@
 double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) {
 	uint_t ne = net->graph->ne;
 	uint_t dim = net->graph->nv;
-	cholmod_sparse *voltcurmat = net->graph->voltcurmat;
+	cholmod_sparse *voltcurmat = net->voltcurmat;
 
 	cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
 
diff --git a/src/net_fracture.c b/lib/net_fracture.c
index e7f18fc..e7f18fc 100644
--- a/src/net_fracture.c
+++ b/lib/net_fracture.c
diff --git a/src/net_voltages.c b/lib/net_voltages.c
index c3537a5..7b07201 100644
--- a/src/net_voltages.c
+++ b/lib/net_voltages.c
@@ -16,7 +16,7 @@ double *net_voltages(const net_t *net, cholmod_common *c) {
 	x->x = NULL;
 	CHOL_F(free_dense)(&x, c);
 
-	graph_t *g = net->graph;
+	const graph_t *g = net->graph;
 
 	if (g->boundary == TORUS_BOUND) {
 		return t_voltages;
diff --git a/src/rand.c b/lib/rand.c
index 75722ac..16aeba4 100644
--- a/src/rand.c
+++ b/lib/rand.c
@@ -1,14 +1,6 @@
 
 #include "fracture.h"
 
-unsigned long int rand_seed() {
-	FILE *f = fopen("/dev/urandom", "r");
-	unsigned long int seed;
-	fread(&seed, sizeof(unsigned long int), 1, f);
-	fclose(f);
-	return seed;
-}
-
 long double rand_dist_pow(const gsl_rng *r, double beta) {
 	long double x = 0;
 
diff --git a/makefile b/makefile
deleted file mode 100644
index 9ea3865..0000000
--- a/makefile
+++ /dev/null
@@ -1,45 +0,0 @@
-
-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 frame_create
-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} ${OBJP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%}
-
-opt:
-	@echo build options:
-	@echo "CC       = ${CC}"
-	@echo "CFLAGS   = ${CFLAGS}"
-	@echo "LDFLAGS  = ${LDFLAGS}"
-
-${OBJ:%=obj/%.o} ${BIN:%=obj/%.o}: obj/%.o: src/%.c
-	@echo CC -c -o $@
-	@${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} ${OBJCCP:%=obj/%.o} src/fracture.h
-	@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/%.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/src/corr_test.c b/src/corr_test.c
index b6f8a18..01b3e11 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -7,7 +7,7 @@ int main() {
 
 	unsigned int width = 64;
 
-	graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false, &c);
+	graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false);
 	net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
 	net_fracture(instance, &c, 1e-10);
 
@@ -19,7 +19,7 @@ int main() {
 	printf("\n");
 
 	net_free(instance, &c);
-	graph_free(network, &c);
+	graph_free(network);
 
 	CHOL_F(finish)(&c);
 
diff --git a/src/correlations.c b/src/correlations.c
deleted file mode 100644
index f403ab8..0000000
--- a/src/correlations.c
+++ /dev/null
@@ -1,299 +0,0 @@
-
-#include "fracture.h"
-
-// I implemented a fibonacci heap because wikipedia said it would be fast
-
-struct fibn {
-	uint_t value;
-	uint_t priority;
-	uint_t rank;
-	bool marked;
-	struct fibn *first_child;
-	struct fibn *parent;
-	struct fibn *prev;
-	struct fibn *next;
-};
-
-typedef struct {
-	struct fibn *min;
-	uint_t rank;
-	uint_t trees;
-	struct fibn **vertex_to_node;
-} fibh;
-
-uint_t fib_findmin(fibh *heap) {
-	return heap->min->value;
-}
-
-void ll_setparent(struct fibn *ll, struct fibn *parent) {
-	struct fibn *curnode = ll->next;
-	ll->parent = parent;
-	while (curnode != ll) {
-		curnode->parent = parent;
-		curnode = curnode->next;
-	}
-}
-
-struct fibn *ll_merge(struct fibn *ll1, struct fibn *ll2) {
-	if (ll1 == NULL)
-		return ll2;
-	if (ll2 == NULL)
-		return ll1;
-
-	// link the beginning of list one to the end of list two and vice versa
-	struct fibn *ll1_beg = ll1;
-	struct fibn *ll1_end = ll1->prev;
-	struct fibn *ll2_beg = ll2;
-	struct fibn *ll2_end = ll2->prev;
-
-	ll1_beg->prev = ll2_end;
-	ll1_end->next = ll2_beg;
-	ll2_beg->prev = ll1_end;
-	ll2_end->next = ll1_beg;
-
-	return ll1;
-}
-
-struct fibn *ll_delroot(struct fibn *ll) {
-	if (ll == NULL)
-		return NULL;
-
-	if (ll->next == ll) {
-		return NULL;
-	}
-
-	struct fibn *ll_beg = ll->next;
-	struct fibn *ll_end = ll->prev;
-
-	ll_beg->prev = ll_end;
-	ll_end->next = ll_beg;
-
-	ll->next = ll;
-	ll->prev = ll;
-
-	return ll_beg;
-}
-
-void fib_insert(fibh *heap, uint_t value, uint_t priority) {
-	struct fibn *newnode = calloc(1, sizeof(struct fibn));
-	newnode->value = value;
-	newnode->priority = priority;
-	newnode->next = newnode;
-	newnode->prev = newnode;
-
-	heap->min = ll_merge(heap->min, newnode);
-
-	if (priority < heap->min->priority) {
-		heap->min = newnode;
-	}
-
-	heap->vertex_to_node[value] = newnode;
-
-	heap->trees++;
-}
-
-void fib_deletemin(fibh *heap) {
-	uint_t min_rank = heap->min->rank;
-	struct fibn *min_children = heap->min->first_child;
-
-	struct fibn *trees = ll_delroot(heap->min);
-	heap->vertex_to_node[heap->min->value] = NULL;
-	free(heap->min);
-
-	if (min_children != NULL)
-		ll_setparent(min_children, NULL);
-	trees = ll_merge(trees, min_children);
-
-	heap->trees += min_rank - 1;
-
-	if (trees == NULL) {
-		// min had no children and was only tree, return empty heap
-		heap->min = NULL;
-		heap->rank = 0;
-	} else {
-		// min had children or there were other trees
-		uint_t min_val = UINT_MAX;
-		struct fibn *min_ptr = NULL;
-		uint_t max_rank = 0;
-		struct fibn *curnode = trees;
-		for (uint_t i = 0; i < heap->trees; i++) {
-			if (curnode->priority < min_val) {
-				min_ptr = curnode;
-				min_val = curnode->priority;
-			}
-			if (curnode->rank > max_rank)
-				max_rank = curnode->rank;
-			curnode = curnode->next;
-		}
-
-		if (min_ptr == NULL)
-			min_ptr = trees;
-		heap->min = min_ptr;
-		heap->rank = max_rank;
-
-		struct fibn **rankslots = calloc(max_rank + 1, sizeof(struct fibn *));
-		curnode = heap->min;
-		while (curnode != rankslots[curnode->rank]) {
-			if (rankslots[curnode->rank] == NULL) {
-				rankslots[curnode->rank] = curnode;
-				curnode = curnode->next;
-			} else {
-				struct fibn *oldnode = rankslots[curnode->rank];
-				rankslots[curnode->rank] = NULL;
-				struct fibn *smaller =
-						curnode->priority < oldnode->priority || curnode == heap->min
-								? curnode
-								: oldnode;
-				struct fibn *larger = smaller == curnode ? oldnode : curnode;
-				ll_delroot(larger);
-				ll_setparent(larger, smaller);
-				struct fibn *smaller_children = smaller->first_child;
-				smaller->first_child = ll_merge(smaller_children, larger);
-				heap->trees--;
-				smaller->rank++;
-				if (smaller->rank > heap->rank) {
-					heap->rank = smaller->rank;
-					rankslots =
-							realloc(rankslots, (heap->rank + 1) * sizeof(struct fibn *));
-					rankslots[heap->rank] = NULL;
-				}
-				curnode = smaller;
-			}
-		}
-		free(rankslots);
-	}
-}
-
-void fib_decreasekey(fibh *heap, uint_t value, uint_t new_priority) {
-	struct fibn *node = heap->vertex_to_node[value];
-	if (node != NULL) {
-		node->priority = new_priority;
-		if (node->parent != NULL) {
-			if (node->priority < node->parent->priority) {
-
-				struct fibn *curnode = node;
-				curnode->marked = true;
-				while (curnode->parent != NULL && curnode->marked) {
-					struct fibn *oldparent = curnode->parent;
-					oldparent->rank--;
-					oldparent->first_child = ll_delroot(curnode);
-					ll_setparent(curnode, NULL);
-					ll_merge(heap->min, curnode);
-					heap->trees++;
-					if (curnode->marked) {
-						curnode->marked = false;
-					}
-
-					if (oldparent->marked) {
-						curnode = oldparent;
-					} else {
-						oldparent->marked = true;
-						break;
-					}
-				}
-			}
-		}
-
-		if (new_priority < heap->min->priority) {
-			heap->min = node;
-		}
-	}
-}
-
-uint_t *dijkstra(const graph_t *network, uint_t source) {
-	uint_t nv = network->dnv;
-	uint_t *vei = network->dvei;
-	uint_t *ve = network->dve;
-	uint_t *ev = network->dev;
-
-	uint_t *dist = (uint_t *)calloc(nv, sizeof(uint_t));
-	fibh *Q = (fibh *)calloc(1, sizeof(fibh));
-	Q->vertex_to_node = (struct fibn **)calloc(nv, sizeof(struct fibn *));
-
-	for (uint_t i = 0; i < nv; i++) {
-		if (i != source) {
-			dist[i] = UINT_MAX;
-		}
-
-		fib_insert(Q, i, dist[i]);
-	}
-
-	while (Q->min != NULL) {
-		uint_t u = fib_findmin(Q);
-		fib_deletemin(Q);
-
-		for (uint_t i = 0; i < vei[u + 1] - vei[u]; i++) {
-			uint_t e = ve[vei[u] + i];
-			uint_t v = ev[2 * e] == u ? ev[2 * e + 1] : ev[2 * e];
-			uint_t alt = dist[u] + 1;
-			if (alt < dist[v]) {
-				dist[v] = alt;
-				fib_decreasekey(Q, v, alt);
-			}
-		}
-	}
-
-	free(Q->vertex_to_node);
-	free(Q);
-
-	return dist;
-}
-
-uint_t **get_dists(const graph_t *network) {
-	uint_t nv = network->dnv;
-	uint_t **dists = (uint_t **)malloc(nv * sizeof(uint_t *));
-
-	for (uint_t i = 0; i < nv; i++) {
-		dists[i] = dijkstra(network, i);
-	}
-
-	return dists;
-}
-
-double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) {
-	uint_t nv = instance->graph->dnv;
-	uint_t ne = instance->graph->ne;
-	uint_t *ev = instance->graph->dev;
-	bool nulldists = false;
-	if (dists == NULL) {
-		dists = get_dists(instance->graph);
-		nulldists = true;
-	}
-	double *corr = calloc(nv, sizeof(double));
-	uint_t *marks = get_clusters(instance, c);
-	uint_t *numat = calloc(nv, sizeof(uint_t));
-
-	for (uint_t i = 0; i < ne; i++) {
-		uint_t v1 = ev[2 * i];
-		uint_t v2 = ev[2 * i + 1];
-		for (uint_t j = 0; j < ne; j++) {
-			uint_t v3 = ev[2 * j];
-			uint_t v4 = ev[2 * j + 1];
-			uint_t dist1 = dists[v1][v3];
-			uint_t dist2 = dists[v1][v4];
-			uint_t dist3 = dists[v2][v3];
-			uint_t dist4 = dists[v2][v4];
-			uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4;
-			corr[dist] += instance->fuses[i] && instance->fuses[j];
-			numat[dist]++;
-		}
-	}
-
-	for (uint_t i = 0; i < nv; i++) {
-		if (numat[i] > 0) {
-			corr[i] /= numat[i];
-		}
-	}
-
-	if (nulldists) {
-		for (int i = 0; i < nv; i++) {
-			free(dists[i]);
-		}
-		free(dists);
-	}
-
-	free(marks);
-
-	return corr;
-}
-
diff --git a/src/fitting_functions.c b/src/fitting_functions.c
deleted file mode 100644
index aa7a305..0000000
--- a/src/fitting_functions.c
+++ /dev/null
@@ -1,290 +0,0 @@
-
-#include "fracture.h"
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_multifit_nlinear.h>
-
-struct model_params {
-	double nu;
-	double ac;
-	double mc;
-	double ac02;
-	double ac03;
-	double ac12;
-	double ac13;
-	double bc02;
-	double bc03;
-	double bc12;
-	double bc13;
-	double Ac20;
-	double Ac21;
-	double Ac22;
-	double Ac23;
-	double Ac30;
-	double Ac31;
-	double Ac32;
-	double Ac33;
-	double Bc20;
-	double Bc21;
-	double Bc22;
-	double Bc23;
-	double Bc30;
-	double Bc31;
-	double Bc32;
-	double Bc33;
-};
-
-struct data {
-	uint_t nc;
-	uint_t *Ls_c;
-	double *betas_c;
-	double *vals_c2;
-	double *vals_c3;
-};
-
-double Jc(uint_t n, uint_t o, struct model_params *par, double x) {
-	double nu, ac, mc, ac0n, bc0n, *Acn, yc, sum;
-
-	nu = par->nu;
-	ac = par->ac;
-	mc = par->mc;
-	ac0n = *(&(par->ac02) + (n - 2) * sizeof(double));
-	bc0n = *(&(par->bc02) + (n - 2) * sizeof(double));
-	Acn = &(par->Ac20) + (n - 2) * o * sizeof(double);
-
-	yc = (log(x) - mc) / ac;
-
-	sum = 0;
-
-	for (uint_t i = 0; i < o; i++) {
-		sum += Acn[i] * gsl_sf_laguerre_n(i, 0, yc);
-	}
-
-	return exp(ac0n + bc0n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum);
-}
-
-double Kc(uint_t n, uint_t o, struct model_params *par, double x) {
-	double nu, ac, mc, ac1n, bc1n, *Bcn, yc, sum;
-
-	nu = par->nu;
-	ac = par->ac;
-	mc = par->mc;
-	ac1n = *(&(par->ac12) + (n - 2) * sizeof(double));
-	bc1n = *(&(par->bc12) + (n - 2) * sizeof(double));
-	Bcn = &(par->Bc20) + (n - 2) * o * sizeof(double);
-
-	yc = (log(x) - mc) / ac;
-
-	sum = 0;
-
-	for (uint_t i = 0; i < o; i++) {
-		sum += Bcn[i] * gsl_sf_laguerre_n(i, 0, yc);
-	}
-
-	return exp(ac1n + bc1n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum);
-}
-
-double sc(uint_t n, uint_t o, struct model_params *par, uint_t L, double beta) {
-	double nu, bLnu, deltanu, sigmanu, tau, Lntau, Ldelta;
-
-	nu = par->nu;
-	deltanu = 1.5;
-	sigmanu = 48. / 91;
-	tau = 187. / 91;
-
-	bLnu = beta * exp(log(L) / nu);
-	Lntau = exp((n + 1 - tau) / sigmanu * log(L));
-	Ldelta = exp(-deltanu * log(L));
-
-	return Lntau * (Jc(n, o, par, bLnu) + Ldelta * Kc(n, 0, par, bLnu));
-}
-
-int f_moms(const gsl_vector *x, void *params, gsl_vector *f) {
-	struct data *dat = (struct data *)params;
-
-	for (uint_t i = 0; i < dat->nc; i++) {
-		double f2i = sc(2, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]);
-		double F2i = dat->vals_c2[i];
-
-		f2i = log(f2i);
-		F2i = log(F2i);
-
-		gsl_vector_set(f, i, f2i - F2i);
-
-		double f3i = sc(3, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]);
-		double F3i = dat->vals_c3[i];
-
-		f3i = log(f3i);
-		F3i = log(F3i);
-
-		gsl_vector_set(f, dat->nc + i, f3i - F3i);
-	}
-
-	return GSL_SUCCESS;
-}
-
-void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) {
-	uint_t total = 0;
-	double mom = 0;
-	double mom_err = 0;
-
-	for (uint_t i = 0; i < len; i++) {
-		uint32_t datai = data[i];
-		double in = pow(i, n);
-
-		total += datai;
-		mom += in * datai;
-		mom_err += gsl_pow_2(in) * datai; 
-	}
-
-	double momf = mom / total;
-	double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
-
-	result[0] = momf;
-	result[1] = momf_err;
-}
-
-void
-callback(const size_t iter, void *params,
-         const gsl_multifit_nlinear_workspace *w)
-{
-  gsl_vector *f = gsl_multifit_nlinear_residual(w);
-
-  fprintf(stderr, "iter %2zu: |f(x)| = %.4f\n",
-          iter,
-          gsl_blas_dnrm2(f));
-}
-
-
-int main(int argc, char *argv[]) {
-	struct data *d = (struct data *)malloc(sizeof(struct data));
-	uint_t nc = argc - 1;
-	uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
-	double *betas_c = (double *)malloc(nc * sizeof(double));
-	double *vals_c2 = (double *)malloc(nc * sizeof(double));
-	double *weights_c2 = (double *)malloc(nc * sizeof(double));
-	double *vals_c3 = (double *)malloc(nc * sizeof(double));
-	double *weights_c3 = (double *)malloc(nc * sizeof(double));
-
-	double chisq, chisq0;
-
-	d->nc = nc;
-	d->Ls_c = Ls_c;
-	d->betas_c = betas_c;
-	d->vals_c2 = vals_c2;
-	d->vals_c3 = vals_c3;
-
-	for (uint_t i = 0; i < nc; i++) {
-		char *fn = argv[1 + i];
-		char *buff = malloc(20 * sizeof(char));
-		uint_t pos = 0; uint_t c = 0;
-		while (c < 4) {
-			if (fn[pos] == '_') {
-				c++;
-			}
-			pos++;
-		}
-		c = 0;
-		while (fn[pos] != '_') {
-			buff[c] = fn[pos];
-			pos++;
-			c++;
-		}
-		buff[c] = '\0';
-		Ls_c[i] = atoi(buff);
-		c = 0;
-		pos++;
-		while (fn[pos] != '_') {
-			buff[c] = fn[pos];
-			pos++;
-			c++;
-		}
-		buff[c] = '\0';
-		betas_c[i] = atof(buff);
-		free(buff);
-
-		uint_t dist_len = 4 * pow(Ls_c[i], 2);
-		uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
-		FILE *dist_file = fopen(fn, "rb");
-		fread(dist, sizeof(uint32_t), dist_len, dist_file);
-		fclose(dist_file);
-		{
-			double mom2[2];
-			mom(dist_len, 2, dist, mom2);
-			vals_c2[i] = mom2[0];
-			weights_c2[i] = mom2[1];
-
-			double mom3[2];
-			mom(dist_len, 3, dist, mom3);
-			vals_c3[i] = mom3[0];
-			weights_c3[i] = mom3[1];
-		}
-		free(dist);
-	}
-
-	const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
-	gsl_multifit_nlinear_workspace *w;
-	gsl_multifit_nlinear_fdf fdf;
-	gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters();
-	uint_t n = 2 * nc;
-	uint_t p = 27;
-
-	double x_init[27] = { 1.56, .3, 2, -6, -10, -10, -10, 10, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
-	double weights[n];
-	for (uint_t i = 0; i < nc; i++) {
-		weights[i] = 1/gsl_pow_2(weights_c2[i] / vals_c2[i]);
-		weights[nc + i] = 1/gsl_pow_2(weights_c3[i] / vals_c3[i]);
-	}
-
-	gsl_vector_view x = gsl_vector_view_array(x_init, p);
-	gsl_vector_view wts = gsl_vector_view_array(weights, n);
-
-	gsl_vector *f;
-
-	fdf.f = f_moms;
-	fdf.df = NULL;
-	fdf.fvv = NULL;
-	fdf.n = n;
-	fdf.p = p;
-	fdf.params = d;
-
-	fdf_params.trs = gsl_multifit_nlinear_trs_lm;
-
-	w = gsl_multifit_nlinear_alloc(T, &fdf_params, n, p);
-	gsl_multifit_nlinear_winit(&x.vector, &wts.vector, &fdf, w);
-	f = gsl_multifit_nlinear_residual(w);
-	gsl_blas_ddot(f, f, &chisq0);
-
-	const double xtol = 0.0;
-	const double gtol = 0.0;
-	const double ftol = 0.0;
-
-	int info;
-	int status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol, callback, NULL, &info, w);
-	gsl_blas_ddot(f, f, &chisq);
-
-	fprintf(stderr, "summary from method '%s/%s'\n",
-			          gsl_multifit_nlinear_name(w),
-								          gsl_multifit_nlinear_trs_name(w));
-
-	fprintf(stderr, "number of iterations: %zu\n",
-			          gsl_multifit_nlinear_niter(w));
-	  fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
-			  fprintf(stderr, "reason for stopping: %s\n",
-						          (info == 1) ? "small step size" : "small gradient");
-				  fprintf(stderr, "initial |f(x)| = %g\n", sqrt(chisq0));
-					  fprintf(stderr, "final   |f(x)| = %g\n", sqrt(chisq));
-
-	free(Ls_c);
-	free(betas_c);
-	free(vals_c2);
-	free(weights_c2);
-	free(vals_c3);
-	free(weights_c3);
-
-	return 0;
-}
-
diff --git a/src/fortune/defs.h b/src/fortune/defs.h
deleted file mode 100644
index 9f1b1db..0000000
--- a/src/fortune/defs.h
+++ /dev/null
@@ -1,134 +0,0 @@
-#ifndef NULL
-#define NULL 0
-#endif
-#define DELETED -2
-#include <stdlib.h>
-#include <limits.h>
-#include <stdbool.h>
-
-extern int triangulate, sorted, plot, debug;
-
-struct Freenode {
-	struct Freenode *nextfree;
-};
-struct Freelist {
-	struct Freenode *head;
-	int nodesize;
-};
-char *getfree();
-char *myalloc();
-
-extern double xmin, xmax, ymin, ymax, deltax, deltay;
-
-struct Point {
-	double x, y;
-};
-
-/* structure used both for sites and for vertices */
-struct Site {
-	struct Point coord;
-	int sitenbr;
-	int refcnt;
-};
-
-extern struct Site *sites;
-extern int nsites;
-extern int siteidx;
-extern int sqrt_nsites;
-extern int nvertices;
-extern struct Freelist sfl;
-extern struct Site *bottomsite;
-
-void **alloclist;
-int allocnum;
-
-struct Edge {
-	double a, b, c;
-	struct Site *ep[2];
-	struct Site *reg[2];
-	int edgenbr;
-};
-#define le 0
-#define re 1
-extern int nedges;
-extern struct Freelist efl;
-
-int has_endpoint(), right_of();
-struct Site *intersect();
-double dist();
-struct Point PQ_min();
-struct Halfedge *PQextractmin();
-struct Edge *bisect();
-
-struct Halfedge {
-	struct Halfedge *ELleft, *ELright;
-	struct Edge *ELedge;
-	int ELrefcnt;
-	char ELpm;
-	struct Site *vertex;
-	double ystar;
-	struct Halfedge *PQnext;
-};
-
-extern struct Freelist hfl;
-extern struct Halfedge *ELleftend, *ELrightend;
-extern int ELhashsize;
-extern struct Halfedge **ELhash;
-struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
-struct Site *leftreg(), *rightreg();
-
-extern int PQhashsize;
-extern struct Halfedge *PQhash;
-extern struct Halfedge *PQfind();
-extern int PQcount;
-extern int PQmin;
-int PQempty();
-
-extern double *site_list;
-extern double *vert_list;
-extern double *line_list;
-extern unsigned int *edge_list;
-extern unsigned int *dual_list;
-
-extern int vert_count;
-extern int edge_count;
-extern int line_count;
-extern int dual_count;
-extern int site_count;
-
-void makefree(struct Freenode *curr, struct Freelist *fl);
-void voronoi(int triangulate, struct Site *(*nextsite)());
-char *getfree(struct Freelist *fl);
-char *myalloc(unsigned n);
-void freeinit(struct Freelist *fl, int size);
-void ELinitialize();
-struct Halfedge *HEcreate(struct Edge *e, int pm);
-void ELinsert(struct Halfedge *lb, struct Halfedge *new);
-struct Halfedge *ELgethash(int b);
-struct Halfedge *ELleftbnd(struct Point *p);
-void ELdelete(struct Halfedge *he);
-struct Halfedge *ELright(struct Halfedge *he);
-struct Halfedge *ELleft(struct Halfedge *he);
-struct Site *leftreg(struct Halfedge *he);
-struct Site *rightreg(struct Halfedge *he);
-void geominit();
-struct Edge *bisect(struct Site *s1, struct Site *s2);
-struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2);
-int right_of(struct Halfedge *el, struct Point *p);
-void endpoint(struct Edge *e, int lr, struct Site *s);
-double dist(struct Site *s, struct Site *t);
-void makevertex(struct Site *v);
-void deref(struct Site *v);
-void ref(struct Site *v);
-void PQinsert(struct Halfedge *he, struct Site *v, double offset);
-void PQdelete(struct Halfedge *he);
-int PQbucket(struct Halfedge *he);
-int PQempty();
-struct Point PQ_min();
-struct Halfedge *PQextractmin();
-void PQinitialize();
-void out_bisector(struct Edge *e);
-void out_ep(struct Edge *e);
-void out_vertex(struct Site *v);
-void out_site(struct Site *s);
-void out_triple(struct Site *s1, struct Site *s2, struct Site *s3);
diff --git a/src/fortune/edgelist.c b/src/fortune/edgelist.c
deleted file mode 100644
index dfd85a7..0000000
--- a/src/fortune/edgelist.c
+++ /dev/null
@@ -1,136 +0,0 @@
-#
-#include "defs.h"
-int nedges;
-struct Freelist efl;
-struct Freelist hfl;
-struct Halfedge *ELleftend, *ELrightend;
-int ELhashsize;
-struct Halfedge **ELhash;
-
-int ntry, totalsearch;
-
-void ELinitialize() {
-	int i;
-	freeinit(&hfl, sizeof **ELhash);
-	ELhashsize = 2 * sqrt_nsites;
-	ELhash = (struct Halfedge **)myalloc(sizeof *ELhash * ELhashsize);
-	for (i = 0; i < ELhashsize; i += 1)
-		ELhash[i] = (struct Halfedge *)NULL;
-	ELleftend = HEcreate((struct Edge *)NULL, 0);
-	ELrightend = HEcreate((struct Edge *)NULL, 0);
-	ELleftend->ELleft = (struct Halfedge *)NULL;
-	ELleftend->ELright = ELrightend;
-	ELrightend->ELleft = ELleftend;
-	ELrightend->ELright = (struct Halfedge *)NULL;
-	ELhash[0] = ELleftend;
-	ELhash[ELhashsize - 1] = ELrightend;
-}
-
-struct Halfedge *HEcreate(e, pm) struct Edge *e;
-int pm;
-{
-	struct Halfedge *answer;
-	answer = (struct Halfedge *)getfree(&hfl);
-	answer->ELedge = e;
-	answer->ELpm = pm;
-	answer->PQnext = (struct Halfedge *)NULL;
-	answer->vertex = (struct Site *)NULL;
-	answer->ELrefcnt = 0;
-	return (answer);
-}
-
-void ELinsert(struct Halfedge *lb, struct Halfedge *new) {
-	new->ELleft = lb;
-	new->ELright = lb->ELright;
-	(lb->ELright)->ELleft = new;
-	lb->ELright = new;
-}
-
-/* Get entry from hash table, pruning any deleted nodes */
-struct Halfedge *ELgethash(b) int b;
-{
-	struct Halfedge *he;
-
-	if (b < 0 || b >= ELhashsize)
-		return ((struct Halfedge *)NULL);
-	he = ELhash[b];
-	if (he == (struct Halfedge *)NULL || he->ELedge != (struct Edge *)DELETED)
-		return (he);
-
-	/* Hash table points to deleted half edge.  Patch as necessary. */
-	ELhash[b] = (struct Halfedge *)NULL;
-	if ((he->ELrefcnt -= 1) == 0)
-		makefree((struct Freenode *)he, &hfl);
-	return ((struct Halfedge *)NULL);
-}
-
-struct Halfedge *ELleftbnd(p) struct Point *p;
-{
-	int i, bucket;
-	struct Halfedge *he;
-
-	/* Use hash table to get close to desired halfedge */
-	bucket = (p->x - xmin) / deltax * ELhashsize;
-	if (bucket < 0)
-		bucket = 0;
-	if (bucket >= ELhashsize)
-		bucket = ELhashsize - 1;
-	he = ELgethash(bucket);
-	if (he == (struct Halfedge *)NULL) {
-		for (i = 1; 1; i += 1) {
-			if ((he = ELgethash(bucket - i)) != (struct Halfedge *)NULL)
-				break;
-			if ((he = ELgethash(bucket + i)) != (struct Halfedge *)NULL)
-				break;
-		};
-		totalsearch += i;
-	};
-	ntry += 1;
-	/* Now search linear list of halfedges for the corect one */
-	if (he == ELleftend || (he != ELrightend && right_of(he, p))) {
-		do {
-			he = he->ELright;
-		} while (he != ELrightend && right_of(he, p));
-		he = he->ELleft;
-	} else
-		do {
-			he = he->ELleft;
-		} while (he != ELleftend && !right_of(he, p));
-
-	/* Update hash table and reference counts */
-	if (bucket > 0 && bucket < ELhashsize - 1) {
-		if (ELhash[bucket] != (struct Halfedge *)NULL)
-			ELhash[bucket]->ELrefcnt -= 1;
-		ELhash[bucket] = he;
-		ELhash[bucket]->ELrefcnt += 1;
-	};
-	return (he);
-}
-
-/* This delete routine can't reclaim node, since pointers from hash
-	 table may be present.   */
-void ELdelete(struct Halfedge *he) {
-	(he->ELleft)->ELright = he->ELright;
-	(he->ELright)->ELleft = he->ELleft;
-	he->ELedge = (struct Edge *)DELETED;
-}
-
-struct Halfedge *ELright(he) struct Halfedge *he;
-{ return (he->ELright); }
-
-struct Halfedge *ELleft(he) struct Halfedge *he;
-{ return (he->ELleft); }
-
-struct Site *leftreg(he) struct Halfedge *he;
-{
-	if (he->ELedge == (struct Edge *)NULL)
-		return (bottomsite);
-	return (he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]);
-}
-
-struct Site *rightreg(he) struct Halfedge *he;
-{
-	if (he->ELedge == (struct Edge *)NULL)
-		return (bottomsite);
-	return (he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]);
-}
diff --git a/src/fortune/geometry.c b/src/fortune/geometry.c
deleted file mode 100644
index 131cacf..0000000
--- a/src/fortune/geometry.c
+++ /dev/null
@@ -1,184 +0,0 @@
-#
-#include "defs.h"
-#include <math.h>
-
-void geominit() {
-	struct Edge e;
-	double sn;
-
-	freeinit(&efl, sizeof e);
-	nvertices = 0;
-	nedges = 0;
-	sn = nsites + 4;
-	sqrt_nsites = sqrt(sn);
-	deltay = ymax - ymin;
-	deltax = xmax - xmin;
-}
-
-struct Edge *bisect(s1, s2) struct Site *s1, *s2;
-{
-	double dx, dy, adx, ady;
-	struct Edge *newedge;
-
-	newedge = (struct Edge *)getfree(&efl);
-
-	newedge->reg[0] = s1;
-	newedge->reg[1] = s2;
-	ref(s1);
-	ref(s2);
-	newedge->ep[0] = (struct Site *)NULL;
-	newedge->ep[1] = (struct Site *)NULL;
-
-	dx = s2->coord.x - s1->coord.x;
-	dy = s2->coord.y - s1->coord.y;
-	adx = dx > 0 ? dx : -dx;
-	ady = dy > 0 ? dy : -dy;
-	newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
-	if (adx > ady) {
-		newedge->a = 1.0;
-		newedge->b = dy / dx;
-		newedge->c /= dx;
-	} else {
-		newedge->b = 1.0;
-		newedge->a = dx / dy;
-		newedge->c /= dy;
-	};
-
-	newedge->edgenbr = nedges;
-	out_bisector(newedge);
-	nedges += 1;
-	return (newedge);
-}
-
-struct Site *intersect(el1, el2) struct Halfedge *el1, *el2;
-{
-	struct Edge *e1, *e2, *e;
-	struct Halfedge *el;
-	double d, xint, yint;
-	int right_of_site;
-	struct Site *v;
-
-	e1 = el1->ELedge;
-	e2 = el2->ELedge;
-	if (e1 == (struct Edge *)NULL || e2 == (struct Edge *)NULL)
-		return ((struct Site *)NULL);
-	if (e1->reg[1] == e2->reg[1])
-		return ((struct Site *)NULL);
-
-	d = e1->a * e2->b - e1->b * e2->a;
-	/* printf("intersect: d=%g\n", d); */
-	if (-1.0e-20 < d && d < 1.0e-20) {
-		return ((struct Site *)NULL);
-	};
-
-	xint = (e1->c * e2->b - e2->c * e1->b) / d;
-	yint = (e2->c * e1->a - e1->c * e2->a) / d;
-
-	if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
-			(e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
-			 e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
-		el = el1;
-		e = e1;
-	} else {
-		el = el2;
-		e = e2;
-	};
-	right_of_site = xint >= e->reg[1]->coord.x;
-	if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re))
-		return ((struct Site *)NULL);
-
-	v = (struct Site *)getfree(&sfl);
-	v->refcnt = 0;
-	v->coord.x = xint;
-	v->coord.y = yint;
-	return (v);
-}
-
-/* returns 1 if p is to right of halfedge e */
-int right_of(el, p) struct Halfedge *el;
-struct Point *p;
-{
-	struct Edge *e;
-	struct Site *topsite;
-	int right_of_site, above, fast;
-	double dxp, dyp, dxs, t1, t2, t3, yl;
-
-	e = el->ELedge;
-	topsite = e->reg[1];
-	right_of_site = p->x > topsite->coord.x;
-	if (right_of_site && el->ELpm == le)
-		return (1);
-	if (!right_of_site && el->ELpm == re)
-		return (0);
-
-	if (e->a == 1.0) {
-		dyp = p->y - topsite->coord.y;
-		dxp = p->x - topsite->coord.x;
-		fast = 0;
-		if ((!right_of_site && e->b < 0.0) | (right_of_site && e->b >= 0.0)) {
-			above = dyp >= e->b * dxp;
-			fast = above;
-		} else {
-			above = p->x + p->y * e->b > e->c;
-			if (e->b < 0.0)
-				above = !above;
-			if (!above)
-				fast = 1;
-		};
-		if (!fast) {
-			dxs = topsite->coord.x - (e->reg[0])->coord.x;
-			above = e->b * (dxp * dxp - dyp * dyp) <
-							dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
-			if (e->b < 0.0)
-				above = !above;
-		};
-	} else /*e->b==1.0 */
-			{
-		yl = e->c - e->a * p->x;
-		t1 = p->y - yl;
-		t2 = p->x - topsite->coord.x;
-		t3 = yl - topsite->coord.y;
-		above = t1 * t1 > t2 * t2 + t3 * t3;
-	};
-	return (el->ELpm == le ? above : !above);
-}
-
-void endpoint(e, lr, s) struct Edge *e;
-int lr;
-struct Site *s;
-{
-	e->ep[lr] = s;
-	ref(s);
-	if (e->ep[re - lr] == (struct Site *)NULL) {
-	} else {
-		out_ep(e);
-		deref(e->reg[le]);
-		deref(e->reg[re]);
-		makefree((struct Freenode *)e, &efl);
-	}
-}
-
-double dist(s, t) struct Site *s, *t;
-{
-	double dx, dy;
-	dx = s->coord.x - t->coord.x;
-	dy = s->coord.y - t->coord.y;
-	return (sqrt(dx * dx + dy * dy));
-}
-
-void makevertex(v) struct Site *v;
-{
-	v->sitenbr = nvertices;
-	nvertices += 1;
-	out_vertex(v);
-}
-
-void deref(v) struct Site *v;
-{
-	v->refcnt -= 1;
-	if (v->refcnt == 0)
-		makefree((struct Freenode *)v, &sfl);
-}
-
-void ref(v) struct Site *v;
-{ v->refcnt += 1; }
diff --git a/src/fortune/heap.c b/src/fortune/heap.c
deleted file mode 100644
index 2ebbd58..0000000
--- a/src/fortune/heap.c
+++ /dev/null
@@ -1,94 +0,0 @@
-#
-#include "defs.h"
-int PQhashsize;
-struct Halfedge *PQhash;
-struct Halfedge *PQfind();
-int PQcount;
-int PQmin;
-int PQempty();
-
-void PQinsert(he, v, offset) struct Halfedge *he;
-struct Site *v;
-double offset;
-{
-	struct Halfedge *last, *next;
-	he->vertex = v;
-	ref(v);
-	he->ystar = v->coord.y + offset;
-	last = &PQhash[PQbucket(he)];
-	while ((next = last->PQnext) != (struct Halfedge *)NULL &&
-				 (he->ystar > next->ystar ||
-					(he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
-		last = next;
-	};
-	he->PQnext = last->PQnext;
-	last->PQnext = he;
-	PQcount += 1;
-}
-
-void PQdelete(he) struct Halfedge *he;
-{
-	struct Halfedge *last;
-
-	if (he->vertex != (struct Site *)NULL) {
-		last = &PQhash[PQbucket(he)];
-		while (last->PQnext != he)
-			last = last->PQnext;
-		last->PQnext = he->PQnext;
-		PQcount -= 1;
-		deref(he->vertex);
-		he->vertex = (struct Site *)NULL;
-	};
-}
-
-int PQbucket(he) struct Halfedge *he;
-{
-	int bucket;
-
-	if (he->ystar < ymin)
-		bucket = 0;
-	else if (he->ystar >= ymax)
-		bucket = PQhashsize - 1;
-	else
-		bucket = (he->ystar - ymin) / deltay * PQhashsize;
-	if (bucket < 0)
-		bucket = 0;
-	if (bucket >= PQhashsize)
-		bucket = PQhashsize - 1;
-	if (bucket < PQmin)
-		PQmin = bucket;
-	return (bucket);
-}
-
-int PQempty() { return (PQcount == 0); }
-
-struct Point PQ_min() {
-	struct Point answer;
-
-	while (PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {
-		PQmin += 1;
-	};
-	answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
-	answer.y = PQhash[PQmin].PQnext->ystar;
-	return (answer);
-}
-
-struct Halfedge *PQextractmin() {
-	struct Halfedge *curr;
-	curr = PQhash[PQmin].PQnext;
-	PQhash[PQmin].PQnext = curr->PQnext;
-	PQcount -= 1;
-	return (curr);
-}
-
-void PQinitialize() {
-	int i;
-//	struct Point *s;
-
-	PQcount = 0;
-	PQmin = 0;
-	PQhashsize = 4 * sqrt_nsites;
-	PQhash = (struct Halfedge *)myalloc(PQhashsize * sizeof *PQhash);
-	for (i = 0; i < PQhashsize; i += 1)
-		PQhash[i].PQnext = (struct Halfedge *)NULL;
-}
diff --git a/src/fortune/main.c b/src/fortune/main.c
deleted file mode 100644
index b6b4476..0000000
--- a/src/fortune/main.c
+++ /dev/null
@@ -1,287 +0,0 @@
-#
-#include <stdio.h>
-#include <stdint.h>
-#include "defs.h"
-#include <math.h>
-struct Site *nextone();
-
-int triangulate, sorted, plot, debug;
-double xmin, xmax, ymin, ymax, deltax, deltay;
-int vert_count, edge_count, line_count, dual_count, site_count;
-double *site_list, *vert_list, *line_list;
-unsigned int *edge_list, *dual_list;
-
-/* sort sites on y, then x, coord */
-int scomp(s1, s2) struct Point *s1, *s2;
-{
-	if (s1->y < s2->y)
-		return (-1);
-	if (s1->y > s2->y)
-		return (1);
-	if (s1->x < s2->x)
-		return (-1);
-	if (s1->x > s2->x)
-		return (1);
-	return (0);
-}
-
-/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
-void readsites(double *lattice) {
-
-	sites = (struct Site *)myalloc(nsites * sizeof *sites);
-	for (unsigned int j = 0; j < nsites; j++) {
-		sites[j].coord.x = lattice[2 * j];
-		sites[j].coord.y = lattice[2 * j + 1];
-		sites[j].sitenbr = j;
-		sites[j].refcnt = 0;
-	}
-	qsort(sites, nsites, sizeof(struct Site), scomp);
-	xmin = sites[0].coord.x;
-	xmax = sites[0].coord.x;
-	for (unsigned int i = 1; i < nsites; i += 1) {
-		if (sites[i].coord.x < xmin)
-			xmin = sites[i].coord.x;
-		if (sites[i].coord.x > xmax)
-			xmax = sites[i].coord.x;
-	};
-	ymin = sites[0].coord.y;
-	ymax = sites[nsites - 1].coord.y;
-}
-
-unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) {
-	unsigned int ndup = 0;
-	bool *duplicates = (bool *)calloc(ne, sizeof(bool));
-	for (unsigned int i = 0; i < ne; i++) {
-		unsigned int v1 = etv[2 * i];
-		unsigned int v2 = etv[2 * i + 1];
-		for (unsigned int j = 0; j < i; j++) {
-			unsigned int w1 = etv[2 * j];
-			unsigned int w2 = etv[2 * j + 1];
-			if (w1 == v1 && w2 == v2) {
-				duplicates[j] = true;
-				ndup++;
-				break;
-			}
-		}
-	}
-	unsigned int newsize = (int)ne - (int)ndup;
-	unsigned int count = 0;
-	for (unsigned int i = 0; i < ne; i++) {
-		if (!duplicates[i]) {
-			etv[2 * count] = etv[2 * i];
-			etv[2 * count + 1] = etv[2 * i + 1];
-			count++;
-		}
-	}
-	free(duplicates);
-	return newsize;
-}
-
-intptr_t *run_voronoi(unsigned int num, double *lattice, bool periodic, double xmin, double xmax, double ymin, double ymax)
-{
-	struct Site *(*next)();
-
-	sorted = 0;
-	triangulate = 0;
-	plot = 0;
-	debug = 0;
-
-	alloclist = (void **)malloc(9 * num * sizeof(void *));
-	allocnum = 0;
-
-	freeinit(&sfl, sizeof *sites);
-
-	unsigned int eff_num = num;
-	double *eff_lattice = lattice;
-
-	if (periodic) {
-		eff_num = 9 * num;
-		eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
-
-		for (unsigned int i = 0; i < num; i++) {
-			// original sites - our baby boys
-			eff_lattice[2*i] = lattice[2*i];
-			eff_lattice[2*i+1] = lattice[2*i+1];
-
-			// sites to the right
-			eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax;
-			eff_lattice[2*(num+i)] = lattice[2*i];
-
-			// sites to the left
-			eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin;
-			eff_lattice[2*(2*num+i)] = lattice[2*i];
-
-			// sites to the top
-			eff_lattice[2*(3*num+i)+1] = lattice[2*i+1];
-			eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax;
-
-			// sites to the bottom
-			eff_lattice[2*(4*num+i)+1] = lattice[2*i+1];
-			eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin;
-
-			// sites to the upper right
-			eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax;
-			eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax;
-
-			// sites to the upper left
-			eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin;
-			eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax;
-
-			// sites to the lower left
-			eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin;
-			eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax;
-
-			// sites to the lower right
-			eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin;
-			eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax;
-		}
-	}
-
-	nsites = eff_num;
-	readsites(eff_lattice);
-	if (periodic) free(eff_lattice);
-	next = nextone;
-
-	siteidx = 0;
-	geominit();
-
-	site_list = malloc(2 * nsites * sizeof(double));
-	vert_list = NULL;
-	line_list = NULL;
-	edge_list = NULL;
-	dual_list = NULL;
-	site_count = 0;
-	vert_count = 0;
-	edge_count = 0;
-	line_count = 0;
-	dual_count = 0;
-
-	voronoi(triangulate, next);
-
-	unsigned int real_vert_count = vert_count;
-	unsigned int real_edge_count = edge_count;
-	double *real_vert_list = vert_list;
-	unsigned int *real_edge_list = edge_list;
-	unsigned int *real_dual_list = dual_list;
-
-	if (periodic) {
-		real_vert_count = 0;
-		real_edge_count = 0;
-		real_vert_list = (double *)malloc(2 * vert_count * sizeof(double));
-		real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int));
-		real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int));
-		unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int));
-		unsigned int q = 0;
-		int *new_vert_ind = (int *)malloc(vert_count * sizeof(int));
-
-		for (unsigned int i = 0; i < vert_count; i++) {
-			unsigned int t1, t2, t3;
-			t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2];
-			if (t1 >= num && t2 >= num && t3 >= num) {
-				new_vert_ind[i] = -1;
-			}
-			else if (t1 < num && t2 < num && t3 < num) {
-				real_vert_list[2*real_vert_count] = vert_list[2*i];
-				real_vert_list[2*real_vert_count+1] = vert_list[2*i+1];
-				real_dual_list[3*real_vert_count] = t1;
-				real_dual_list[3*real_vert_count+1] = t2;
-				real_dual_list[3*real_vert_count+2] = t3;
-				new_vert_ind[i] = real_vert_count;
-				real_vert_count++;
-			}
-			else {
-				unsigned int tt1, tt2, tt3;
-				tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num;
-				unsigned int s1, s2, s3;
-				s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3);
-				s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3);
-				s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2));
-
-				bool matched = false;
-				unsigned int ass_vert;
-
-				for (unsigned int j = 0; j < q; j++) {
-					if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) {
-						matched = true;
-						ass_vert = triangle_analyzed[4*j+3];
-						break;
-					}
-				}
-
-				if (matched) {
-					new_vert_ind[i] = ass_vert;
-				} else {
-					triangle_analyzed[4*q] = s1;
-					triangle_analyzed[4*q+1] = s2;
-					triangle_analyzed[4*q+2] = s3;
-					triangle_analyzed[4*q+3] = real_vert_count;
-					q++;
-					real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax);
-					real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax);
-					real_dual_list[3*real_vert_count] = t1 % num;
-					real_dual_list[3*real_vert_count+1] = t2 % num;
-					real_dual_list[3*real_vert_count+2] = t3 % num;
-					new_vert_ind[i] = real_vert_count;
-					real_vert_count++;
-				}
-			}
-		}
-		for (unsigned int i = 0; i < edge_count; i++) {
-			unsigned int v1, v2;
-			v1 = edge_list[2 * i];
-			v2 = edge_list[2 * i+1];
-			if (v1 != UINT_MAX && v2 != UINT_MAX) {
-				if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) {
-					unsigned int nv1 = new_vert_ind[v1];
-					unsigned int nv2 = new_vert_ind[v2];
-					real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2;
-					real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1;
-					real_edge_count++;
-				}
-			}
-		}
-		unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list);
-		real_edge_count = new_edge_count;
-		free(triangle_analyzed);
-		free(new_vert_ind);
-		free(dual_list);
-		free(edge_list);
-		free(vert_list);
-	}
-
-
-	intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t));
-	output[0] = (intptr_t)malloc(4 * sizeof(unsigned int));
-	((unsigned int *)output[0])[0] = real_vert_count;
-	((unsigned int *)output[0])[1] = real_edge_count;
-	((unsigned int *)output[0])[2] = line_count;
-	((unsigned int *)output[0])[3] = real_vert_count;
-	output[1] = (intptr_t)site_list;
-	output[2] = (intptr_t)real_vert_list;
-	output[3] = (intptr_t)real_edge_list;
-	output[4] = (intptr_t)line_list;
-	output[5] = (intptr_t)real_dual_list;
-
-	free(ELhash);
-	free(PQhash);
-	free(sites);
-	for (unsigned int i = 0; i < allocnum; i++) {
-		free(alloclist[i]);
-	}
-	free(alloclist);
-
-	return output;
-}
-
-/* return a single in-storage site */
-struct Site *nextone() {
-	for (; siteidx < nsites; siteidx += 1) {
-		if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x ||
-				sites[siteidx].coord.y != sites[siteidx - 1].coord.y) {
-			siteidx += 1;
-			return (&sites[siteidx - 1]);
-		};
-	};
-	return ((struct Site *)NULL);
-}
-
diff --git a/src/fortune/memory.c b/src/fortune/memory.c
deleted file mode 100644
index 3d62a92..0000000
--- a/src/fortune/memory.c
+++ /dev/null
@@ -1,44 +0,0 @@
-#
-#include "defs.h"
-#include <stdio.h>
-
-void freeinit(struct Freelist *fl, int size) {
-	fl->head = (struct Freenode *)NULL;
-	fl->nodesize = size;
-}
-
-char *getfree(fl) struct Freelist *fl;
-{
-	int i;
-	struct Freenode *t;
-	if (fl->head == (struct Freenode *)NULL) {
-		t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize);
-		alloclist[allocnum] = (void *)t;
-		allocnum++;
-		for (i = 0; i < sqrt_nsites; i += 1)
-			makefree((struct Freenode *)((char *)t + i * fl->nodesize), fl);
-	};
-	t = fl->head;
-	fl->head = (fl->head)->nextfree;
-	return ((char *)t);
-}
-
-void makefree(curr, fl) struct Freenode *curr;struct Freelist *fl;
-{
-	curr->nextfree = fl->head;
-	fl->head = curr;
-}
-
-int total_alloc;
-char *myalloc(n) unsigned n;
-{
-	char *t;
-	if ((t = malloc(n)) == (char *)0) {
-		fprintf(stderr,
-						"Insufficient memory processing site %d (%d bytes in use)\n",
-						siteidx, total_alloc);
-		exit(1);
-	};
-	total_alloc += n;
-	return (t);
-}
diff --git a/src/fortune/output.c b/src/fortune/output.c
deleted file mode 100644
index d496feb..0000000
--- a/src/fortune/output.c
+++ /dev/null
@@ -1,46 +0,0 @@
-#
-#include "defs.h"
-#include <stdio.h>
-double pxmin, pxmax, pymin, pymax, cradius;
-
-void out_bisector(e) struct Edge *e;
-{
-	if (line_count % nsites == 0) line_list = realloc(line_list, 3 * (nsites + line_count) * sizeof(double));
-	line_list[3*line_count] = e->a;
-	line_list[3*line_count+1] = e->b;
-	line_list[3*line_count+2] = e->c;
-	line_count++;
-}
-
-void out_ep(e) struct Edge *e;
-{
-	if (edge_count % nsites == 0) edge_list = realloc(edge_list, 2 * (nsites + edge_count) * sizeof(unsigned int));
-	edge_list[2 * edge_count] = (e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : UINT_MAX);
-	edge_list[2 * edge_count + 1] = (e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : UINT_MAX);
-	edge_count++;
-}
-
-void out_vertex(v) struct Site *v;
-{
-	if (vert_count % nsites == 0) vert_list = realloc(vert_list, 2 * (nsites + vert_count) * sizeof(double));
-	vert_list[2 * vert_count] = v->coord.x;
-	vert_list[2 * vert_count + 1] = v->coord.y;
-	vert_count++;
-}
-
-void out_site(s) struct Site *s;
-{
-	site_list[2 * site_count] = s->coord.x;
-	site_list[2 * site_count + 1] = s->coord.y;
-	site_count++;
-}
-
-void out_triple(s1, s2, s3) struct Site *s1, *s2, *s3;
-{
-	if (dual_count % nsites == 0) dual_list = realloc(dual_list, 3 * (nsites + dual_count) * sizeof(unsigned int));
-	dual_list[dual_count * 3] = s1->sitenbr;
-	dual_list[dual_count * 3 + 1] = s2->sitenbr;
-	dual_list[dual_count * 3 + 2] = s3->sitenbr;
-	dual_count++;
-}
-
diff --git a/src/fortune/voronoi.c b/src/fortune/voronoi.c
deleted file mode 100644
index dc9945f..0000000
--- a/src/fortune/voronoi.c
+++ /dev/null
@@ -1,103 +0,0 @@
-#
-#include "defs.h"
-
-struct Site *sites;
-int nsites;
-int siteidx;
-int sqrt_nsites;
-int nvertices;
-struct Freelist sfl;
-struct Site *bottomsite;
-
-/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
-	 deltax, deltay (can all be estimates).
-	 Performance suffers if they are wrong; better to make nsites,
-	 deltax, and deltay too big than too small.  (?) */
-
-void voronoi(triangulate, nextsite) int triangulate;
-struct Site *(*nextsite)();
-{
-	struct Site *newsite, *bot, *top, *temp, *p;
-	struct Site *v;
-	struct Point newintstar;
-	int pm;
-	struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
-	struct Edge *e;
-
-	PQinitialize();
-	bottomsite = (*nextsite)();
-	out_site(bottomsite);
-	ELinitialize();
-
-	newsite = (*nextsite)();
-	while (1) {
-		if (!PQempty())
-			newintstar = PQ_min();
-
-		if (newsite != (struct Site *)NULL &&
-				(PQempty() || newsite->coord.y < newintstar.y ||
-				 (newsite->coord.y == newintstar.y &&
-					newsite->coord.x < newintstar.x))) {/* new site is smallest */
-			out_site(newsite);
-			lbnd = ELleftbnd(&(newsite->coord));
-			rbnd = ELright(lbnd);
-			bot = rightreg(lbnd);
-			e = bisect(bot, newsite);
-			bisector = HEcreate(e, le);
-			ELinsert(lbnd, bisector);
-			if ((p = intersect(lbnd, bisector)) != (struct Site *)NULL) {
-				PQdelete(lbnd);
-				PQinsert(lbnd, p, dist(p, newsite));
-			};
-			lbnd = bisector;
-			bisector = HEcreate(e, re);
-			ELinsert(lbnd, bisector);
-			if ((p = intersect(bisector, rbnd)) != (struct Site *)NULL) {
-				PQinsert(bisector, p, dist(p, newsite));
-			};
-			newsite = (*nextsite)();
-		} else if (!PQempty())
-				/* intersection is smallest */
-				{
-			lbnd = PQextractmin();
-			llbnd = ELleft(lbnd);
-			rbnd = ELright(lbnd);
-			rrbnd = ELright(rbnd);
-			bot = leftreg(lbnd);
-			top = rightreg(rbnd);
-			out_triple(bot, top, rightreg(lbnd));
-			v = lbnd->vertex;
-			makevertex(v);
-			endpoint(lbnd->ELedge, lbnd->ELpm, v);
-			endpoint(rbnd->ELedge, rbnd->ELpm, v);
-			ELdelete(lbnd);
-			PQdelete(rbnd);
-			ELdelete(rbnd);
-			pm = le;
-			if (bot->coord.y > top->coord.y) {
-				temp = bot;
-				bot = top;
-				top = temp;
-				pm = re;
-			}
-			e = bisect(bot, top);
-			bisector = HEcreate(e, pm);
-			ELinsert(llbnd, bisector);
-			endpoint(e, re - pm, v);
-			deref(v);
-			if ((p = intersect(llbnd, bisector)) != (struct Site *)NULL) {
-				PQdelete(llbnd);
-				PQinsert(llbnd, p, dist(p, bot));
-			};
-			if ((p = intersect(bisector, rrbnd)) != (struct Site *)NULL) {
-				PQinsert(bisector, p, dist(p, bot));
-			};
-		} else
-			break;
-	};
-
-	for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
-		e = lbnd->ELedge;
-		out_ep(e);
-	};
-}
diff --git a/src/fracture.c b/src/fracture.c
index 6ad2f26..fd9de61 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -261,7 +261,7 @@ int main(int argc, char *argv[]) {
 	for (uint32_t i = 0; i < N; i++) {
 		printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
 
-		graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
+		graph_t *g = graph_create(lattice, boundary, L, use_dual);
 		net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
 		net_t *tmp_net = net_copy(net, &c);
 		data_t *data = net_fracture(tmp_net, &c, cutoff);
@@ -394,7 +394,7 @@ int main(int argc, char *argv[]) {
 		}
 
 		if (save_cluster_dist) {
-			uint_t *tmp_cluster_dist = get_cluster_dist(net, &c);
+			uint_t *tmp_cluster_dist = get_cluster_dist(net);
 			for (uint_t j = 0; j < g->dnv; j++) {
 				cluster_size_dist[j] += tmp_cluster_dist[j];
 			}
@@ -403,7 +403,7 @@ int main(int argc, char *argv[]) {
 
 
 		net_free(net, &c);
-		graph_free(g, &c);
+		graph_free(g);
 	}
 
 	printf("\033[F\033[JFRACTURE: COMPLETE\n");
diff --git a/src/fracture.h b/src/fracture.h
index 0a3a687..b1114fb 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -23,6 +23,10 @@
 #include <sys/types.h>
 #include <unistd.h>
 
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+
 // these defs allow me to switch to long int cholmod in a sitch
 #define int_t int
 #define uint_t unsigned int
@@ -31,54 +35,6 @@
 
 #define GSL_RAND_GEN gsl_rng_mt19937
 
-typedef enum lattice_t {
-	VORONOI_LATTICE,
-	SQUARE_LATTICE,
-	VORONOI_HYPERUNIFORM_LATTICE
-} lattice_t;
-
-typedef enum bound_t {
-	FREE_BOUND,
-	CYLINDER_BOUND,
-	TORUS_BOUND,
-	EMBEDDED_BOUND
-} bound_t;
-
-typedef struct {
-	uint_t ne;
-	uint_t nv;
-	uint_t dnv;
-	uint_t *ev;
-	uint_t *dev;
-	double *vx;
-	double *dvx;
-} frame_t;
-
-typedef struct {
-	uint_t L;
-	bound_t boundary;
-	uint_t ne;
-	uint_t nv;
-	uint_t dnv;
-	uint_t nb;
-	uint_t *ev; 
-	uint_t *dev;
-	double *vx;
-	double *dvx;
-	uint_t *vei;
-	uint_t *ve;
-	uint_t *dvei;
-	uint_t *dve;
-	uint_t *bi;
-	uint_t *b;
-	uint_t *nbi;
-	uint_t *bni;
-	bool *bq;
-	uint_t *spanning_edges;
-	uint_t num_spanning_edges;
-	cholmod_sparse *voltcurmat;
-} graph_t;
-
 typedef struct {
 	const graph_t *graph;
 	bool *fuses;
@@ -91,6 +47,7 @@ typedef struct {
 	uint_t dim;
 	uint_t nep;
 	uint_t *evp;
+	cholmod_sparse *voltcurmat;
 } net_t;
 
 typedef struct {
@@ -102,11 +59,6 @@ typedef struct {
 
 intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
 
-int update_components(const cholmod_sparse *laplacian, uint_t *marks,
-											int old_num_components, int v1, int v2, int exclude);
-
-uint_t *find_components(const cholmod_sparse *laplacian, uint_t skip);
-
 cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
 
 cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
@@ -121,6 +73,7 @@ double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
 													bool index);
 
 void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
+void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c);
 
 void net_notch(net_t *net, double notch_len, cholmod_common *c);
 data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
@@ -142,35 +95,19 @@ cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts,
 
 net_t *net_copy(const net_t *net, cholmod_common *c);
 
-graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
-												 cholmod_common *c);
-
-void graph_free(graph_t *network, cholmod_common *c);
 void net_free(net_t *instance, cholmod_common *c);
 
 net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
 
-graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
-													double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
-													cholmod_common *c);
-
 bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
 
-void finish_instance(net_t *instance, cholmod_common *c);
-
-net_t *coursegrain_square(net_t *instance, graph_t *network_p, cholmod_common *c);
-
-uint_t *get_clusters(net_t *instance, cholmod_common *c);
+components_t *get_clusters(net_t *instance);
 
-uint_t *get_cluster_dist(net_t *instance, cholmod_common *c);
+uint_t *get_cluster_dist(net_t *instance);
 
-double *genfunc_uniform(uint_t num, gsl_rng *r);
-double *genfunc_hyperuniform(uint_t num, gsl_rng *r);
 void randfunc_flat(gsl_rng *r, double *x, double *y);
 void randfunc_gaus(gsl_rng *r, double *x, double *y);
 
-uint_t **get_dists(const graph_t *network);
-uint_t *dijkstra(const graph_t *g, uint_t source);
 double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c);
 
 double *bin_values(graph_t *network, uint_t width, double *values);
@@ -181,18 +118,6 @@ data_t *data_create(uint_t num_edges);
 void data_free(data_t *data);
 void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
 
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
-
-uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const uint_t *vei, const uint_t *ve, int **cycles);
-
-bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, int label, int stop_at, int exclude);
-
-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);
-
-frame_t *frame_create(lattice_t lattice, uint_t L, bool dual);
-void frame_free(frame_t *frame);
 bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/src/frame_create.c b/src/frame_create.c
deleted file mode 100644
index 89ce69d..0000000
--- a/src/frame_create.c
+++ /dev/null
@@ -1,167 +0,0 @@
-
-#include "fracture.h"
-
-uint_t *get_voro_dual_edges(uint_t num_edges,
-																	uint_t num_verts, uint_t *edges,
-																	uint_t *triangles) {
-	uint_t *dual_edges =
-			(uint_t *)malloc(2 * num_edges * sizeof(uint_t));
-	uint_t place = 0;
-	for (uint_t i = 0; i < num_edges; i++) {
-		uint_t v1, v2;
-		v1 = edges[2 * i];
-		v2 = edges[2 * i + 1];
-		if (v1 < num_verts && v2 < num_verts) {
-			bool found_match = false;
-			for (uint_t j = 0; j < 3; j++) {
-				for (uint_t k = 0; k < 3; k++) {
-					uint_t t11, t12, t21, t22;
-					t11 = triangles[3 * v1 + j];
-					t12 = triangles[3 * v1 + ((j + 1) % 3)];
-					t21 = triangles[3 * v2 + k];
-					t22 = triangles[3 * v2 + ((k + 1) % 3)];
-					if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
-						dual_edges[2 * place] = t11 < t12 ? t11 : t12;
-						dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
-						place++;
-						found_match = true;
-						break;
-					}
-				}
-				if (found_match)
-					break;
-			}
-		}
-	}
-
-	return dual_edges;
-}
-
-frame_t *frame_create_voronoi(uint_t L, bool dual, bool hyperuniform) {
-	double *dvx = NULL;
-	uint_t dnv = 2 * pow(L / 2, 2);
-
-	{
-		gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
-		gsl_rng_set(r, rand_seed());
-
-		if (hyperuniform) {
-			dvx = genfunc_hyperuniform(dnv, r);
-		} else {
-			dvx = genfunc_uniform(dnv, r);
-		}
-
-		gsl_rng_free(r);
-	}
-
-	// retrieve a periodic voronoi tesselation of the lattice
-	intptr_t *vout = run_voronoi(dnv, dvx, true, 0, 1, 0, 1);
-
-	uint_t nv = ((uint_t *)vout[0])[0];
-	uint_t ne = ((uint_t *)vout[0])[1];
-	double *vx = (double *)vout[2];
-	uint_t *ev = (uint_t *)vout[3];
-	uint_t *voro_tris = (uint_t *)vout[5];
-
-	free((void *)vout[0]);
-	free((void *)vout[1]);
-	free((void *)vout[4]);
-	free(vout);
-
-	// get dual edges of the fully periodic graph
-	uint_t *dev = get_voro_dual_edges(ne, nv, ev, voro_tris);
-
-	frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
-	frame->ne = ne;
-
-	// when use_dual is specificed, the edge and vertex sets are swapped with the
-	// dual edge and dual vertex sets.  once formally relabelled, everything
-	// works the same way
-	if (dual) {
-		frame->nv = dnv;
-		frame->dnv = nv;
-		frame->ev = dev;
-		frame->dev = ev;
-		frame->vx = dvx;
-		frame->dvx = vx;
-	} else {
-		frame->nv = nv;
-		frame->dnv = dnv;
-		frame->ev = ev;
-		frame->dev = dev;
-		frame->vx = vx;
-		frame->dvx = dvx;
-	}
-
-	return frame;
-}
-
-frame_t *frame_create_square(uint_t L, bool dual) {
-	uint_t nv = 2 * pow(L / 2, 2);
-	uint_t ne = pow(L, 2);
-
-	uint_t *ev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
-	uint_t *dev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
-	double *vx = (double *)malloc(2 * nv * sizeof(double));
-	double *dvx = (double *)malloc(2 * nv * sizeof(double));
-
-	for (uint_t i = 0; i < ne; i++) {
-		uint_t x = i / L;
-		uint_t y = i % L;
-
-		ev[2 * i] = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2);
-		ev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv;
-
-		dev[2 * i] = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2);
-		dev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv;
-	}
-
-	double dx = 1. / L;
-
-	for (uint_t i = 0; i < nv; i++) {
-		vx[2 * i] = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2)));
-		vx[2 * i + 1] = dx * (i / (L / 2));
-
-		dvx[2 * i] = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2)));
-		dvx[2 * i + 1] = dx * (i / (L / 2));
-	}
-
-	frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
-	frame->ne = ne;
-	frame->nv = nv;
-	frame->dnv = nv;
-
-	if (dual) {
-		frame->ev = dev;
-		frame->dev = ev;
-		frame->vx = dvx;
-		frame->dvx = vx;
-	} else {
-		frame->ev = ev;
-		frame->dev = dev;
-		frame->vx = vx;
-		frame->dvx = dvx;
-	}
-
-	return frame;
-}
-
-frame_t *frame_create(lattice_t lattice, uint_t L, bool dual) {
-	switch (lattice) {
-		case (SQUARE_LATTICE):
-			return frame_create_square(L, dual);
-		case (VORONOI_LATTICE):
-			return frame_create_voronoi(L, dual, false);
-		case (VORONOI_HYPERUNIFORM_LATTICE):
-			return frame_create_voronoi(L, dual, true);
-	}
-}
-
-void frame_free(frame_t *frame) {
-	free(frame->ev);
-	free(frame->dev);
-	free(frame->vx);
-	free(frame->dvx);
-	free(frame);
-}
-
diff --git a/src/get_dual_clusters.c b/src/get_dual_clusters.c
deleted file mode 100644
index 78bf185..0000000
--- a/src/get_dual_clusters.c
+++ /dev/null
@@ -1,36 +0,0 @@
-
-#include "fracture.h"
-
-unsigned int *get_clusters(net_t *instance, cholmod_common *c) {
-	cholmod_sparse *s_dual = gen_adjacency(instance, true, false, true, c);
-
-	unsigned int *dual_marks = find_components(s_dual, 0);
-	CHOL_F(free_sparse)(&s_dual, c);
-
-	return dual_marks;
-}
-
-unsigned int *get_cluster_dist(net_t *instance, cholmod_common *c) {
-	unsigned int *clusters = get_clusters(instance, c);
-	unsigned int *cluster_dist = (unsigned int *)calloc(
-			instance->graph->dnv, sizeof(unsigned int));
-
-	unsigned int cur_mark = 0;
-	while (true) {
-		cur_mark++;
-		unsigned int num_in_cluster = 0;
-		for (unsigned int i = 0; i < instance->graph->dnv; i++) {
-			if (clusters[i] == cur_mark)
-				num_in_cluster++;
-		}
-
-		if (num_in_cluster == 0)
-			break;
-
-		cluster_dist[num_in_cluster - 1]++;
-	}
-
-	free(clusters);
-
-	return cluster_dist;
-}
diff --git a/src/graph_components.c b/src/graph_components.c
deleted file mode 100644
index 9680675..0000000
--- a/src/graph_components.c
+++ /dev/null
@@ -1,185 +0,0 @@
-
-#include "fracture.h"
-
-bool set_connected(const cholmod_sparse *laplacian, uint_t *marks,
-									 int vertex, int label, int stop_at, int exclude) {
-	unsigned int num_verts = (int_t)laplacian->nrow;
-
-	int_t *ia = (int_t *)laplacian->p;
-	int_t *ja = (int_t *)laplacian->i;
-	double *a = (double *)laplacian->x;
-
-	bool stopped = false;
-
-	int *queue = (int *)malloc((num_verts) * sizeof(int));
-	assert(queue != NULL);
-
-	int queue_pos = 0;
-
-	int *component_list = (int *)malloc((num_verts) * sizeof(int));
-	assert(component_list != NULL);
-
-	int list_pos = 0;
-	queue[0] = vertex;
-	component_list[0] = vertex;
-
-	unsigned int old_label = marks[vertex];
-	marks[vertex] = label;
-
-	while (queue_pos >= 0) {
-		int next_vertex = queue[queue_pos];
-		queue_pos--;
-		for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) {
-			if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) {
-				// if we run into the other vertex, the components haven't changed
-				stopped = true; queue_pos = -1;
-				break;
-			}
-			if (marks[ja[ia[next_vertex] + j]] != label &&
-					a[ia[next_vertex] + j] != 0 &&
-					ja[ia[next_vertex] + j] < num_verts - exclude) {
-				marks[ja[ia[next_vertex] + j]] = label;
-				queue_pos++;
-				queue[queue_pos] = ja[ia[next_vertex] + j];
-				list_pos++;
-				component_list[list_pos] = ja[ia[next_vertex] + j];
-			}
-		}
-	}
-
-	if (stopped) {
-		for (unsigned int i = 0; i <= list_pos; i++) {
-			marks[component_list[i]] = old_label;
-		}
-	}
-
-	free(queue);
-	free(component_list);
-
-	if (stopped)
-		return false;
-	else
-		return true;
-}
-
-int update_components(const cholmod_sparse *laplacian, unsigned int *marks,
-											int old_num_components, int v1, int v2, int exclude) {
-	assert(laplacian != NULL);
-	assert(marks != NULL);
-
-	if (set_connected(laplacian, marks, v1, old_num_components + 1, v2,
-										exclude)) {
-		return old_num_components + 1;
-	} else
-		return old_num_components;
-}
-
-unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) {
-	size_t num_verts = laplacian->nrow;
-	unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int));
-
-	int num_components = 0;
-
-	for (int i = 0; i < num_verts; i++) {
-		if (marks[i] == 0) {
-			num_components++;
-			set_connected(laplacian, marks, i, num_components, -1, skip);
-		}
-	}
-
-	return marks;
-}
-
-uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const unsigned
-		int *vei, const uint_t *ve, int **cycles) {
-	uint_t num_cycles = 0;
-	uint_t *elist = (uint_t *)malloc(num_edges * sizeof(uint_t));
-	uint_t *side = (uint_t *)calloc(num_edges, sizeof(uint_t));
-	uint_t *num = (uint_t *)malloc(num_edges * sizeof(uint_t));
-	for (uint_t i = 0; i < num_edges; i++) {
-		if (fuses[i]) {
-		bool in_cycle = false;
-		for (uint_t j = 0; j < num_cycles; j++) {
-			uint_t k = 0;
-			int e;
-			while ((e = cycles[j][k]) >= 0) {
-				if (e == i) {
-					in_cycle = true;
-					break;
-				}
-				k++;
-			}
-			if (in_cycle) {
-				break;
-			}
-		}
-		if (!in_cycle) {
-			uint_t pos = 0;
-			uint_t cur_e = i;
-			bool *included = (bool *)calloc(num_edges, sizeof(bool));
-			included[cur_e] = true;
-			elist[0] = cur_e;
-			side[0] = 1;
-			num[0] = 0;
-			while (true) {
-				uint_t cv = ev[2 * cur_e + side[pos]];
-				uint_t new_e = cur_e;
-				uint_t j;
-				for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) {
-					new_e = ve[vei[cv] + j];
-					if (new_e != cur_e && fuses[new_e]) break;
-				}
-				if (new_e != cur_e && fuses[new_e]) {
-					if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) {
-						pos++;
-						cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int));
-						cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int));
-						for (uint_t i = 0; i < pos; i++) {
-							cycles[2*num_cycles][i] = elist[i];
-							cycles[2*num_cycles+1][i] = side[i];
-						}
-						cycles[2*num_cycles][pos] = -1;
-						cycles[2*num_cycles+1][pos] = -1;
-						num_cycles++;
-						pos--;
-						num[pos]++;
-						included[pos] = false;
-						cur_e = elist[pos];
-					} else if (included[new_e]) {
-						if (pos == 0) break;
-						included[cur_e] = false;
-						pos--;
-						num[pos]++;
-						cur_e = elist[pos];
-					} else {
-						num[pos] = j;
-						pos++;
-						elist[pos] = new_e;
-						included[new_e] = true;
-						uint_t nv1 = ev[2 * new_e];
-						uint_t nv2 = ev[2 * new_e + 1];
-						uint_t v = nv1 == cv ? nv2 : nv1;
-						side[pos] = v == nv1 ? 0 : 1;
-						num[pos] = 0;
-						cur_e = new_e;
-					}
-				} else {
-					if (pos == 0) break;
-					included[cur_e] = false;
-					pos--;
-					num[pos]++;
-					cur_e = elist[pos];
-				}
-			}
-			free(included);
-		}
-		}
-	}
-
-	free(elist);
-	free(side);
-	free(num);
-
-	return num_cycles;
-}
-
diff --git a/src/graph_create.c b/src/graph_create.c
deleted file mode 100644
index 09a3aed..0000000
--- a/src/graph_create.c
+++ /dev/null
@@ -1,299 +0,0 @@
-
-#include "fracture.h"
-
-uint_t *get_spanning_edges(uint_t num_edges, uint_t *edges_to_verts, double *vert_coords, double cut, uint_t *n) {
-	uint_t *spanning_edges = (uint_t *)malloc(num_edges * sizeof(uint_t));
-	(*n) = 0;
-	for (uint_t i = 0; i < num_edges; i++) {
-		uint_t v1, v2;
-		v1 = edges_to_verts[2 * i];
-		v2 = edges_to_verts[2 * i + 1];
-		double v1y, v2y;
-		v1y = vert_coords[2 * v1 + 1];
-		v2y = vert_coords[2 * v2 + 1];
-		if ((fabs(v1y - v2y) < 0.5) && ((v1y <= cut && v2y > cut) || (v1y >= cut && v2y < cut))) {
-			spanning_edges[*n] = i;
-			(*n)++;
-		}
-	}
-	return spanning_edges;
-}
-
-double *get_edge_coords(uint_t num_edges, double *vert_coords,
-												uint_t *edges_to_verts) {
-	double *output = (double *)malloc(2 * num_edges * sizeof(double));
-
-	for (uint_t i = 0; i < num_edges; i++) {
-		uint_t v1, v2;
-		double v1x, v1y, v2x, v2y, dx, dy;
-		v1 = edges_to_verts[2 * i];
-		v2 = edges_to_verts[2 * i + 1];
-		output[2 * i] = 0;
-		output[2 * i + 1] = 0;
-		v1x = vert_coords[2 * v1];
-		v1y = vert_coords[2 * v1 + 1];
-		v2x = vert_coords[2 * v2];
-		v2y = vert_coords[2 * v2 + 1];
-		dx = v1x - v2x;
-		dy = v1y - v2y;
-		if (fabs(dx) > 0.5) {
-			if (dx > 0) {
-				v2x = v2x + 1;
-			} else {
-				v1x = v1x + 1;
-			}
-		}
-		if (fabs(dy) > 0.5) {
-			if (dy > 0) {
-				v2y = v2y + 1;
-			} else {
-				v1y = v1y + 1;
-			}
-		}
-		output[2 * i] = (v1x + v2x) / 2;
-		output[2 * i + 1] = (v1y + v2y) / 2;
-	}
-
-	return output;
-}
-
-uint_t *get_verts_to_edges_ind(uint_t num_verts,
-																		 uint_t num_edges,
-																		 const uint_t *edges_to_verts) {
-	uint_t *output =
-			(uint_t *)calloc(num_verts + 1, sizeof(uint_t));
-	assert(output != NULL);
-
-	for (uint_t i = 0; i < 2 * num_edges; i++) {
-		if (edges_to_verts[i] < num_verts) {
-			output[edges_to_verts[i] + 1]++;
-		}
-	}
-
-	for (uint_t i = 0; i < num_verts; i++) {
-		output[i + 1] += output[i];
-	}
-
-	return output;
-}
-
-uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
-																 const uint_t *edges_to_verts,
-																 const uint_t *verts_to_edges_ind) {
-	uint_t *output = (uint_t *)calloc(verts_to_edges_ind[num_verts],
-																								sizeof(uint_t));
-	uint_t *counts =
-			(uint_t *)calloc(num_verts, sizeof(uint_t));
-	for (int i = 0; i < 2 * num_edges; i++) {
-		if (edges_to_verts[i] < num_verts) {
-			output[verts_to_edges_ind[edges_to_verts[i]] +
-						 counts[edges_to_verts[i]]] = i / 2;
-			counts[edges_to_verts[i]]++;
-		}
-	}
-
-	free(counts);
-
-	return output;
-}
-
-uint_t get_cut_edges(uint_t ne, const uint_t *ev, const double *vx, bool both, uint_t *ce) {
-	uint_t nce = 0;
-
-	for (uint_t i = 0; i < ne; i++) {
-		uint_t v1 = ev[2 * i];
-		uint_t v2 = ev[2 * i + 1];
-
-		double v1y = vx[2 * v1 + 1];
-		double v2y = vx[2 * v2 + 1];
-
-		if (fabs(v1y - v2y) > 0.5) {
-			ce[nce] = i;
-			nce++;
-		} else if (both) {
-			double v1x = vx[2 * v1];
-			double v2x = vx[2 * v2];
-
-			if (fabs(v1x - v2x) > 0.5) {
-				ce[nce] = i;
-				nce++;
-			}
-		}
-	}
-
-	return nce;
-}
-
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
-	graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
-	frame_t *f = frame_create(lattice, L, dual);
-
-	g->L = L;
-	g->boundary = bound;
-	g->ne = f->ne;
-
-	if (bound == TORUS_BOUND) {
-		g->nv = f->nv;
-		g->dnv = f->dnv;
-		g->nb = 1;
-
-		g->ev = f->ev;
-		f->ev = NULL;
-		g->dev = f->dev;
-		f->dev = NULL;
-		g->vx = f->vx;
-		f->vx = NULL;
-		g->dvx = f->dvx;
-		f->dvx = NULL;
-
-		// the boundary for the torus consists of a list of edges required to cut
-		// the torus into a cylinder
-		g->bi = (uint_t *)calloc(2, sizeof(uint_t));
-		g->b = (uint_t *)malloc(g->ne * sizeof(uint_t));
-		g->bi[1] = get_cut_edges(g->ne, g->ev, g->vx, false, g->b);
-		g->bq = (bool *)calloc(g->ne, sizeof(bool));
-		for (uint_t i = 0; i < g->bi[1]; i++) {
-			g->bq[g->b[i]] = true;
-		}
-	} else {
-		uint_t *ce = (uint_t *)malloc(g->ne * sizeof(uint_t));
-		uint_t nce = 0;
-
-		if (bound == CYLINDER_BOUND) {
-			g->nb = 2;
-			nce = get_cut_edges(f->ne, f->ev, f->vx, false, ce);
-		} else {
-			g->nb = 4;
-			nce = get_cut_edges(f->ne, f->ev, f->vx, true, ce);
-		}
-
-		g->nv = f->nv;
-		g->dnv = f->dnv;
-		g->vx = (double *)malloc(2 * (f->nv + nce) * sizeof(double));
-		g->dvx = (double *)malloc(2 * (f->dnv + nce) * sizeof(double));
-		g->ev = f->ev;
-		g->dev = f->dev;
-		f->ev = NULL;
-		f->dev = NULL;
-		memcpy(g->vx, f->vx, 2 * f->nv * sizeof(double));
-		memcpy(g->dvx, f->dvx, 2 * f->dnv * sizeof(double));
-
-		uint_t nbv = 0;
-		uint_t *bv = (uint_t *)calloc(f->nv, sizeof(uint_t));
-		uint_t *dbv = (uint_t *)calloc(f->dnv, sizeof(uint_t));
-		uint_t nside = 0;
-		bool *side = (bool *)calloc(f->nv, sizeof(bool));
-
-		for (uint_t i = 0; i < nce; i++) {
-			uint_t cv1 = g->ev[2 * ce[i]];
-			uint_t cv2 = g->ev[2 * ce[i] + 1];
-			uint_t dv1 = g->dev[2 * ce[i]];
-			uint_t dv2 = g->dev[2 * ce[i] + 1];
-
-			uint_t cin = 1;
-
-			if (bound == FREE_BOUND && (f->vx[2 * cv2] - f->vx[2 * cv1]) > 0.5) {
-				cin = 0; 
-			}
-
-			uint_t vin = f->vx[2 * cv1 + cin] < f->vx[2 * cv2 + cin] ? 0 : 1;
-			uint_t dvin = f->dvx[2 * dv1 + cin] < f->dvx[2 * dv2 + cin] ? 0 : 1;
-
-			if (bv[g->ev[2 * ce[i] + vin]] == 0) {
-				nbv++;
-				if (cin == 0) {
-					nside++;
-					side[g->ev[2 * ce[i] + vin]] = true;
-				}
-
-				bv[g->ev[2 * ce[i] + vin]] = g->nv;
-
-				g->vx[2 * g->nv + cin] = 1 + f->vx[2 * g->ev[2 * ce[i] + vin] + cin];
-				g->vx[2 * g->nv + !cin] = f->vx[2 * g->ev[2 * ce[i] + vin] + !cin];
-				g->ev[2 * ce[i] + vin] = g->nv;
-
-				g->nv++;
-			} else {
-				g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
-			}
-			if (dbv[g->dev[2 * ce[i] + dvin]] == 0) {
-				dbv[g->dev[2 * ce[i] + dvin]] = g->dnv;
-
-				if (f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin] < 0.5) {
-					g->dvx[2 * g->dnv + cin] = 1 + f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
-				} else {
-					g->dvx[2 * g->dnv + cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
-				}
-				g->dvx[2 * g->dnv + !cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + !cin];
-				g->dev[2 * ce[i] + dvin] = g->dnv;
-
-				g->dnv++;
-			} else {
-				g->dev[2 * ce[i] + dvin] = dbv[g->dev[2 * ce[i] + dvin]];
-			}
-		}
-
-		g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
-		g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));
-
-		g->bi[1] = ((int_t)nbv - (int_t)nside);
-		g->bi[g->nb] = 2 * nbv;
-
-		if (bound == FREE_BOUND) {
-			g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
-			g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
-		}
-
-		uint_t n_ins0 = 0;
-		uint_t n_ins1 = 0;
-
-		g->nbi = (uint_t *)malloc(((int_t)g->nv - (int_t)g->bi[g->nb]) * sizeof(uint_t));
-		g->bni = (uint_t *)malloc((g->nv + 1) * sizeof(uint_t));
-		g->bq = (bool *)calloc(g->nv, sizeof(bool));
-		uint_t n_tmp = 0;
-
-		for (uint_t i = 0; i < f->nv; i++) {
-			g->bni[i] = n_tmp;
-			if (bv[i] != 0) {
-				g->bq[i] = true;
-				g->bq[bv[i]] = true;
-				if (side[i]) {
-					g->b[g->bi[2] + n_ins1] = i;
-					g->b[g->bi[3] + n_ins1] = bv[i];
-					n_ins1++;
-				} else {
-					g->b[g->bi[0] + n_ins0] = i;
-					g->b[g->bi[1] + n_ins0] = bv[i];
-					n_ins0++;
-				}
-			} else {
-				g->nbi[n_tmp] = i;
-				n_tmp++;
-			}
-		}
-
-		for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
-			g->bni[f->nv + i] = n_tmp;
-		}
-
-		free(bv);
-		free(dbv);
-		free(side);
-	}
-
-	g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
-	g->ve = get_verts_to_edges(g->nv, g->ne, g->ev, g->vei);
-	g->dvei = get_verts_to_edges_ind(g->dnv, g->ne, g->dev);
-	g->dve = get_verts_to_edges(g->dnv, g->ne, g->dev, g->dvei);
-
-	g->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
-	uint_t num_spanning_edges;
-	g->spanning_edges = get_spanning_edges(g->ne, g->ev, g->vx, 0.5, &num_spanning_edges);
-	g->num_spanning_edges = num_spanning_edges;
-
-	frame_free(f);
-
-	return g;
-}
-
-
diff --git a/src/graph_free.c b/src/graph_free.c
deleted file mode 100644
index ee30c99..0000000
--- a/src/graph_free.c
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#include "fracture.h"
-
-void graph_free(graph_t *g, cholmod_common *c) {
-	free(g->ev);
-	free(g->vei);
-	free(g->ve);
-	free(g->bi);
-	free(g->b);
-	free(g->vx);
-	free(g->dvx);
-	free(g->dev);
-	free(g->dvei);
-	free(g->dve);
-	free(g->nbi);
-	free(g->bq);
-	CHOL_F(free_sparse)(&(g->voltcurmat), c);
-	free(g);
-}
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
deleted file mode 100644
index d396206..0000000
--- a/src/graph_genfunc.c
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#include "fracture.h"
-
-double *genfunc_uniform(uint_t num, gsl_rng *r) {
-	double *lattice = (double *)malloc(2 * num * sizeof(double));
-
-	for (uint_t i = 0; i < 2 * num; i++) {
-		lattice[i] = gsl_rng_uniform(r);
-	}
-
-	return lattice;
-}
-
-double *genfunc_hyperuniform(uint_t num, gsl_rng *r) {
-	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
deleted file mode 100644
index 716090d..0000000
--- a/src/spheres_poly/box.cpp
+++ /dev/null
@@ -1,1157 +0,0 @@
-/* 
-   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
deleted file mode 100644
index 69418f4..0000000
--- a/src/spheres_poly/box.h
+++ /dev/null
@@ -1,169 +0,0 @@
-/* 
-   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
deleted file mode 100644
index c8c104d..0000000
--- a/src/spheres_poly/event.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#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
deleted file mode 100644
index 55dd1fc..0000000
--- a/src/spheres_poly/event.h
+++ /dev/null
@@ -1,58 +0,0 @@
-#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
deleted file mode 100644
index 10430ed..0000000
--- a/src/spheres_poly/grid_field.h
+++ /dev/null
@@ -1,148 +0,0 @@
-/* 
-   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
deleted file mode 100644
index c83865c..0000000
--- a/src/spheres_poly/heap.cpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#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
deleted file mode 100644
index 91180e3..0000000
--- a/src/spheres_poly/heap.h
+++ /dev/null
@@ -1,42 +0,0 @@
-//---------------------------------------------------------------------------
-// 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
deleted file mode 100644
index 920c099..0000000
--- a/src/spheres_poly/neighbor.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#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
deleted file mode 100644
index 0fbbd68..0000000
--- a/src/spheres_poly/sphere.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#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
deleted file mode 100644
index 28c64b4..0000000
--- a/src/spheres_poly/sphere.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#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
deleted file mode 100644
index 9b14b01..0000000
--- a/src/spheres_poly/spheres.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-//===========================================================
-//===========================================================
-//===========================================================
-//
-//  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
deleted file mode 100644
index 44d8052..0000000
--- a/src/spheres_poly/spheres.h
+++ /dev/null
@@ -1,4 +0,0 @@
-
-#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
deleted file mode 100644
index 9ee481f..0000000
--- a/src/spheres_poly/vector.h
+++ /dev/null
@@ -1,471 +0,0 @@
-#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