ran clang-format on all files

23 files changed, 1879 insertions, 1849 deletions

diff --git a/lib/bound_set.c b/lib/bound_set.c
index 65f1e6f..32d0fb7 100644
--- a/lib/bound_set.c
+++ b/lib/bound_set.c
@@ -2,101 +2,109 @@
 #include "fracture.h"
 
 double th_p(double x, double y, double th) {
-	if (x >= 0 && y >= 0) return th;
-	else if (x <  0 && y >= 0) return M_PI - th;
-	else if (x <  0 && y <  0) return th - M_PI;
-	else return -th;
-
+  if (x >= 0 && y >= 0)
+    return th;
+  else if (x < 0 && y >= 0)
+    return M_PI - th;
+  else if (x < 0 && y < 0)
+    return th - M_PI;
+  else
+    return -th;
 }
 
 double u_y(double x, double y) {
-	double r = sqrt(pow(x, 2) + pow(y, 2));
-	double th = th_p(x, y, atan(fabs(y / x)));
+  double r = sqrt(pow(x, 2) + pow(y, 2));
+  double th = th_p(x, y, atan(fabs(y / x)));
 
-	return sqrt(r) * sin(th / 2);
+  return sqrt(r) * sin(th / 2);
 }
 
 void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
-	uint_t L = g->L;
-
-	for (uint_t i = 0; i < L / 2; i++) {
-		double x1, y1, x2, y2, x3, y3, x4, y4;
-		x1 = (2. * i + 1.) / L - notch_len; y1 = 0.5 - 1.;
-		x2 = (2. * i + 1.) / L - notch_len; y2 = 0.5 - 0.;
-		y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
-		y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
-
-		bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
-		bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
-		bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
-		bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
-	}
+  uint_t L = g->L;
+
+  for (uint_t i = 0; i < L / 2; i++) {
+    double x1, y1, x2, y2, x3, y3, x4, y4;
+    x1 = (2. * i + 1.) / L - notch_len;
+    y1 = 0.5 - 1.;
+    x2 = (2. * i + 1.) / L - notch_len;
+    y2 = 0.5 - 0.;
+    y3 = (2. * i + 1.) / L - 0.5;
+    x3 = 0.5 - 1.;
+    y4 = (2. * i + 1.) / L - 0.5;
+    x4 = 0.5 - 0.;
+
+    bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
+    bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
+    bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
+    bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
+  }
 }
 
 bool is_in(uint_t len, uint_t *list, uint_t element) {
-	for (uint_t i = 0; i < len; i++) {
-		if (list[i] == element) { 
-			return true;
-		}
-	}
-	return false;
+  for (uint_t i = 0; i < len; i++) {
+    if (list[i] == element) {
+      return true;
+    }
+  }
+  return false;
 }
 
-cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) {
-
-	uint_t dim = g->nv;
-
-	if (vb && g->boundary != TORUS_BOUND) {
-		dim -= g->bi[g->nb];
-	} else if (!vb) {
-		dim += 2;
-	}
-
-	cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
-	double *bound = (double *)boundary->x;
-
-	switch (g->boundary) {
-		case TORUS_BOUND:
-			for (uint_t i = 0; i < g->bi[1]; i++) {
-				uint_t be = g->b[i];
-				uint_t v1 = g->ev[2 * be];
-				uint_t v2 = g->ev[2 * be + 1];
-				double v1y = g->vx[2 * v1 + 1];
-				double v2y = g->vx[2 * v2 + 1];
-
-				uint_t ind = v1y < v2y ? 0 : 1;
-
-				bound[g->ev[2 * be + ind]] += 1;
-				bound[g->ev[2 * be + !ind]] -= 1;
-			}
-			break;
-			/*
-		case EMBEDDED_BOUND:
-			bound_set_embedded(bound, g, notch_len);
-			break;
-			*/
-		default:
-			if (vb) {
-				for (uint_t i = 0; i < dim; i++) {
-					uint_t v = g->nbi[i];
-					for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
-						uint_t e = g->ve[g->vei[v] + j];
-						uint_t v0 = g->ev[2 * e];
-						uint_t v1 = g->ev[2 * e + 1];
-
-						if (g->bq[v0] || g->bq[v1]) {
-							uint_t vv = v0 == v ? v1 : v0;
-							if (is_in(g->bi[1], g->b, vv)) {
-								bound[i]++;
-							}
-						}
-					}
-				}
-			} else {
-				bound[g->nv] = 1;
-				bound[g->nv + 1] = -1;
-			}
-	}
-
-	return boundary;
+cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len,
+                         cholmod_common *c) {
+
+  uint_t dim = g->nv;
+
+  if (vb && g->boundary != TORUS_BOUND) {
+    dim -= g->bi[g->nb];
+  } else if (!vb) {
+    dim += 2;
+  }
+
+  cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
+  double *bound = (double *)boundary->x;
+
+  switch (g->boundary) {
+  case TORUS_BOUND:
+    for (uint_t i = 0; i < g->bi[1]; i++) {
+      uint_t be = g->b[i];
+      uint_t v1 = g->ev[2 * be];
+      uint_t v2 = g->ev[2 * be + 1];
+      double v1y = g->vx[2 * v1 + 1];
+      double v2y = g->vx[2 * v2 + 1];
+
+      uint_t ind = v1y < v2y ? 0 : 1;
+
+      bound[g->ev[2 * be + ind]] += 1;
+      bound[g->ev[2 * be + !ind]] -= 1;
+    }
+    break;
+  /*
+case EMBEDDED_BOUND:
+  bound_set_embedded(bound, g, notch_len);
+  break;
+  */
+  default:
+    if (vb) {
+      for (uint_t i = 0; i < dim; i++) {
+        uint_t v = g->nbi[i];
+        for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) {
+          uint_t e = g->ve[g->vei[v] + j];
+          uint_t v0 = g->ev[2 * e];
+          uint_t v1 = g->ev[2 * e + 1];
+
+          if (g->bq[v0] || g->bq[v1]) {
+            uint_t vv = v0 == v ? v1 : v0;
+            if (is_in(g->bi[1], g->b, vv)) {
+              bound[i]++;
+            }
+          }
+        }
+      }
+    } else {
+      bound[g->nv] = 1;
+      bound[g->nv + 1] = -1;
+    }
+  }
+
+  return boundary;
 }
diff --git a/lib/break_edge.c b/lib/break_edge.c
index 2f112c2..e53f7c1 100644
--- a/lib/break_edge.c
+++ b/lib/break_edge.c
@@ -1,51 +1,50 @@
 
 #include "fracture.h"
 
-
 bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
-	net->fuses[edge] = true;
-	net->num_broken++;
-
-	double *x = (double *)net->boundary_cond->x;
-	const graph_t *g = net->graph;
-
-	uint_t v1 = net->graph->ev[2 * edge];
-	uint_t v2 = net->graph->ev[2 * edge + 1];
-
-	uint_t s1 = v1 < v2 ? v1 : v2;
-	uint_t s2 = v1 < v2 ? v2 : v1;
-
-	if (net->graph->boundary == TORUS_BOUND) {
-		if (net->graph->bq[edge]) {
-			double v1y = g->vx[2 * v1 + 1];
-			double v2y = g->vx[2 * v2 + 1];
-			uint_t ind = v1y < v2y ? 0 : 1;
-			x[g->ev[2 * edge + ind]] -= 1;
-			x[g->ev[2 * edge + !ind]] += 1;
-		}
-		if (net->factor != NULL) {
-			factor_update(net->factor, s1, s2, c);
-		}
-	} else if (net->voltage_bound) {
-		if (g->bq[v1] || g->bq[v2]) {
-			uint_t vv = g->bq[v1] ? v2 : v1;
-			uint_t uu = v1 == vv ? v2 : v1;
-			if (is_in(g->bi[1], g->b, uu)) {
-				x[g->bni[vv]] -= 1;
-			}
-			if (net->factor != NULL) {
-				factor_update2(net->factor, g->bni[vv], c);
-			}
-		} else {
-			if (net->factor != NULL) {
-				factor_update(net->factor, g->bni[s1], g->bni[s2], c);
-			}
-		}
-	} else {
-		if (net->factor != NULL) {
-			factor_update(net->factor, s1, s2, c);
-		}
-	}
-
-	return true;
+  net->fuses[edge] = true;
+  net->num_broken++;
+
+  double *x = (double *)net->boundary_cond->x;
+  const graph_t *g = net->graph;
+
+  uint_t v1 = net->graph->ev[2 * edge];
+  uint_t v2 = net->graph->ev[2 * edge + 1];
+
+  uint_t s1 = v1 < v2 ? v1 : v2;
+  uint_t s2 = v1 < v2 ? v2 : v1;
+
+  if (net->graph->boundary == TORUS_BOUND) {
+    if (net->graph->bq[edge]) {
+      double v1y = g->vx[2 * v1 + 1];
+      double v2y = g->vx[2 * v2 + 1];
+      uint_t ind = v1y < v2y ? 0 : 1;
+      x[g->ev[2 * edge + ind]] -= 1;
+      x[g->ev[2 * edge + !ind]] += 1;
+    }
+    if (net->factor != NULL) {
+      factor_update(net->factor, s1, s2, c);
+    }
+  } else if (net->voltage_bound) {
+    if (g->bq[v1] || g->bq[v2]) {
+      uint_t vv = g->bq[v1] ? v2 : v1;
+      uint_t uu = v1 == vv ? v2 : v1;
+      if (is_in(g->bi[1], g->b, uu)) {
+        x[g->bni[vv]] -= 1;
+      }
+      if (net->factor != NULL) {
+        factor_update2(net->factor, g->bni[vv], c);
+      }
+    } else {
+      if (net->factor != NULL) {
+        factor_update(net->factor, g->bni[s1], g->bni[s2], c);
+      }
+    }
+  } else {
+    if (net->factor != NULL) {
+      factor_update(net->factor, s1, s2, c);
+    }
+  }
+
+  return true;
 }
diff --git a/lib/correlations.c b/lib/correlations.c
index 63532ec..98c5767 100644
--- a/lib/correlations.c
+++ b/lib/correlations.c
@@ -2,46 +2,45 @@
 #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));
+  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 < 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];
-		}
-	}
+  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);
-	}
+  if (nulldists) {
+    for (int i = 0; i < nv; i++) {
+      free(dists[i]);
+    }
+    free(dists);
+  }
 
-	return corr;
+  return corr;
 }
-
diff --git a/lib/data.c b/lib/data.c
index 2047c44..e39c450 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -2,34 +2,34 @@
 #include "fracture.h"
 
 data_t *data_create(uint_t ne) {
-	data_t *data = malloc(1 * sizeof(data_t));
-	assert(data != NULL);
+  data_t *data = malloc(1 * sizeof(data_t));
+  assert(data != NULL);
 
-	data->num_broken = 0;
+  data->num_broken = 0;
 
-	data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
-	assert(data->break_list != NULL);
+  data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
+  assert(data->break_list != NULL);
 
-	data->extern_field = (long double *)malloc(ne * sizeof(long double));
-	assert(data->extern_field != NULL);
+  data->extern_field = (long double *)malloc(ne * sizeof(long double));
+  assert(data->extern_field != NULL);
 
-	data->conductivity = (double *)malloc(ne * sizeof(double));
-	assert(data->conductivity != NULL);
+  data->conductivity = (double *)malloc(ne * sizeof(double));
+  assert(data->conductivity != NULL);
 
-	return data;
+  return data;
 }
 
 void data_free(data_t *data) {
-	free(data->break_list);
-	free(data->extern_field);
-	free(data->conductivity);
-	free(data);
+  free(data->break_list);
+  free(data->extern_field);
+  free(data->conductivity);
+  free(data);
 }
 
-void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity) {
-	data->break_list[data->num_broken] = last_broke;
-	data->extern_field[data->num_broken] = strength;
-	data->conductivity[data->num_broken] = conductivity;
-	data->num_broken++;
+void data_update(data_t *data, uint_t last_broke, long double strength,
+                 double conductivity) {
+  data->break_list[data->num_broken] = last_broke;
+  data->extern_field[data->num_broken] = strength;
+  data->conductivity[data->num_broken] = conductivity;
+  data->num_broken++;
 }
-
diff --git a/lib/factor_update.c b/lib/factor_update.c
index a51705a..ceaa01a 100644
--- a/lib/factor_update.c
+++ b/lib/factor_update.c
@@ -1,68 +1,69 @@
 
 #include "fracture.h"
 
-void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c) {
-	uint_t n = factor->n;
+void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2,
+                   cholmod_common *c) {
+  uint_t n = factor->n;
 
-	cholmod_sparse *update_mat =
-			CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
+  cholmod_sparse *update_mat =
+      CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
 
-	uint_t s1, s2;
-	s1 = v1 < v2 ? v1 : v2;
-	s2 = v1 > v2 ? v1 : v2;
+  uint_t s1, s2;
+  s1 = v1 < v2 ? v1 : v2;
+  s2 = v1 > v2 ? v1 : v2;
 
-	int_t *pp = (int_t *)update_mat->p;
-	int_t *ii = (int_t *)update_mat->i;
-	double *xx = (double *)update_mat->x;
+  int_t *pp = (int_t *)update_mat->p;
+  int_t *ii = (int_t *)update_mat->i;
+  double *xx = (double *)update_mat->x;
 
-	for (uint_t i = 0; i <= s1; i++) {
-		pp[i] = 0;
-	}
+  for (uint_t i = 0; i <= s1; i++) {
+    pp[i] = 0;
+  }
 
-	for (uint_t i = s1 + 1; i <= n; i++) {
-		pp[i] = 2;
-	}
+  for (uint_t i = s1 + 1; i <= n; i++) {
+    pp[i] = 2;
+  }
 
-	ii[0] = s1;
-	ii[1] = s2;
-	xx[0] = 1;
-	xx[1] = -1;
+  ii[0] = s1;
+  ii[1] = s2;
+  xx[0] = 1;
+  xx[1] = -1;
 
-	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+      update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
 
-	CHOL_F(updown)(false, perm_update_mat, factor, c);
+  CHOL_F(updown)(false, perm_update_mat, factor, c);
 
-	CHOL_F(free_sparse)(&perm_update_mat, c);
-	CHOL_F(free_sparse)(&update_mat, c);
+  CHOL_F(free_sparse)(&perm_update_mat, c);
+  CHOL_F(free_sparse)(&update_mat, c);
 }
 
 void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
-	uint_t n = factor->n;
+  uint_t n = factor->n;
 
-	cholmod_sparse *update_mat =
-			CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
+  cholmod_sparse *update_mat =
+      CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
 
-	int_t *pp = (int_t *)update_mat->p;
-	int_t *ii = (int_t *)update_mat->i;
-	double *xx = (double *)update_mat->x;
+  int_t *pp = (int_t *)update_mat->p;
+  int_t *ii = (int_t *)update_mat->i;
+  double *xx = (double *)update_mat->x;
 
-	for (uint_t i = 0; i <= v; i++) {
-		pp[i] = 0;
-	}
+  for (uint_t i = 0; i <= v; i++) {
+    pp[i] = 0;
+  }
 
-	for (uint_t i = v + 1; i <= n; i++) {
-		pp[i] = 1;
-	}
+  for (uint_t i = v + 1; i <= n; i++) {
+    pp[i] = 1;
+  }
 
-	ii[0] = v;
-	xx[0] = 1;
+  ii[0] = v;
+  xx[0] = 1;
 
-	cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
-			update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
+  cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
+      update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
 
-	CHOL_F(updown)(false, perm_update_mat, factor, c);
+  CHOL_F(updown)(false, perm_update_mat, factor, c);
 
-	CHOL_F(free_sparse)(&perm_update_mat, c);
-	CHOL_F(free_sparse)(&update_mat, c);
+  CHOL_F(free_sparse)(&perm_update_mat, c);
+  CHOL_F(free_sparse)(&update_mat, c);
 }
diff --git a/lib/fracture.h b/lib/fracture.h
index b1114fb..7b4f146 100644
--- a/lib/fracture.h
+++ b/lib/fracture.h
@@ -17,7 +17,6 @@
 #include <stdbool.h>
 #include <stdint.h>
 #include <stdio.h>
-#include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <sys/types.h>
@@ -26,7 +25,6 @@
 #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
@@ -36,68 +34,69 @@
 #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;
+  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;
+  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);
+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_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 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);
+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);
+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_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_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);
+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);
+                   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);
+                               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);
+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);
 
@@ -112,11 +111,13 @@ 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);
+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);
+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);
 
diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c
index 1cc93a4..75dccce 100644
--- a/lib/gen_laplacian.c
+++ b/lib/gen_laplacian.c
@@ -1,137 +1,142 @@
 
 #include "fracture.h"
 
-cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) {
-	const graph_t *g = net->graph;
-	uint_t nv;
-	uint_t ne;
-	uint_t nre;
-	uint_t *ev;
-
-	if (use_gp) {
-		nv = net->dim;
-		ne = net->nep;
-		nre = (int_t)net->nep - (int_t)net->num_broken;
-		ev = net->evp;
-	} else if (dual) {
-		nv = g->dnv;
-		ne = g->ne;
-		nre = net->num_broken;
-		ev = g->dev;
-	} else {
-		nv = g->nv;
-		ne = g->ne;
-		nre = (int_t)g->ne - (int_t)net->num_broken;
-		ev = g->ev;
-	}
-
-	uint_t nnz = nre;
-
-	cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
-	int_t *ri = (int_t *)t->i;
-	int_t *ci = (int_t *)t->j;
-	double *ai = (double *)t->x;
-
-	t->nnz = nnz;
-
-	uint_t a = 0;
-
-	for (uint_t i = 0; i < ne; i++) {
-		if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
-			uint_t v1 = ev[2 * i];
-			uint_t v2 = ev[2 * i + 1];
-
-			uint_t s1 = v1 < v2 ? v1 : v2;
-			uint_t s2 = v1 < v2 ? v2 : v1;
-
-			ri[a] = s2;
-			ci[a] = s1;
-			ai[a] = 1;
-
-			a++;
-		}
-	}
-
-	cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
-	CHOL_F(free_triplet)(&t, c);
-
-	if (!symmetric) {
-		cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
-		CHOL_F(free_sparse)(&s, c);
-		s = tmp_s;
-	}
-
-	return s;
+cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp,
+                              bool symmetric, cholmod_common *c) {
+  const graph_t *g = net->graph;
+  uint_t nv;
+  uint_t ne;
+  uint_t nre;
+  uint_t *ev;
+
+  if (use_gp) {
+    nv = net->dim;
+    ne = net->nep;
+    nre = (int_t)net->nep - (int_t)net->num_broken;
+    ev = net->evp;
+  } else if (dual) {
+    nv = g->dnv;
+    ne = g->ne;
+    nre = net->num_broken;
+    ev = g->dev;
+  } else {
+    nv = g->nv;
+    ne = g->ne;
+    nre = (int_t)g->ne - (int_t)net->num_broken;
+    ev = g->ev;
+  }
+
+  uint_t nnz = nre;
+
+  cholmod_triplet *t =
+      CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
+
+  int_t *ri = (int_t *)t->i;
+  int_t *ci = (int_t *)t->j;
+  double *ai = (double *)t->x;
+
+  t->nnz = nnz;
+
+  uint_t a = 0;
+
+  for (uint_t i = 0; i < ne; i++) {
+    if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
+      uint_t v1 = ev[2 * i];
+      uint_t v2 = ev[2 * i + 1];
+
+      uint_t s1 = v1 < v2 ? v1 : v2;
+      uint_t s2 = v1 < v2 ? v2 : v1;
+
+      ri[a] = s2;
+      ci[a] = s1;
+      ai[a] = 1;
+
+      a++;
+    }
+  }
+
+  cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
+  CHOL_F(free_triplet)(&t, c);
+
+  if (!symmetric) {
+    cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
+    CHOL_F(free_sparse)(&s, c);
+    s = tmp_s;
+  }
+
+  return s;
 }
 
 cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
-	const graph_t *g = net->graph;
-	uint_t nv = net->dim;
-	uint_t ne = net->nep;
-	uint_t *ev = net->evp;
-
-	uint_t nnz = nv;
-
-	cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
-	int_t *rowind = (int_t *)temp_m->i;
-	int_t *colind = (int_t *)temp_m->j;
-	double *acoo = (double *)temp_m->x;
-
-	temp_m->nnz = nnz;
-
-	for (uint_t i = 0; i < nv; i++) {
-		rowind[i] = i;
-		colind[i] = i;
-		acoo[i] = 0;
-	}
-
-	cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
-
-	for (uint_t i = 0; i < ne; i++) {
-		if (!net->fuses[i]){
-			uint_t v1 = ev[2 * i];
-			uint_t v2 = ev[2 * i + 1];
-
-			acoo[v1]++;
-			acoo[v2]++;
-		}
-	}
-
-	if (net->voltage_bound && g->boundary != TORUS_BOUND) {
-		for (uint_t i = 0; i < net->dim; i++) {
-			uint_t v = g->nbi[i];
-			for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
-				uint_t e = g->ve[g->vei[v] + j];
-				uint_t v0 = g->ev[2 * e];
-				uint_t v1 = g->ev[2 * e + 1];
-
-				if (g->bq[v0] || g->bq[v1]) {
-					acoo[i]++;
-				}
-			}
-		}
-	} else {
-		acoo[0]++;
-	}
-
-	for (uint_t i = 0; i < nv; i++) {
-		if (acoo[i] == 0) acoo[i]++;
-	}
-
-	//assert(CHOL_F(check_triplet)(temp_m, c));
-
-	cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
-	//assert(CHOL_F(check_sparse)(t_out, c));
-
-	double alpha[2] = {1, 0};
-	double beta[2] = {-1, 0};
-	cholmod_sparse *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
-
-	CHOL_F(free_sparse)(&t_out, c);
-	CHOL_F(free_sparse)(&adjacency, c);
-	CHOL_F(free_triplet)(&temp_m, c);
-
-	return laplacian;
+  const graph_t *g = net->graph;
+  uint_t nv = net->dim;
+  uint_t ne = net->nep;
+  uint_t *ev = net->evp;
+
+  uint_t nnz = nv;
+
+  cholmod_triplet *temp_m =
+      CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
+
+  int_t *rowind = (int_t *)temp_m->i;
+  int_t *colind = (int_t *)temp_m->j;
+  double *acoo = (double *)temp_m->x;
+
+  temp_m->nnz = nnz;
+
+  for (uint_t i = 0; i < nv; i++) {
+    rowind[i] = i;
+    colind[i] = i;
+    acoo[i] = 0;
+  }
+
+  cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
+
+  for (uint_t i = 0; i < ne; i++) {
+    if (!net->fuses[i]) {
+      uint_t v1 = ev[2 * i];
+      uint_t v2 = ev[2 * i + 1];
+
+      acoo[v1]++;
+      acoo[v2]++;
+    }
+  }
+
+  if (net->voltage_bound && g->boundary != TORUS_BOUND) {
+    for (uint_t i = 0; i < net->dim; i++) {
+      uint_t v = g->nbi[i];
+      for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) {
+        uint_t e = g->ve[g->vei[v] + j];
+        uint_t v0 = g->ev[2 * e];
+        uint_t v1 = g->ev[2 * e + 1];
+
+        if (g->bq[v0] || g->bq[v1]) {
+          acoo[i]++;
+        }
+      }
+    }
+  } else {
+    acoo[0]++;
+  }
+
+  for (uint_t i = 0; i < nv; i++) {
+    if (acoo[i] == 0)
+      acoo[i]++;
+  }
+
+  // assert(CHOL_F(check_triplet)(temp_m, c));
+
+  cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
+  // assert(CHOL_F(check_sparse)(t_out, c));
+
+  double alpha[2] = {1, 0};
+  double beta[2] = {-1, 0};
+  cholmod_sparse *laplacian =
+      CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
+
+  CHOL_F(free_sparse)(&t_out, c);
+  CHOL_F(free_sparse)(&adjacency, c);
+  CHOL_F(free_triplet)(&temp_m, c);
+
+  return laplacian;
 }
diff --git a/lib/gen_voltcurmat.c b/lib/gen_voltcurmat.c
index e870140..fede836 100644
--- a/lib/gen_voltcurmat.c
+++ b/lib/gen_voltcurmat.c
@@ -2,35 +2,35 @@
 #include "fracture.h"
 
 cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
-															 unsigned int *edges_to_verts,
-															 cholmod_common *c) {
+                               unsigned int *edges_to_verts,
+                               cholmod_common *c) {
 
-	cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
-			num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
-	assert(t_m != NULL);
+  cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
+      num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
+  assert(t_m != NULL);
 
-	int_t *rowind = (int_t *)t_m->i;
-	int_t *colind = (int_t *)t_m->j;
-	double *acoo = (double *)t_m->x;
+  int_t *rowind = (int_t *)t_m->i;
+  int_t *colind = (int_t *)t_m->j;
+  double *acoo = (double *)t_m->x;
 
-	for (unsigned int i = 0; i < num_edges; i++) {
-		unsigned int v1 = edges_to_verts[2 * i];
-		unsigned int v2 = edges_to_verts[2 * i + 1];
-		rowind[2 * i] = i;
-		rowind[2 * i + 1] = i;
-		colind[2 * i] = v1;
-		colind[2 * i + 1] = v2;
-		acoo[2 * i] = 1;
-		acoo[2 * i + 1] = -1;
-	}
+  for (unsigned int i = 0; i < num_edges; i++) {
+    unsigned int v1 = edges_to_verts[2 * i];
+    unsigned int v2 = edges_to_verts[2 * i + 1];
+    rowind[2 * i] = i;
+    rowind[2 * i + 1] = i;
+    colind[2 * i] = v1;
+    colind[2 * i + 1] = v2;
+    acoo[2 * i] = 1;
+    acoo[2 * i + 1] = -1;
+  }
 
-	t_m->nnz = num_edges * 2;
+  t_m->nnz = num_edges * 2;
 
-	assert(CHOL_F(check_triplet)(t_m, c));
+  assert(CHOL_F(check_triplet)(t_m, c));
 
-	cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
+  cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
 
-	CHOL_F(free_triplet)(&t_m, c);
+  CHOL_F(free_triplet)(&t_m, c);
 
-	return m;
+  return m;
 }
diff --git a/lib/geometry.c b/lib/geometry.c
index ec788f1..2c1987b 100644
--- a/lib/geometry.c
+++ b/lib/geometry.c
@@ -2,54 +2,54 @@
 #include "fracture.h"
 
 int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
-									bool index) {
-	assert(edge < pow(width, 2));
-
-	int x = edge / width + 1;
-	int y = edge % width + 1;
-
-	if (periodic) {
-		return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
-								width +
-						(x - index) * width) /
-					 2;
-	} else {
-		return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
-						(x - index) * (width + 1) - 1) /
-					 2;
-	}
+                  bool index) {
+  assert(edge < pow(width, 2));
+
+  int x = edge / width + 1;
+  int y = edge % width + 1;
+
+  if (periodic) {
+    return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
+                width +
+            (x - index) * width) /
+           2;
+  } else {
+    return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
+            (x - index) * (width + 1) - 1) /
+           2;
+  }
 }
 
 int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
-											 bool index) {
-	assert(edge < pow(width, 2));
-
-	int x = edge / width + 1;
-	int y = edge % width + 1;
-
-	if (periodic) {
-		return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
-								width +
-						(x - index) * width) /
-					 2;
-	} else {
-		return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
-						(x - index) * (width + 1) - 1) /
-					 2;
-	}
+                       bool index) {
+  assert(edge < pow(width, 2));
+
+  int x = edge / width + 1;
+  int y = edge % width + 1;
+
+  if (periodic) {
+    return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
+                width +
+            (x - index) * width) /
+           2;
+  } else {
+    return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
+            (x - index) * (width + 1) - 1) /
+           2;
+  }
 }
 
 double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
-													bool index) {
-	if (periodic) {
-		if (index)
-			return (2 * vert) % width + (2 * vert / width) % 2;
-		else
-			return 2 * vert / width;
-	} else {
-		if (index)
-			return (2 * vert) % (width + 1);
-		else
-			return (2 * vert) / (width + 1);
-	}
+                          bool index) {
+  if (periodic) {
+    if (index)
+      return (2 * vert) % width + (2 * vert / width) % 2;
+    else
+      return 2 * vert / width;
+  } else {
+    if (index)
+      return (2 * vert) % (width + 1);
+    else
+      return (2 * vert) / (width + 1);
+  }
 }
diff --git a/lib/get_dual_clusters.c b/lib/get_dual_clusters.c
index eaf7562..9336106 100644
--- a/lib/get_dual_clusters.c
+++ b/lib/get_dual_clusters.c
@@ -2,29 +2,30 @@
 #include "fracture.h"
 
 components_t *get_clusters(net_t *instance) {
-	components_t *c = graph_components_get(instance->graph, instance->fuses, true);
-	return c;
+  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));
+  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++;
-		}
+  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;
+    if (num_in_cluster == 0)
+      break;
 
-		cluster_dist[num_in_cluster - 1]++;
-	}
+    cluster_dist[num_in_cluster - 1]++;
+  }
 
-	graph_components_free(c);
+  graph_components_free(c);
 
-	return cluster_dist;
+  return cluster_dist;
 }
diff --git a/lib/net.c b/lib/net.c
index d26b22c..ac89415 100644
--- a/lib/net.c
+++ b/lib/net.c
@@ -2,144 +2,145 @@
 #include "fracture.h"
 
 long double *get_thres(uint_t ne, double beta) {
-	long double *thres = (long double *)malloc(ne * sizeof(long double));
-	assert(thres != NULL);
+  long double *thres = (long double *)malloc(ne * sizeof(long double));
+  assert(thres != NULL);
 
-	gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
-	gsl_rng_set(r, jst_rand_seed());
+  gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
+  gsl_rng_set(r, jst_rand_seed());
 
-	for (uint_t i = 0; i < ne; i++) {
-		thres[i] = rand_dist_pow(r, beta);
-	}
+  for (uint_t i = 0; i < ne; i++) {
+    thres[i] = rand_dist_pow(r, beta);
+  }
 
-	gsl_rng_free(r);
+  gsl_rng_free(r);
 
-	return thres;
+  return thres;
 }
 
 void net_notch(net_t *net, double notch_len, cholmod_common *c) {
-	for (uint_t i = 0; i < net->graph->ne; i++) {
-		uint_t v1, v2;
-		double v1x, v1y, v2x, v2y, dy;
-		bool crosses_center, not_wrapping, correct_length;
-
-		v1 = net->graph->ev[2 * i];
-		v2 = net->graph->ev[2 * i + 1];
-
-		v1x = net->graph->vx[2 * v1];
-		v1y = net->graph->vx[2 * v1 + 1];
-		v2x = net->graph->vx[2 * v2];
-		v2y = net->graph->vx[2 * v2 + 1];
-
-		dy = v1y - v2y;
-
-		crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
-		not_wrapping = fabs(dy) < 0.5;
-		//correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
-		correct_length = v1x < notch_len && v2x < notch_len;
-
-		if (crosses_center && not_wrapping && correct_length) {
-			break_edge(net, i, c);
-		}
-	}
+  for (uint_t i = 0; i < net->graph->ne; i++) {
+    uint_t v1, v2;
+    double v1x, v1y, v2x, v2y, dy;
+    bool crosses_center, not_wrapping, correct_length;
+
+    v1 = net->graph->ev[2 * i];
+    v2 = net->graph->ev[2 * i + 1];
+
+    v1x = net->graph->vx[2 * v1];
+    v1y = net->graph->vx[2 * v1 + 1];
+    v2x = net->graph->vx[2 * v2];
+    v2y = net->graph->vx[2 * v2 + 1];
+
+    dy = v1y - v2y;
+
+    crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
+    not_wrapping = fabs(dy) < 0.5;
+    // correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
+    correct_length = v1x < notch_len && v2x < notch_len;
+
+    if (crosses_center && not_wrapping && correct_length) {
+      break_edge(net, i, c);
+    }
+  }
 }
 
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
-	net_t *net = (net_t *)calloc(1, sizeof(net_t));
-	assert(net != NULL);
-
-	net->graph = g;
-	net->num_broken = 0;
-
-	net->fuses = (bool *)calloc(g->ne, sizeof(bool));
-	assert(net->fuses != NULL);
-
-	net->thres = get_thres(g->ne, beta);
-	net->inf = inf;
-
-	net->dim = g->nv;
-
-	if (g->boundary == TORUS_BOUND) {
-		net->nep = g->ne;
-		net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
-		memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
-	} else {
-		if (vb) {
-			net->dim -= g->bi[g->nb];
-			net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
-			net->nep = 0;
-			for (uint_t i = 0; i < g->ne; i++) {
-				if (!(g->bq[g->ev[2*i]] || g->bq[g->ev[2*i+1]])) {
-					net->evp[2*net->nep] = g->bni[g->ev[2*i]];
-					net->evp[2*net->nep + 1] = g->bni[g->ev[2*i + 1]];
-					net->nep++;
-				}
-			}
-		} else {
-			net->dim += 2;
-			net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
-			memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
-			net->nep = g->ne + g->bi[2];
-
-			for (uint_t i = 0; i < 2; i++) {
-				for (uint_t j = 0; j < g->bi[i+1] - g->bi[i]; j++){ 
-					net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
-					net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
-				}
-			}
-		}
-	}
-
-	net->voltage_bound = vb;
-	net->boundary_cond = bound_set(g, vb, notch_len, c);
-
-	net_notch(net, notch_len, c);
-
-	{
-		cholmod_sparse *laplacian = gen_laplacian(net, c);
-		net->factor = CHOL_F(analyze)(laplacian, c);
-		CHOL_F(factorize)(laplacian, net->factor, c);
-		CHOL_F(free_sparse)(&laplacian, c);
-	}
-
-	net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
-
-	return net;
+net_t *net_create(const graph_t *g, double inf, double beta, double notch_len,
+                  bool vb, cholmod_common *c) {
+  net_t *net = (net_t *)calloc(1, sizeof(net_t));
+  assert(net != NULL);
+
+  net->graph = g;
+  net->num_broken = 0;
+
+  net->fuses = (bool *)calloc(g->ne, sizeof(bool));
+  assert(net->fuses != NULL);
+
+  net->thres = get_thres(g->ne, beta);
+  net->inf = inf;
+
+  net->dim = g->nv;
+
+  if (g->boundary == TORUS_BOUND) {
+    net->nep = g->ne;
+    net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+    memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+  } else {
+    if (vb) {
+      net->dim -= g->bi[g->nb];
+      net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
+      net->nep = 0;
+      for (uint_t i = 0; i < g->ne; i++) {
+        if (!(g->bq[g->ev[2 * i]] || g->bq[g->ev[2 * i + 1]])) {
+          net->evp[2 * net->nep] = g->bni[g->ev[2 * i]];
+          net->evp[2 * net->nep + 1] = g->bni[g->ev[2 * i + 1]];
+          net->nep++;
+        }
+      }
+    } else {
+      net->dim += 2;
+      net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
+      memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
+      net->nep = g->ne + g->bi[2];
+
+      for (uint_t i = 0; i < 2; i++) {
+        for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
+          net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
+          net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
+        }
+      }
+    }
+  }
+
+  net->voltage_bound = vb;
+  net->boundary_cond = bound_set(g, vb, notch_len, c);
+
+  net_notch(net, notch_len, c);
+
+  {
+    cholmod_sparse *laplacian = gen_laplacian(net, c);
+    net->factor = CHOL_F(analyze)(laplacian, c);
+    CHOL_F(factorize)(laplacian, net->factor, c);
+    CHOL_F(free_sparse)(&laplacian, c);
+  }
+
+  net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
+
+  return net;
 }
 
 net_t *net_copy(const net_t *net, cholmod_common *c) {
-	net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
-	assert(net_copy != NULL);
-	memcpy(net_copy, net, sizeof(net_t));
-
-	size_t fuses_size = (net->graph)->ne * sizeof(bool);
-	net_copy->fuses = (bool *)malloc(fuses_size);
-	assert(net_copy->fuses != NULL);
-	memcpy(net_copy->fuses, net->fuses, fuses_size);
-
-	size_t thres_size = (net->graph)->ne * sizeof(long double);
-	net_copy->thres = (long double *)malloc(thres_size);
-	assert(net_copy->thres != NULL);
-	memcpy(net_copy->thres, net->thres, thres_size);
-
-	size_t evp_size = 2 * net->nep * sizeof(uint_t);
-	net_copy->evp = (uint_t *)malloc(thres_size);
-	assert(net_copy->evp != NULL);
-	memcpy(net_copy->evp, net->evp, evp_size);
-
-	net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
-	net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
-	net_copy->voltcurmat = CHOL_F(copy_sparse)(net->voltcurmat, c);
-
-	return net_copy;
+  net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
+  assert(net_copy != NULL);
+  memcpy(net_copy, net, sizeof(net_t));
+
+  size_t fuses_size = (net->graph)->ne * sizeof(bool);
+  net_copy->fuses = (bool *)malloc(fuses_size);
+  assert(net_copy->fuses != NULL);
+  memcpy(net_copy->fuses, net->fuses, fuses_size);
+
+  size_t thres_size = (net->graph)->ne * sizeof(long double);
+  net_copy->thres = (long double *)malloc(thres_size);
+  assert(net_copy->thres != NULL);
+  memcpy(net_copy->thres, net->thres, thres_size);
+
+  size_t evp_size = 2 * net->nep * sizeof(uint_t);
+  net_copy->evp = (uint_t *)malloc(thres_size);
+  assert(net_copy->evp != NULL);
+  memcpy(net_copy->evp, net->evp, evp_size);
+
+  net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
+  net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
+  net_copy->voltcurmat = CHOL_F(copy_sparse)(net->voltcurmat, c);
+
+  return net_copy;
 }
 
 void net_free(net_t *net, cholmod_common *c) {
-	free(net->fuses);
-	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);
+  free(net->fuses);
+  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/lib/net_conductivity.c b/lib/net_conductivity.c
index e9325bb..61148da 100644
--- a/lib/net_conductivity.c
+++ b/lib/net_conductivity.c
@@ -2,34 +2,34 @@
 #include "fracture.h"
 
 double net_conductivity(const net_t *net, const double *voltages) {
-	if (net->voltage_bound) {
-		// the voltage drop across the network is fixed to one with voltage
-		// boundary conditions, so the conductivity is the total current flowing
-		double tot_cur = 0;
-		for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) {
-			uint_t e = net->graph->spanning_edges[i];
-
-			if (!net->fuses[e]) {
-				uint_t v1, v2, s1, s2;
-				double v1y, v2y;
-
-				v1 = net->graph->ev[2 * e];
-				v2 = net->graph->ev[2 * e + 1];
-
-				v1y = net->graph->vx[2 * v1 + 1];
-				v2y = net->graph->vx[2 * v2 + 1];
-
-				s1 = v1y < v2y ? v1 : v2;
-				s2 = v1y < v2y ? v2 : v1;
-
-				tot_cur += voltages[s1] - voltages[s2];
-			}
-		}
-
-		return fabs(tot_cur);
-	} else {
-		// the current across the network is fixed to one with current boundary
-		// conditions, so the conductivity is the inverse of the total voltage drop
-		return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]);
-	}
+  if (net->voltage_bound) {
+    // the voltage drop across the network is fixed to one with voltage
+    // boundary conditions, so the conductivity is the total current flowing
+    double tot_cur = 0;
+    for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) {
+      uint_t e = net->graph->spanning_edges[i];
+
+      if (!net->fuses[e]) {
+        uint_t v1, v2, s1, s2;
+        double v1y, v2y;
+
+        v1 = net->graph->ev[2 * e];
+        v2 = net->graph->ev[2 * e + 1];
+
+        v1y = net->graph->vx[2 * v1 + 1];
+        v2y = net->graph->vx[2 * v2 + 1];
+
+        s1 = v1y < v2y ? v1 : v2;
+        s2 = v1y < v2y ? v2 : v1;
+
+        tot_cur += voltages[s1] - voltages[s2];
+      }
+    }
+
+    return fabs(tot_cur);
+  } else {
+    // the current across the network is fixed to one with current boundary
+    // conditions, so the conductivity is the inverse of the total voltage drop
+    return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]);
+  }
 }
diff --git a/lib/net_currents.c b/lib/net_currents.c
index a078336..32dba04 100644
--- a/lib/net_currents.c
+++ b/lib/net_currents.c
@@ -1,51 +1,52 @@
 
 #include "fracture.h"
 
-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->voltcurmat;
-
-	cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
-
-	double *tmp_x = x->x;
-	x->x = (void *)voltages;
-
-	cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);
-
-	double alpha[2] = {1, 0};
-	double beta[2] = {0, 0};
-	CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
-	double *currents = (double *)malloc(ne * sizeof(double));
-
-	for (int i = 0; i < ne; i++) {
-		if (net->fuses[i]) {
-			currents[i] = 0;
-		} else {
-			currents[i] = ((double *)y->x)[i];
-		}
-	}
-
-	if (net->graph->boundary == TORUS_BOUND) {
-		for (uint_t i = 0; i < net->graph->bi[1]; i++) {
-			uint_t e = net->graph->b[i];
-			uint_t v1 = net->graph->ev[2 * e];
-			uint_t v2 = net->graph->ev[2 * e + 1];
-			double v1y = net->graph->vx[2 * v1 + 1];
-			double v2y = net->graph->vx[2 * v2 + 1];
-
-			if (v1y > v2y) {
-				currents[e] += 1;
-			} else {
-				currents[e] -= 1;
-			}
-		}
-	}
-
-	x->x = tmp_x;
-	CHOL_F(free_dense)(&x, c);
-	CHOL_F(free_dense)(&y, c);
-
-	return currents;
+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->voltcurmat;
+
+  cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
+
+  double *tmp_x = x->x;
+  x->x = (void *)voltages;
+
+  cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);
+
+  double alpha[2] = {1, 0};
+  double beta[2] = {0, 0};
+  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+  double *currents = (double *)malloc(ne * sizeof(double));
+
+  for (int i = 0; i < ne; i++) {
+    if (net->fuses[i]) {
+      currents[i] = 0;
+    } else {
+      currents[i] = ((double *)y->x)[i];
+    }
+  }
+
+  if (net->graph->boundary == TORUS_BOUND) {
+    for (uint_t i = 0; i < net->graph->bi[1]; i++) {
+      uint_t e = net->graph->b[i];
+      uint_t v1 = net->graph->ev[2 * e];
+      uint_t v2 = net->graph->ev[2 * e + 1];
+      double v1y = net->graph->vx[2 * v1 + 1];
+      double v2y = net->graph->vx[2 * v2 + 1];
+
+      if (v1y > v2y) {
+        currents[e] += 1;
+      } else {
+        currents[e] -= 1;
+      }
+    }
+  }
+
+  x->x = tmp_x;
+  CHOL_F(free_dense)(&x, c);
+  CHOL_F(free_dense)(&y, c);
+
+  return currents;
 }
diff --git a/lib/net_fracture.c b/lib/net_fracture.c
index e7f18fc..65ede9b 100644
--- a/lib/net_fracture.c
+++ b/lib/net_fracture.c
@@ -2,66 +2,67 @@
 #include "fracture.h"
 
 uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
-	uint_t max_pos = UINT_MAX;
-	long double max_val = 0;
+  uint_t max_pos = UINT_MAX;
+  long double max_val = 0;
 
-	for (uint_t i = 0; i < net->graph->ne; i++) {
-		long double current = fabs(currents[i]);
-		bool broken = net->fuses[i];
+  for (uint_t i = 0; i < net->graph->ne; i++) {
+    long double current = fabs(currents[i]);
+    bool broken = net->fuses[i];
 
-		if (!broken && current > cutoff) {
-			long double val = current / net->thres[i];
+    if (!broken && current > cutoff) {
+      long double val = current / net->thres[i];
 
-			if (val > max_val) {
-				max_val = val;
-				max_pos = i;
-			}
-		}
-	}
+      if (val > max_val) {
+        max_val = val;
+        max_pos = i;
+      }
+    }
+  }
 
-	if (max_pos == UINT_MAX) {
-		printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
-		exit(EXIT_FAILURE);
-	}
+  if (max_pos == UINT_MAX) {
+    printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
+    exit(EXIT_FAILURE);
+  }
 
-	return max_pos;
+  return max_pos;
 }
 
 data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
-	data_t *data = data_create(net->graph->ne);
+  data_t *data = data_create(net->graph->ne);
 
-	uint_t n = 0;
-	while (true) {
-		n++;
-		double *voltages = net_voltages(net, c);
-		double *currents = net_currents(net, voltages, c);
+  uint_t n = 0;
+  while (true) {
+    n++;
+    double *voltages = net_voltages(net, c);
+    double *currents = net_currents(net, voltages, c);
 
-		double conductivity = net_conductivity(net, voltages);
+    double conductivity = net_conductivity(net, voltages);
 
-		if (conductivity < cutoff) {
-			free(voltages);
-			free(currents);
-			break;
-		}
+    if (conductivity < cutoff) {
+      free(voltages);
+      free(currents);
+      break;
+    }
 
-		uint_t last_broke = get_next_broken(net, currents, cutoff);
+    uint_t last_broke = get_next_broken(net, currents, cutoff);
 
-		long double sim_current;
+    long double sim_current;
 
-		if (net->voltage_bound) {
-			sim_current = conductivity;
-		} else {
-			sim_current = 1;
-		}
+    if (net->voltage_bound) {
+      sim_current = conductivity;
+    } else {
+      sim_current = 1;
+    }
 
-		data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity);
+    data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] /
+                                        ((long double)currents[last_broke])),
+                conductivity);
 
-		free(voltages);
-		free(currents);
+    free(voltages);
+    free(currents);
 
-		break_edge(net, last_broke, c);
-	}
+    break_edge(net, last_broke, c);
+  }
 
-	return data;
+  return data;
 }
-
diff --git a/lib/net_voltages.c b/lib/net_voltages.c
index 7b07201..d456a65 100644
--- a/lib/net_voltages.c
+++ b/lib/net_voltages.c
@@ -2,39 +2,38 @@
 #include "fracture.h"
 
 double *net_voltages(const net_t *net, cholmod_common *c) {
-	cholmod_dense *b = net->boundary_cond;
-	cholmod_factor *factor = net->factor;
-
-	cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
-	if (((double *)x->x)[0] != ((double *)x->x)[0]) {
-		printf("GET_VOLTAGE: value is NaN\n");
-		exit(EXIT_FAILURE);
-	}
-
-	double *t_voltages = (double *)x->x;
-	x->x = NULL;
-	CHOL_F(free_dense)(&x, c);
-
-	const graph_t *g = net->graph;
-
-	if (g->boundary == TORUS_BOUND) {
-		return t_voltages;
-	} else if (net->voltage_bound) {
-		double *voltages = (double *)malloc(g->nv * sizeof(double));
-		for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
-			voltages[g->nbi[i]] = t_voltages[i];
-		}
-		for (uint_t i = 0; i < 2; i++) {
-			for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
-				voltages[g->b[g->bi[i] + j]] = 1 - i;
-			}
-		}
-
-		free(t_voltages);
-		return voltages;
-	} else {
-		return t_voltages;
-	}
+  cholmod_dense *b = net->boundary_cond;
+  cholmod_factor *factor = net->factor;
+
+  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+
+  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+    printf("GET_VOLTAGE: value is NaN\n");
+    exit(EXIT_FAILURE);
+  }
+
+  double *t_voltages = (double *)x->x;
+  x->x = NULL;
+  CHOL_F(free_dense)(&x, c);
+
+  const graph_t *g = net->graph;
+
+  if (g->boundary == TORUS_BOUND) {
+    return t_voltages;
+  } else if (net->voltage_bound) {
+    double *voltages = (double *)malloc(g->nv * sizeof(double));
+    for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
+      voltages[g->nbi[i]] = t_voltages[i];
+    }
+    for (uint_t i = 0; i < 2; i++) {
+      for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
+        voltages[g->b[g->bi[i] + j]] = 1 - i;
+      }
+    }
+
+    free(t_voltages);
+    return voltages;
+  } else {
+    return t_voltages;
+  }
 }
-
diff --git a/lib/rand.c b/lib/rand.c
index 16aeba4..8ba2b2b 100644
--- a/lib/rand.c
+++ b/lib/rand.c
@@ -2,14 +2,14 @@
 #include "fracture.h"
 
 long double rand_dist_pow(const gsl_rng *r, double beta) {
-	long double x = 0;
+  long double x = 0;
 
-	// underflow means that for very small beta x is sometimes identically zero,
-	// which causes problems
-	while (x == 0.0) {
-		long double y = logl(gsl_rng_uniform_pos(r)) / beta;
-		x = expl(y);
-	}
+  // underflow means that for very small beta x is sometimes identically zero,
+  // which causes problems
+  while (x == 0.0) {
+    long double y = logl(gsl_rng_uniform_pos(r)) / beta;
+    x = expl(y);
+  }
 
-	return x;
+  return x;
 }
diff --git a/src/anal_process.c b/src/anal_process.c
index 0ebec84..de27571 100644
--- a/src/anal_process.c
+++ b/src/anal_process.c
@@ -1,134 +1,135 @@
 
 #include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.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>
 
 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;
+  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);
+  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; 
-	}
+    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);
+  double momf = mom / total;
+  double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
 
-	result[0] = momf;
-	result[1] = momf_err;
+  result[0] = momf;
+  result[1] = momf_err;
 }
 
 int main(int argc, char *argv[]) {
-	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_c1 = (double *)malloc(nc * sizeof(double));
-	double *errors_c1 = (double *)malloc(nc * sizeof(double));
-	double *vals_c2 = (double *)malloc(nc * sizeof(double));
-	double *errors_c2 = (double *)malloc(nc * sizeof(double));
-	double *vals_c3 = (double *)malloc(nc * sizeof(double));
-	double *errors_c3 = (double *)malloc(nc * sizeof(double));
-	double *vals_c4 = (double *)malloc(nc * sizeof(double));
-	double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
-	char *out_filename = (char *)malloc(100 * sizeof(char));
-	uint_t out_filename_len = 0;
-
-	for (uint_t i = 0; i < nc; i++) {
-		char *fn = argv[1 + i];
-		char *buff = (char *)malloc(20 * sizeof(char));
-		uint_t pos = 0; uint_t c = 0;
-		while (c < 5) {
-			if (fn[pos] == '_') {
-				c++;
-			}
-			if (i == 0) {
-				out_filename[pos] = fn[pos];
-				out_filename_len++;
-			}
-			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 mom1[2];
-			mom(dist_len, 1, dist, mom1);
-			vals_c1[i] = mom1[0];
-			errors_c1[i] = mom1[1];
-
-			double mom2[2];
-			mom(dist_len, 2, dist, mom2);
-			vals_c2[i] = mom2[0];
-			errors_c2[i] = mom2[1];
-
-			double mom3[2];
-			mom(dist_len, 3, dist, mom3);
-			vals_c3[i] = mom3[0];
-			errors_c3[i] = mom3[1];
-
-			double mom4[2];
-			mom(dist_len, 4, dist, mom4);
-			vals_c4[i] = mom4[0];
-			errors_c4[i] = mom4[1];
-		}
-		free(dist);
-	}
-
-	out_filename[out_filename_len-1] = '.';
-	out_filename[out_filename_len] = 't';
-	out_filename[out_filename_len+1] = 'x';
-	out_filename[out_filename_len+2] = 't';
-	out_filename[out_filename_len+3] = '\0';
-
-	FILE *out_file = fopen(out_filename, "w");
-
-	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
-	}
-
-	fclose(out_file);
-
-
-	free(Ls_c);
-	free(betas_c);
-	free(vals_c1);
-	free(errors_c1);
-	free(vals_c2);
-	free(errors_c2);
-	free(vals_c3);
-	free(errors_c3);
-
-	return 0;
+  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_c1 = (double *)malloc(nc * sizeof(double));
+  double *errors_c1 = (double *)malloc(nc * sizeof(double));
+  double *vals_c2 = (double *)malloc(nc * sizeof(double));
+  double *errors_c2 = (double *)malloc(nc * sizeof(double));
+  double *vals_c3 = (double *)malloc(nc * sizeof(double));
+  double *errors_c3 = (double *)malloc(nc * sizeof(double));
+  double *vals_c4 = (double *)malloc(nc * sizeof(double));
+  double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+  char *out_filename = (char *)malloc(100 * sizeof(char));
+  uint_t out_filename_len = 0;
+
+  for (uint_t i = 0; i < nc; i++) {
+    char *fn = argv[1 + i];
+    char *buff = (char *)malloc(20 * sizeof(char));
+    uint_t pos = 0;
+    uint_t c = 0;
+    while (c < 5) {
+      if (fn[pos] == '_') {
+        c++;
+      }
+      if (i == 0) {
+        out_filename[pos] = fn[pos];
+        out_filename_len++;
+      }
+      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 mom1[2];
+      mom(dist_len, 1, dist, mom1);
+      vals_c1[i] = mom1[0];
+      errors_c1[i] = mom1[1];
+
+      double mom2[2];
+      mom(dist_len, 2, dist, mom2);
+      vals_c2[i] = mom2[0];
+      errors_c2[i] = mom2[1];
+
+      double mom3[2];
+      mom(dist_len, 3, dist, mom3);
+      vals_c3[i] = mom3[0];
+      errors_c3[i] = mom3[1];
+
+      double mom4[2];
+      mom(dist_len, 4, dist, mom4);
+      vals_c4[i] = mom4[0];
+      errors_c4[i] = mom4[1];
+    }
+    free(dist);
+  }
+
+  out_filename[out_filename_len - 1] = '.';
+  out_filename[out_filename_len] = 't';
+  out_filename[out_filename_len + 1] = 'x';
+  out_filename[out_filename_len + 2] = 't';
+  out_filename[out_filename_len + 3] = '\0';
+
+  FILE *out_file = fopen(out_filename, "w");
+
+  for (uint_t i = 0; i < nc; i++) {
+    fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+            vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+            errors_c3[i], vals_c4[i], errors_c4[i]);
+  }
+
+  fclose(out_file);
+
+  free(Ls_c);
+  free(betas_c);
+  free(vals_c1);
+  free(errors_c1);
+  free(vals_c2);
+  free(errors_c2);
+  free(vals_c3);
+  free(errors_c3);
+
+  return 0;
 }
-
diff --git a/src/big_anal_process.c b/src/big_anal_process.c
index 0f7106f..8c1d8ba 100644
--- a/src/big_anal_process.c
+++ b/src/big_anal_process.c
@@ -1,156 +1,158 @@
 
 #include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.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 <sys/stat.h>
 
 void get_mean(uint_t len, double *data, double result[2]) {
-	double total = 0;
+  double total = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total += data[i];
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total += data[i];
+  }
 
-	double mean = total / len;
-	double total_sq_diff = 0;
+  double mean = total / len;
+  double total_sq_diff = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total_sq_diff += pow(mean - data[i], 2);
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total_sq_diff += pow(mean - data[i], 2);
+  }
 
-	double mean_err = sqrt(total_sq_diff) / len;
+  double mean_err = sqrt(total_sq_diff) / len;
 
-	result[0] = mean;
-	result[1] = mean_err;
+  result[0] = mean;
+  result[1] = mean_err;
 }
 
-void get_mom(uint_t len, uint_t n, double *data, double mean[2], double result[2]) {
-	double total_n_diff = 0;
-	double total_n2_diff = 0;
+void get_mom(uint_t len, uint_t n, double *data, double mean[2],
+             double result[2]) {
+  double total_n_diff = 0;
+  double total_n2_diff = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total_n_diff += pow(fabs(mean[0] - data[i]), n);
-		total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n);
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total_n_diff += pow(fabs(mean[0] - data[i]), n);
+    total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n);
+  }
 
-	double mom = total_n_diff / len;
-	double mom_err = sqrt(total_n2_diff) / len;
+  double mom = total_n_diff / len;
+  double mom_err = sqrt(total_n2_diff) / len;
 
-	result[0] = mom;
-	result[1] = mom_err;
+  result[0] = mom;
+  result[1] = mom_err;
 }
 
 int main(int argc, char *argv[]) {
-	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_c1 = (double *)malloc(nc * sizeof(double));
-	double *errors_c1 = (double *)malloc(nc * sizeof(double));
-	double *vals_c2 = (double *)malloc(nc * sizeof(double));
-	double *errors_c2 = (double *)malloc(nc * sizeof(double));
-	double *vals_c3 = (double *)malloc(nc * sizeof(double));
-	double *errors_c3 = (double *)malloc(nc * sizeof(double));
-	double *vals_c4 = (double *)malloc(nc * sizeof(double));
-	double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
-	char *out_filename = (char *)malloc(100 * sizeof(char));
-	uint_t out_filename_len = 0;
-
-	for (uint_t i = 0; i < nc; i++) {
-		char *fn = argv[1 + i];
-		char *buff = (char *)malloc(20 * sizeof(char));
-		uint_t pos = 0; uint_t c = 0;
-		while (c < 5) {
-			if (fn[pos] == '_') {
-				c++;
-			}
-			if (i == 0) {
-				out_filename[pos] = fn[pos];
-				out_filename_len++;
-			}
-			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);
-
-		struct stat info;
-		stat(fn, &info);
-		uint_t num_bytes = info.st_size;
-		uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(double);
-
-		double *dist = malloc(dist_len * sizeof(double));
-		FILE *dist_file = fopen(fn, "rb");
-		fread(dist, sizeof(double), dist_len, dist_file);
-		fclose(dist_file);
-		{
-			double mom1[2];
-			get_mean(dist_len, dist, mom1);
-			vals_c1[i] = mom1[0];
-			errors_c1[i] = mom1[1];
-
-			double mom2[2];
-			get_mom(dist_len, 2, dist, mom1, mom2);
-			vals_c2[i] = mom2[0];
-			errors_c2[i] = mom2[1];
-
-			double mom3[2];
-			get_mom(dist_len, 3, dist, mom1, mom3);
-			vals_c3[i] = mom3[0];
-			errors_c3[i] = mom3[1];
-
-			double mom4[2];
-			get_mom(dist_len, 4, dist, mom1, mom4);
-			vals_c4[i] = mom4[0];
-			errors_c4[i] = mom4[1];
-		}
-		free(dist);
-	}
-
-	out_filename[out_filename_len-1] = '.';
-	out_filename[out_filename_len] = 't';
-	out_filename[out_filename_len+1] = 'x';
-	out_filename[out_filename_len+2] = 't';
-	out_filename[out_filename_len+3] = '\0';
-
-	FILE *out_file = fopen(out_filename, "w");
-
-	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
-	}
-
-	fclose(out_file);
-
-
-	free(Ls_c);
-	free(betas_c);
-	free(vals_c1);
-	free(errors_c1);
-	free(vals_c2);
-	free(errors_c2);
-	free(vals_c3);
-	free(errors_c3);
-	free(vals_c4);
-	free(errors_c4);
-
-	return 0;
+  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_c1 = (double *)malloc(nc * sizeof(double));
+  double *errors_c1 = (double *)malloc(nc * sizeof(double));
+  double *vals_c2 = (double *)malloc(nc * sizeof(double));
+  double *errors_c2 = (double *)malloc(nc * sizeof(double));
+  double *vals_c3 = (double *)malloc(nc * sizeof(double));
+  double *errors_c3 = (double *)malloc(nc * sizeof(double));
+  double *vals_c4 = (double *)malloc(nc * sizeof(double));
+  double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+  char *out_filename = (char *)malloc(100 * sizeof(char));
+  uint_t out_filename_len = 0;
+
+  for (uint_t i = 0; i < nc; i++) {
+    char *fn = argv[1 + i];
+    char *buff = (char *)malloc(20 * sizeof(char));
+    uint_t pos = 0;
+    uint_t c = 0;
+    while (c < 5) {
+      if (fn[pos] == '_') {
+        c++;
+      }
+      if (i == 0) {
+        out_filename[pos] = fn[pos];
+        out_filename_len++;
+      }
+      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);
+
+    struct stat info;
+    stat(fn, &info);
+    uint_t num_bytes = info.st_size;
+    uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(double);
+
+    double *dist = malloc(dist_len * sizeof(double));
+    FILE *dist_file = fopen(fn, "rb");
+    fread(dist, sizeof(double), dist_len, dist_file);
+    fclose(dist_file);
+    {
+      double mom1[2];
+      get_mean(dist_len, dist, mom1);
+      vals_c1[i] = mom1[0];
+      errors_c1[i] = mom1[1];
+
+      double mom2[2];
+      get_mom(dist_len, 2, dist, mom1, mom2);
+      vals_c2[i] = mom2[0];
+      errors_c2[i] = mom2[1];
+
+      double mom3[2];
+      get_mom(dist_len, 3, dist, mom1, mom3);
+      vals_c3[i] = mom3[0];
+      errors_c3[i] = mom3[1];
+
+      double mom4[2];
+      get_mom(dist_len, 4, dist, mom1, mom4);
+      vals_c4[i] = mom4[0];
+      errors_c4[i] = mom4[1];
+    }
+    free(dist);
+  }
+
+  out_filename[out_filename_len - 1] = '.';
+  out_filename[out_filename_len] = 't';
+  out_filename[out_filename_len + 1] = 'x';
+  out_filename[out_filename_len + 2] = 't';
+  out_filename[out_filename_len + 3] = '\0';
+
+  FILE *out_file = fopen(out_filename, "w");
+
+  for (uint_t i = 0; i < nc; i++) {
+    fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+            vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+            errors_c3[i], vals_c4[i], errors_c4[i]);
+  }
+
+  fclose(out_file);
+
+  free(Ls_c);
+  free(betas_c);
+  free(vals_c1);
+  free(errors_c1);
+  free(vals_c2);
+  free(errors_c2);
+  free(vals_c3);
+  free(errors_c3);
+  free(vals_c4);
+  free(errors_c4);
+
+  return 0;
 }
-
diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c
index 3bf388c..ee2b029 100644
--- a/src/cen_anal_process.c
+++ b/src/cen_anal_process.c
@@ -1,154 +1,157 @@
 
 #include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.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>
 
 void get_mean(uint_t len, uint32_t *data, double result[2]) {
-	uint_t total = 0;
-	double mean = 0;
-	double mean_err = 0;
+  uint_t total = 0;
+  double mean = 0;
+  double mean_err = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		uint32_t datai = data[i];
+  for (uint_t i = 0; i < len; i++) {
+    uint32_t datai = data[i];
 
-		total += datai;
-		mean += i * datai;
-		mean_err += gsl_pow_2(i) * datai; 
-	}
+    total += datai;
+    mean += i * datai;
+    mean_err += gsl_pow_2(i) * datai;
+  }
 
-	double meanf = mean / total;
-	double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
+  double meanf = mean / total;
+  double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total);
 
-	result[0] = meanf;
-	result[1] = meanf_err;
+  result[0] = meanf;
+  result[1] = meanf_err;
 }
 
-void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], double result[2]) {
-	uint_t total = 0;
-	double mom = 0;
-	double mom_err = 0;
+void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2],
+             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(fabs(((double)i) - mean[0]), n);
+  for (uint_t i = 0; i < len; i++) {
+    uint32_t datai = data[i];
+    double in = pow(fabs(((double)i) - mean[0]), n);
 
-		total += datai;
-		mom += in * datai;
-		mom_err += gsl_pow_2(in) * datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0])));
-	}
+    total += datai;
+    mom += in * datai;
+    mom_err += gsl_pow_2(in) *
+               datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0])));
+  }
 
-	double momf = mom / total;
-	double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
+  double momf = mom / total;
+  double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
 
-	result[0] = momf;
-	result[1] = momf_err;
+  result[0] = momf;
+  result[1] = momf_err;
 }
 
 int main(int argc, char *argv[]) {
-	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_c1 = (double *)malloc(nc * sizeof(double));
-	double *errors_c1 = (double *)malloc(nc * sizeof(double));
-	double *vals_c2 = (double *)malloc(nc * sizeof(double));
-	double *errors_c2 = (double *)malloc(nc * sizeof(double));
-	double *vals_c3 = (double *)malloc(nc * sizeof(double));
-	double *errors_c3 = (double *)malloc(nc * sizeof(double));
-	double *vals_c4 = (double *)malloc(nc * sizeof(double));
-	double *errors_c4 = (double *)malloc(nc * sizeof(double));
-
-	char *out_filename = (char *)malloc(100 * sizeof(char));
-	uint_t out_filename_len = 0;
-
-	for (uint_t i = 0; i < nc; i++) {
-		char *fn = argv[1 + i];
-		char *buff = (char *)malloc(20 * sizeof(char));
-		uint_t pos = 0; uint_t c = 0;
-		while (c < 5) {
-			if (fn[pos] == '_') {
-				c++;
-			}
-			if (i == 0) {
-				out_filename[pos] = fn[pos];
-				out_filename_len++;
-			}
-			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 mom1[2];
-			get_mean(dist_len, dist, mom1);
-			vals_c1[i] = mom1[0];
-			errors_c1[i] = mom1[1];
-
-			double mom2[2];
-			get_mom(dist_len, 2, dist, mom1, mom2);
-			vals_c2[i] = mom2[0];
-			errors_c2[i] = mom2[1];
-
-			double mom3[2];
-			get_mom(dist_len, 3, dist, mom1, mom3);
-			vals_c3[i] = mom3[0];
-			errors_c3[i] = mom3[1];
-
-			double mom4[2];
-			get_mom(dist_len, 4, dist, mom1, mom4);
-			vals_c4[i] = mom4[0];
-			errors_c4[i] = mom4[1];
-		}
-		free(dist);
-	}
-
-	out_filename[out_filename_len-1] = '.';
-	out_filename[out_filename_len] = 't';
-	out_filename[out_filename_len+1] = 'x';
-	out_filename[out_filename_len+2] = 't';
-	out_filename[out_filename_len+3] = '\0';
-
-	FILE *out_file = fopen(out_filename, "w");
-
-	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
-	}
-
-	fclose(out_file);
-
-
-	free(Ls_c);
-	free(betas_c);
-	free(vals_c1);
-	free(errors_c1);
-	free(vals_c2);
-	free(errors_c2);
-	free(vals_c3);
-	free(errors_c3);
-
-	return 0;
+  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_c1 = (double *)malloc(nc * sizeof(double));
+  double *errors_c1 = (double *)malloc(nc * sizeof(double));
+  double *vals_c2 = (double *)malloc(nc * sizeof(double));
+  double *errors_c2 = (double *)malloc(nc * sizeof(double));
+  double *vals_c3 = (double *)malloc(nc * sizeof(double));
+  double *errors_c3 = (double *)malloc(nc * sizeof(double));
+  double *vals_c4 = (double *)malloc(nc * sizeof(double));
+  double *errors_c4 = (double *)malloc(nc * sizeof(double));
+
+  char *out_filename = (char *)malloc(100 * sizeof(char));
+  uint_t out_filename_len = 0;
+
+  for (uint_t i = 0; i < nc; i++) {
+    char *fn = argv[1 + i];
+    char *buff = (char *)malloc(20 * sizeof(char));
+    uint_t pos = 0;
+    uint_t c = 0;
+    while (c < 5) {
+      if (fn[pos] == '_') {
+        c++;
+      }
+      if (i == 0) {
+        out_filename[pos] = fn[pos];
+        out_filename_len++;
+      }
+      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 mom1[2];
+      get_mean(dist_len, dist, mom1);
+      vals_c1[i] = mom1[0];
+      errors_c1[i] = mom1[1];
+
+      double mom2[2];
+      get_mom(dist_len, 2, dist, mom1, mom2);
+      vals_c2[i] = mom2[0];
+      errors_c2[i] = mom2[1];
+
+      double mom3[2];
+      get_mom(dist_len, 3, dist, mom1, mom3);
+      vals_c3[i] = mom3[0];
+      errors_c3[i] = mom3[1];
+
+      double mom4[2];
+      get_mom(dist_len, 4, dist, mom1, mom4);
+      vals_c4[i] = mom4[0];
+      errors_c4[i] = mom4[1];
+    }
+    free(dist);
+  }
+
+  out_filename[out_filename_len - 1] = '.';
+  out_filename[out_filename_len] = 't';
+  out_filename[out_filename_len + 1] = 'x';
+  out_filename[out_filename_len + 2] = 't';
+  out_filename[out_filename_len + 3] = '\0';
+
+  FILE *out_file = fopen(out_filename, "w");
+
+  for (uint_t i = 0; i < nc; i++) {
+    fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i],
+            vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i],
+            errors_c3[i], vals_c4[i], errors_c4[i]);
+  }
+
+  fclose(out_file);
+
+  free(Ls_c);
+  free(betas_c);
+  free(vals_c1);
+  free(errors_c1);
+  free(vals_c2);
+  free(errors_c2);
+  free(vals_c3);
+  free(errors_c3);
+
+  return 0;
 }
-
diff --git a/src/corr_test.c b/src/corr_test.c
index 01b3e11..cb00212 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -2,26 +2,26 @@
 #include "fracture.h"
 
 int main() {
-	cholmod_common c;
-	CHOL_F(start)(&c);
+  cholmod_common c;
+  CHOL_F(start)(&c);
 
-	unsigned int width = 64;
+  unsigned int width = 64;
 
-	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);
+  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);
 
-	double *corr = get_corr(instance, NULL, &c);
+  double *corr = get_corr(instance, NULL, &c);
 
-	for (int i = 0; i < 2 * width; i++) {
-		printf("%g ", corr[i]);
-	}
-	printf("\n");
+  for (int i = 0; i < 2 * width; i++) {
+    printf("%g ", corr[i]);
+  }
+  printf("\n");
 
-	net_free(instance, &c);
-	graph_free(network);
+  net_free(instance, &c);
+  graph_free(network);
 
-	CHOL_F(finish)(&c);
+  CHOL_F(finish)(&c);
 
-	return 0;
+  return 0;
 }
diff --git a/src/fracture.c b/src/fracture.c
index 73059ff..ede7a24 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -3,510 +3,515 @@
 #include "fracture.h"
 
 int main(int argc, char *argv[]) {
-	int opt;
-
-	// defining variables to be (potentially) set by command line flags
-	uint8_t filename_len;
-	uint32_t N;
-	uint_t L;
-	double beta, inf, cutoff, crack_len;
-	bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
-			 save_crit_stress, save_energy, save_conductivity,
-			 save_damage, save_threshold, save_current_load;
-	bound_t boundary;
-	lattice_t lattice;
-
-
-	// assume filenames will be less than 100 characters
-
-	filename_len = 100;
-
-
-	// set default values
-
-	N = 100;
-	L = 16;
-	crack_len = 0.;
-	beta = .3;
-	inf = 1e10;
-	cutoff = 1e-9;
-	boundary = FREE_BOUND;
-	lattice = VORONOI_LATTICE;
-	save_data = false;
-	save_cluster_dist = false;
-	use_voltage_boundaries = false;
-	use_dual = false;
-	save_network = false;
-	save_crit_stress = false;
-	save_damage = false;
-	save_conductivity = false;
-	save_energy = false;
-	save_threshold = false;
-	save_current_load = false;
-
-
-	uint8_t bound_i;
-	char boundc2 = 'f';
-	uint8_t lattice_i;
-	char lattice_c = 'v';
-	char dual_c = 'o';
-
-
-	// get commandline options
-
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
-		switch (opt) {
-			case 'n':
-				N = atoi(optarg);
-				break;
-			case 'L':
-				L = atoi(optarg);
-				break;
-			case 'b':
-				beta = atof(optarg);
-				break;
-			case 'l':
-				crack_len = atof(optarg);
-				break;
-			case 'B':
-				bound_i = atoi(optarg);
-				switch (bound_i) {
-					case 0:
-						boundary = FREE_BOUND;
-						boundc2 = 'f';
-						break;
-					case 1:
-						boundary = CYLINDER_BOUND;
-						boundc2 = 'c';
-						break;
-					case 2:
-						boundary = TORUS_BOUND;
-						use_voltage_boundaries = true;
-						boundc2 = 't';
-						break;
-					case 3:
-						boundary = EMBEDDED_BOUND;
-						boundc2 = 'e';
-						use_dual = true;
-						use_voltage_boundaries = true;
-						break;
-					default:
-						printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), or 2 (TORUS_BOUND).\n");
-						exit(EXIT_FAILURE);
-				}
-				break;
-			case 'q':
-				lattice_i = atoi(optarg);
-				switch (lattice_i) {
-					case 0:
-						lattice = VORONOI_LATTICE;
-						lattice_c = 'v';
-						break;
-					case 1:
-						lattice = DIAGONAL_LATTICE;
-						lattice_c = 'd';
-						break;
-					case 2:
-						lattice = VORONOI_HYPERUNIFORM_LATTICE;
-						lattice_c = 'h';
-						break;
-					case 3:
-						lattice = TRIANGULAR_LATTICE;
-						lattice_c = 't';
-						break;
-          case 4:
-            lattice = SQUARE_LATTICE;
-            lattice_c = 's';
-            break;
-					default:
-						printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
-						exit(EXIT_FAILURE);
-				}
-				break;
-			case 'd':
-				save_damage = true;
-				break;
-			case 'V':
-				use_voltage_boundaries = true;
-				break;
-			case 'D':
-				use_dual = true;
-				dual_c = 'd';
-				break;
-			case 'c':
-				save_cluster_dist = true;
-				break;
-			case 'o':
-				save_data = true;
-				break;
-			case 'N':
-				save_network = true;
-				break;
-			case 's':
-				save_crit_stress = true;
-				break;
-			case 'r':
-				save_conductivity = true;
-				break;
-			case 'E':
-				save_energy = true;
-				break;
-			case 'T':
-				save_threshold = true;
-				break;
-			case 'C':
-				save_current_load = true;
-				break;
-			default: /* '?' */
-				exit(EXIT_FAILURE);
-		}
-	}
-
-
-	char boundc;
-	if (use_voltage_boundaries) boundc = 'v';
-	else boundc = 'c';
-
-	FILE *data_out;
-	if (save_data) {
-		char *data_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		data_out = fopen(data_filename, "a");
-		free(data_filename);
-	}
-
-	uint_t max_verts, max_edges;
-
-	// these are very liberal estimates
-	max_verts = 4 * pow(L, 2);
-	max_edges = 4 * pow(L, 2);
-
-	if (max_verts > CINT_MAX) {
-		exit(EXIT_FAILURE);
-	}
-
-	// define arrays for saving cluster and avalanche distributions
-	uint32_t *cluster_size_dist;
-	uint32_t *avalanche_size_dist;
-	char *c_filename;
-	char *a_filename;
-	if (save_cluster_dist) {
-		cluster_size_dist =
-				(uint32_t *)calloc(max_verts, sizeof(uint32_t));
-		avalanche_size_dist =
-				(uint32_t *)calloc(max_edges, sizeof(uint32_t));
-
-		c_filename = (char *)malloc(filename_len * sizeof(char));
-		a_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-
-		FILE *cluster_out = fopen(c_filename, "rb");
-		FILE *avalanche_out = fopen(a_filename, "rb");
-
-		if (cluster_out != NULL) {
-			fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
-			fclose(cluster_out);
-		}
-		if (avalanche_out != NULL) {
-			fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
-			fclose(avalanche_out);
-		}
-	}
-
-	long double *crit_stress;
-	if (save_crit_stress) {
-		crit_stress = (long double *)malloc(N * sizeof(long double));
-	}
-
-	double *conductivity;
-	if (save_conductivity) {
-		conductivity = (double *)malloc(N * sizeof(double));
-	}
-
-	// define arrays for saving damage distributions
-	uint32_t *damage;
-	char *d_filename;
-	if (save_damage) {
-		damage =
-				(uint32_t *)calloc(max_edges, sizeof(uint32_t));
-
-		d_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-
-		FILE *damage_out = fopen(d_filename, "rb");
-
-		if (damage_out != NULL) {
-			fread(damage, sizeof(uint32_t), max_edges, damage_out);
-			fclose(damage_out);
-		}
-	}
-
-	long double *energy;
-	if (save_energy) {
-		energy = (long double *)malloc(N * sizeof(long double));
-	}
-
-	long double *thresholds;
-	if (save_threshold) {
-		thresholds = (long double *)malloc(N * sizeof(long double));
-	}
-
-	long double *loads;
-	if (save_current_load) {
-		loads = (long double *)malloc(N * sizeof(long double));
-	}
-
-
-	// start cholmod
-	cholmod_common c;
-	CHOL_F(start)(&c);
-
-	/* if we use voltage boundary conditions, the laplacian matrix is positive
-	 * definite and we can use a supernodal LL decomposition.  otherwise we need
-	 * to use the simplicial LDL decomposition
-	 */
-	if (use_voltage_boundaries) {
-		//(&c)->supernodal = CHOLMOD_SUPERNODAL;
-		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
-	} else {
-		(&c)->supernodal = CHOLMOD_SIMPLICIAL;
-	}
-
-
-	printf("\n");
-	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);
-		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);
-		net_free(tmp_net, &c);
-
-		uint_t max_pos = 0;
-		long double max_val = 0;
-
-		double cond0;
-		{
-			double *tmp_voltages = net_voltages(net, &c);
-			cond0 = net_conductivity(net, tmp_voltages);
-			free(tmp_voltages);
-		}
-
-		for (uint_t j = 0; j < data->num_broken; j++) {
-			long double val = data->extern_field[j];
-
-			if (val > max_val) {
-				max_pos = j;
-				max_val = val;
-			}
-		}
-
-		uint_t av_size = 0;
-		long double cur_val = 0;
-
-		for (uint_t j = 0; j < max_pos; j++) {
-			uint_t next_broken = data->break_list[j];
-
-			break_edge(net, next_broken, &c);
-
-			long double val = data->extern_field[j];
-			if (save_cluster_dist) {
-				if (val < cur_val) {
-					av_size++;
-				}
-
-				if (val > cur_val) {
-					avalanche_size_dist[av_size]++;
-					av_size = 0;
-					cur_val = val;
-				}
-			}
-		}
-
-		if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
-
-		if (save_conductivity) {
-			if (max_pos > 0) {
-				conductivity[i] = data->conductivity[max_pos - 1];
-			} else {
-				conductivity[i] = cond0;
-			}
-		}
-
-		if (save_damage) {
-			uint_t would_break = 0;
-			double *tmp_voltage = net_voltages(net, &c);
-			double *tmp_current = net_currents(net, tmp_voltage, &c);
-			free(tmp_voltage);
-			for (uint_t j = 0; j < g->ne; j++) {
-				bool broken = net->fuses[j];
-				bool under_thres = net->thres[j] < net->thres[data->break_list[max_pos]];
-				bool zero_field = fabs(tmp_current[j]) < cutoff;
-				if (!broken && under_thres && zero_field) {
-					break_edge(net, j, &c);
-				}
-			}
-			damage[net->num_broken]++;
-			free(tmp_current);
-		}
-
-		if (save_energy) {
-			long double tmp_energy = 0;
-			if (max_pos > 0) {
-				long double sigma1 = data->extern_field[0];
-				double cond1 = cond0;
-				for (uint_t j = 0; j < max_pos - 1; j++) {
-					long double sigma2 = data->extern_field[j+1];
-					double cond2 = data->conductivity[j];
-					if (sigma2 > sigma1) {
-						tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
-						sigma1 = sigma2;
-						cond1 = cond2;
-					}
-				}
-			}
-			energy[i] = tmp_energy;
-		}
-
-		if (save_threshold) {
-			thresholds[i] = net->thres[data->break_list[max_pos]];
-		}
-
-		if (save_current_load) {
-			loads[i] = data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
-		}
-
-		if (save_data) {
-			for (uint_t j = 0; j < data->num_broken; j++) {
-				fprintf(data_out, "%u %Lg %g ", data->break_list[j],
-								data->extern_field[j], data->conductivity[j]);
-			}
-			fprintf(data_out, "\n");
-		}
-
-		data_free(data);
-		if (save_network) {
-			FILE *net_out = fopen("network.txt", "w");
-			for (uint_t j = 0; j < g->nv; j++) {
-				fprintf(net_out, "%f %f ", g->vx[2 * j],
-						g->vx[2 * j + 1]);
-			}
-			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < g->ne; j++) {
-				fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
-			}
-			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < g->dnv; j++) {
-				fprintf(net_out, "%f %f ", g->dvx[2 * j],
-						g->dvx[2 * j + 1]);
-			}
-			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < g->ne; j++) {
-				fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
-			}
-			fprintf(net_out, "\n");
-			for (uint_t j = 0; j < g->ne; j++) {
-				fprintf(net_out, "%d ", net->fuses[j]);
-			}
-			fclose(net_out);
-		}
-
-		if (save_cluster_dist) {
-			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];
-			}
-			free(tmp_cluster_dist);
-		}
-
-
-		net_free(net, &c);
-		graph_free(g);
-	}
-
-	printf("\033[F\033[JFRACTURE: COMPLETE\n");
-
-	if (save_cluster_dist) {
-		FILE *cluster_out = fopen(c_filename, "wb");
-		FILE *avalanche_out = fopen(a_filename, "wb");
-
-		fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
-		fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
-
-		fclose(cluster_out);
-		fclose(avalanche_out);
-
-		free(c_filename);
-		free(a_filename);
-		free(cluster_size_dist);
-		free(avalanche_size_dist);
-	}
-
-	if (save_conductivity) {
-		char *cond_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *cond_file = fopen(cond_filename, "ab");
-		fwrite(conductivity, sizeof(double), N, cond_file);
-		fclose(cond_file);
-		free(cond_filename);
-		free(conductivity);
-	}
-
-	if (save_energy) {
-		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *tough_file = fopen(tough_filename, "ab");
-		fwrite(energy, sizeof(long double), N, tough_file);
-		fclose(tough_file);
-		free(tough_filename);
-		free(energy);
-	}
-
-	if (save_threshold) {
-		char *thres_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *thres_file = fopen(thres_filename, "ab");
-		fwrite(thresholds, sizeof(long double), N, thres_file);
-		fclose(thres_file);
-		free(thres_filename);
-		free(thresholds);
-	}
-
-	if (save_current_load) {
-		char *load_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *load_file = fopen(load_filename, "ab");
-		fwrite(loads, sizeof(long double), N, load_file);
-		fclose(load_file);
-		free(load_filename);
-		free(loads);
-	}
-
-	if (save_damage) {
-		FILE *hdam_file = fopen(d_filename, "wb");
-		fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
-		fclose(hdam_file);
-		free(d_filename);
-		free(damage);
-	}
-
-	if (save_data) {
-		fclose(data_out);
-	}
-
-	if (save_crit_stress) {
-		char *str_filename = (char *)malloc(filename_len * sizeof(char));
-		snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
-		FILE *str_file = fopen(str_filename, "ab");
-		fwrite(crit_stress, sizeof(long double), N, str_file);
-		fclose(str_file);
-		free(str_filename);
-		free(crit_stress);
-	}
-
-	CHOL_F(finish)(&c);
-
-	return 0;
+  int opt;
+
+  // defining variables to be (potentially) set by command line flags
+  uint8_t filename_len;
+  uint32_t N;
+  uint_t L;
+  double beta, inf, cutoff, crack_len;
+  bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual,
+      save_network, save_crit_stress, save_energy, save_conductivity,
+      save_damage, save_threshold, save_current_load;
+  bound_t boundary;
+  lattice_t lattice;
+
+  // assume filenames will be less than 100 characters
+
+  filename_len = 100;
+
+  // set default values
+
+  N = 100;
+  L = 16;
+  crack_len = 0.;
+  beta = .3;
+  inf = 1e10;
+  cutoff = 1e-9;
+  boundary = FREE_BOUND;
+  lattice = VORONOI_LATTICE;
+  save_data = false;
+  save_cluster_dist = false;
+  use_voltage_boundaries = false;
+  use_dual = false;
+  save_network = false;
+  save_crit_stress = false;
+  save_damage = false;
+  save_conductivity = false;
+  save_energy = false;
+  save_threshold = false;
+  save_current_load = false;
+
+  uint8_t bound_i;
+  char boundc2 = 'f';
+  uint8_t lattice_i;
+  char lattice_c = 'v';
+  char dual_c = 'o';
+
+  // get commandline options
+
+  while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TE")) != -1) {
+    switch (opt) {
+    case 'n':
+      N = atoi(optarg);
+      break;
+    case 'L':
+      L = atoi(optarg);
+      break;
+    case 'b':
+      beta = atof(optarg);
+      break;
+    case 'l':
+      crack_len = atof(optarg);
+      break;
+    case 'B':
+      bound_i = atoi(optarg);
+      switch (bound_i) {
+      case 0:
+        boundary = FREE_BOUND;
+        boundc2 = 'f';
+        break;
+      case 1:
+        boundary = CYLINDER_BOUND;
+        boundc2 = 'c';
+        break;
+      case 2:
+        boundary = TORUS_BOUND;
+        use_voltage_boundaries = true;
+        boundc2 = 't';
+        break;
+      case 3:
+        boundary = EMBEDDED_BOUND;
+        boundc2 = 'e';
+        use_dual = true;
+        use_voltage_boundaries = true;
+        break;
+      default:
+        printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), "
+               "or 2 (TORUS_BOUND).\n");
+        exit(EXIT_FAILURE);
+      }
+      break;
+    case 'q':
+      lattice_i = atoi(optarg);
+      switch (lattice_i) {
+      case 0:
+        lattice = VORONOI_LATTICE;
+        lattice_c = 'v';
+        break;
+      case 1:
+        lattice = DIAGONAL_LATTICE;
+        lattice_c = 'd';
+        break;
+      case 2:
+        lattice = VORONOI_HYPERUNIFORM_LATTICE;
+        lattice_c = 'h';
+        break;
+      case 3:
+        lattice = TRIANGULAR_LATTICE;
+        lattice_c = 't';
+        break;
+      case 4:
+        lattice = SQUARE_LATTICE;
+        lattice_c = 's';
+        break;
+      default:
+        printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 "
+               "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n");
+        exit(EXIT_FAILURE);
+      }
+      break;
+    case 'd':
+      save_damage = true;
+      break;
+    case 'V':
+      use_voltage_boundaries = true;
+      break;
+    case 'D':
+      use_dual = true;
+      dual_c = 'd';
+      break;
+    case 'c':
+      save_cluster_dist = true;
+      break;
+    case 'o':
+      save_data = true;
+      break;
+    case 'N':
+      save_network = true;
+      break;
+    case 's':
+      save_crit_stress = true;
+      break;
+    case 'r':
+      save_conductivity = true;
+      break;
+    case 'E':
+      save_energy = true;
+      break;
+    case 'T':
+      save_threshold = true;
+      break;
+    case 'C':
+      save_current_load = true;
+      break;
+    default: /* '?' */
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  char boundc;
+  if (use_voltage_boundaries)
+    boundc = 'v';
+  else
+    boundc = 'c';
+
+  FILE *data_out;
+  if (save_data) {
+    char *data_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    data_out = fopen(data_filename, "a");
+    free(data_filename);
+  }
+
+  uint_t max_verts, max_edges;
+
+  // these are very liberal estimates
+  max_verts = 4 * pow(L, 2);
+  max_edges = 4 * pow(L, 2);
+
+  if (max_verts > CINT_MAX) {
+    exit(EXIT_FAILURE);
+  }
+
+  // define arrays for saving cluster and avalanche distributions
+  uint32_t *cluster_size_dist;
+  uint32_t *avalanche_size_dist;
+  char *c_filename;
+  char *a_filename;
+  if (save_cluster_dist) {
+    cluster_size_dist = (uint32_t *)calloc(max_verts, sizeof(uint32_t));
+    avalanche_size_dist = (uint32_t *)calloc(max_edges, sizeof(uint32_t));
+
+    c_filename = (char *)malloc(filename_len * sizeof(char));
+    a_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+
+    FILE *cluster_out = fopen(c_filename, "rb");
+    FILE *avalanche_out = fopen(a_filename, "rb");
+
+    if (cluster_out != NULL) {
+      fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
+      fclose(cluster_out);
+    }
+    if (avalanche_out != NULL) {
+      fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
+      fclose(avalanche_out);
+    }
+  }
+
+  long double *crit_stress;
+  if (save_crit_stress) {
+    crit_stress = (long double *)malloc(N * sizeof(long double));
+  }
+
+  double *conductivity;
+  if (save_conductivity) {
+    conductivity = (double *)malloc(N * sizeof(double));
+  }
+
+  // define arrays for saving damage distributions
+  uint32_t *damage;
+  char *d_filename;
+  if (save_damage) {
+    damage = (uint32_t *)calloc(max_edges, sizeof(uint32_t));
+
+    d_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+
+    FILE *damage_out = fopen(d_filename, "rb");
+
+    if (damage_out != NULL) {
+      fread(damage, sizeof(uint32_t), max_edges, damage_out);
+      fclose(damage_out);
+    }
+  }
+
+  long double *energy;
+  if (save_energy) {
+    energy = (long double *)malloc(N * sizeof(long double));
+  }
+
+  long double *thresholds;
+  if (save_threshold) {
+    thresholds = (long double *)malloc(N * sizeof(long double));
+  }
+
+  long double *loads;
+  if (save_current_load) {
+    loads = (long double *)malloc(N * sizeof(long double));
+  }
+
+  // start cholmod
+  cholmod_common c;
+  CHOL_F(start)(&c);
+
+  /* if we use voltage boundary conditions, the laplacian matrix is positive
+   * definite and we can use a supernodal LL decomposition.  otherwise we need
+   * to use the simplicial LDL decomposition
+   */
+  if (use_voltage_boundaries) {
+    //(&c)->supernodal = CHOLMOD_SUPERNODAL;
+    (&c)->supernodal = CHOLMOD_SIMPLICIAL;
+  } else {
+    (&c)->supernodal = CHOLMOD_SIMPLICIAL;
+  }
+
+  printf("\n");
+  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);
+    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);
+    net_free(tmp_net, &c);
+
+    uint_t max_pos = 0;
+    long double max_val = 0;
+
+    double cond0;
+    {
+      double *tmp_voltages = net_voltages(net, &c);
+      cond0 = net_conductivity(net, tmp_voltages);
+      free(tmp_voltages);
+    }
+
+    for (uint_t j = 0; j < data->num_broken; j++) {
+      long double val = data->extern_field[j];
+
+      if (val > max_val) {
+        max_pos = j;
+        max_val = val;
+      }
+    }
+
+    uint_t av_size = 0;
+    long double cur_val = 0;
+
+    for (uint_t j = 0; j < max_pos; j++) {
+      uint_t next_broken = data->break_list[j];
+
+      break_edge(net, next_broken, &c);
+
+      long double val = data->extern_field[j];
+      if (save_cluster_dist) {
+        if (val < cur_val) {
+          av_size++;
+        }
+
+        if (val > cur_val) {
+          avalanche_size_dist[av_size]++;
+          av_size = 0;
+          cur_val = val;
+        }
+      }
+    }
+
+    if (save_crit_stress)
+      crit_stress[i] = data->extern_field[max_pos];
+
+    if (save_conductivity) {
+      if (max_pos > 0) {
+        conductivity[i] = data->conductivity[max_pos - 1];
+      } else {
+        conductivity[i] = cond0;
+      }
+    }
+
+    if (save_damage) {
+      uint_t would_break = 0;
+      double *tmp_voltage = net_voltages(net, &c);
+      double *tmp_current = net_currents(net, tmp_voltage, &c);
+      free(tmp_voltage);
+      for (uint_t j = 0; j < g->ne; j++) {
+        bool broken = net->fuses[j];
+        bool under_thres =
+            net->thres[j] < net->thres[data->break_list[max_pos]];
+        bool zero_field = fabs(tmp_current[j]) < cutoff;
+        if (!broken && under_thres && zero_field) {
+          break_edge(net, j, &c);
+        }
+      }
+      damage[net->num_broken]++;
+      free(tmp_current);
+    }
+
+    if (save_energy) {
+      long double tmp_energy = 0;
+      if (max_pos > 0) {
+        long double sigma1 = data->extern_field[0];
+        double cond1 = cond0;
+        for (uint_t j = 0; j < max_pos - 1; j++) {
+          long double sigma2 = data->extern_field[j + 1];
+          double cond2 = data->conductivity[j];
+          if (sigma2 > sigma1) {
+            tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
+            sigma1 = sigma2;
+            cond1 = cond2;
+          }
+        }
+      }
+      energy[i] = tmp_energy;
+    }
+
+    if (save_threshold) {
+      thresholds[i] = net->thres[data->break_list[max_pos]];
+    }
+
+    if (save_current_load) {
+      loads[i] =
+          data->extern_field[max_pos] / net->thres[data->break_list[max_pos]];
+    }
+
+    if (save_data) {
+      for (uint_t j = 0; j < data->num_broken; j++) {
+        fprintf(data_out, "%u %Lg %g ", data->break_list[j],
+                data->extern_field[j], data->conductivity[j]);
+      }
+      fprintf(data_out, "\n");
+    }
+
+    data_free(data);
+    if (save_network) {
+      FILE *net_out = fopen("network.txt", "w");
+      for (uint_t j = 0; j < g->nv; j++) {
+        fprintf(net_out, "%f %f ", g->vx[2 * j], g->vx[2 * j + 1]);
+      }
+      fprintf(net_out, "\n");
+      for (uint_t j = 0; j < g->ne; j++) {
+        fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]);
+      }
+      fprintf(net_out, "\n");
+      for (uint_t j = 0; j < g->dnv; j++) {
+        fprintf(net_out, "%f %f ", g->dvx[2 * j], g->dvx[2 * j + 1]);
+      }
+      fprintf(net_out, "\n");
+      for (uint_t j = 0; j < g->ne; j++) {
+        fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]);
+      }
+      fprintf(net_out, "\n");
+      for (uint_t j = 0; j < g->ne; j++) {
+        fprintf(net_out, "%d ", net->fuses[j]);
+      }
+      fclose(net_out);
+    }
+
+    if (save_cluster_dist) {
+      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];
+      }
+      free(tmp_cluster_dist);
+    }
+
+    net_free(net, &c);
+    graph_free(g);
+  }
+
+  printf("\033[F\033[JFRACTURE: COMPLETE\n");
+
+  if (save_cluster_dist) {
+    FILE *cluster_out = fopen(c_filename, "wb");
+    FILE *avalanche_out = fopen(a_filename, "wb");
+
+    fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out);
+    fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out);
+
+    fclose(cluster_out);
+    fclose(avalanche_out);
+
+    free(c_filename);
+    free(a_filename);
+    free(cluster_size_dist);
+    free(avalanche_size_dist);
+  }
+
+  if (save_conductivity) {
+    char *cond_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    FILE *cond_file = fopen(cond_filename, "ab");
+    fwrite(conductivity, sizeof(double), N, cond_file);
+    fclose(cond_file);
+    free(cond_filename);
+    free(conductivity);
+  }
+
+  if (save_energy) {
+    char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    FILE *tough_file = fopen(tough_filename, "ab");
+    fwrite(energy, sizeof(long double), N, tough_file);
+    fclose(tough_file);
+    free(tough_filename);
+    free(energy);
+  }
+
+  if (save_threshold) {
+    char *thres_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    FILE *thres_file = fopen(thres_filename, "ab");
+    fwrite(thresholds, sizeof(long double), N, thres_file);
+    fclose(thres_file);
+    free(thres_filename);
+    free(thresholds);
+  }
+
+  if (save_current_load) {
+    char *load_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    FILE *load_file = fopen(load_filename, "ab");
+    fwrite(loads, sizeof(long double), N, load_file);
+    fclose(load_file);
+    free(load_filename);
+    free(loads);
+  }
+
+  if (save_damage) {
+    FILE *hdam_file = fopen(d_filename, "wb");
+    fwrite(damage, sizeof(uint32_t), max_edges, hdam_file);
+    fclose(hdam_file);
+    free(d_filename);
+    free(damage);
+  }
+
+  if (save_data) {
+    fclose(data_out);
+  }
+
+  if (save_crit_stress) {
+    char *str_filename = (char *)malloc(filename_len * sizeof(char));
+    snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat",
+             lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+    FILE *str_file = fopen(str_filename, "ab");
+    fwrite(crit_stress, sizeof(long double), N, str_file);
+    fclose(str_file);
+    free(str_filename);
+    free(crit_stress);
+  }
+
+  CHOL_F(finish)(&c);
+
+  return 0;
 }
diff --git a/src/fracture.h b/src/fracture.h
index b1114fb..7b4f146 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -17,7 +17,6 @@
 #include <stdbool.h>
 #include <stdint.h>
 #include <stdio.h>
-#include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <sys/types.h>
@@ -26,7 +25,6 @@
 #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
@@ -36,68 +34,69 @@
 #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;
+  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;
+  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);
+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_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 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);
+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);
+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_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_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);
+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);
+                   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);
+                               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);
+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);
 
@@ -112,11 +111,13 @@ 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);
+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);
+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);
 
diff --git a/src/long_anal_process.c b/src/long_anal_process.c
index ba29152..d4ec4e0 100644
--- a/src/long_anal_process.c
+++ b/src/long_anal_process.c
@@ -1,156 +1,158 @@
 
 #include "fracture.h"
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_matrix.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 <sys/stat.h>
 
 void get_mean(uint_t len, long double *data, long double result[2]) {
-	long double total = 0;
+  long double total = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total += data[i];
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total += data[i];
+  }
 
-	long double mean = total / len;
-	long double total_sq_diff = 0;
+  long double mean = total / len;
+  long double total_sq_diff = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total_sq_diff += pow(mean - data[i], 2);
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total_sq_diff += pow(mean - data[i], 2);
+  }
 
-	long double mean_err = sqrt(total_sq_diff) / len;
+  long double mean_err = sqrt(total_sq_diff) / len;
 
-	result[0] = mean;
-	result[1] = mean_err;
+  result[0] = mean;
+  result[1] = mean_err;
 }
 
-void get_mom(uint_t len, uint_t n, long double *data, long double mean[2], long double result[2]) {
-	long double total_n_diff = 0;
-	long double total_n2_diff = 0;
+void get_mom(uint_t len, uint_t n, long double *data, long double mean[2],
+             long double result[2]) {
+  long double total_n_diff = 0;
+  long double total_n2_diff = 0;
 
-	for (uint_t i = 0; i < len; i++) {
-		total_n_diff += pow(fabsl(mean[0] - data[i]), n);
-		total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n);
-	}
+  for (uint_t i = 0; i < len; i++) {
+    total_n_diff += pow(fabsl(mean[0] - data[i]), n);
+    total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n);
+  }
 
-	long double mom = total_n_diff / len;
-	long double mom_err = sqrt(total_n2_diff) / len;
+  long double mom = total_n_diff / len;
+  long double mom_err = sqrt(total_n2_diff) / len;
 
-	result[0] = mom;
-	result[1] = mom_err;
+  result[0] = mom;
+  result[1] = mom_err;
 }
 
 int main(int argc, char *argv[]) {
-	uint_t nc = argc - 1;
-	uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
-	double *betas_c = (double *)malloc(nc * sizeof(double));
-	long double *vals_c1 = (long double *)malloc(nc * sizeof(long double));
-	long double *errors_c1 = (long double *)malloc(nc * sizeof(long double));
-	long double *vals_c2 = (long double *)malloc(nc * sizeof(long double));
-	long double *errors_c2 = (long double *)malloc(nc * sizeof(long double));
-	long double *vals_c3 = (long double *)malloc(nc * sizeof(long double));
-	long double *errors_c3 = (long double *)malloc(nc * sizeof(long double));
-	long double *vals_c4 = (long double *)malloc(nc * sizeof(long double));
-	long double *errors_c4 = (long double *)malloc(nc * sizeof(long double));
-
-	char *out_filename = (char *)malloc(100 * sizeof(char));
-	uint_t out_filename_len = 0;
-
-	for (uint_t i = 0; i < nc; i++) {
-		char *fn = argv[1 + i];
-		char *buff = (char *)malloc(20 * sizeof(char));
-		uint_t pos = 0; uint_t c = 0;
-		while (c < 5) {
-			if (fn[pos] == '_') {
-				c++;
-			}
-			if (i == 0) {
-				out_filename[pos] = fn[pos];
-				out_filename_len++;
-			}
-			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);
-
-		struct stat info;
-		stat(fn, &info);
-		uint_t num_bytes = info.st_size;
-		uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double);
-
-		long double *dist = malloc(dist_len * sizeof(long double));
-		FILE *dist_file = fopen(fn, "rb");
-		fread(dist, sizeof(long double), dist_len, dist_file);
-		fclose(dist_file);
-		{
-			long double mom1[2];
-			get_mean(dist_len, dist, mom1);
-			vals_c1[i] = mom1[0];
-			errors_c1[i] = mom1[1];
-
-			long double mom2[2];
-			get_mom(dist_len, 2, dist, mom1, mom2);
-			vals_c2[i] = mom2[0];
-			errors_c2[i] = mom2[1];
-
-			long double mom3[2];
-			get_mom(dist_len, 3, dist, mom1, mom3);
-			vals_c3[i] = mom3[0];
-			errors_c3[i] = mom3[1];
-
-			long double mom4[2];
-			get_mom(dist_len, 4, dist, mom1, mom4);
-			vals_c4[i] = mom4[0];
-			errors_c4[i] = mom4[1];
-		}
-		free(dist);
-	}
-
-	out_filename[out_filename_len-1] = '.';
-	out_filename[out_filename_len] = 't';
-	out_filename[out_filename_len+1] = 'x';
-	out_filename[out_filename_len+2] = 't';
-	out_filename[out_filename_len+3] = '\0';
-
-	FILE *out_file = fopen(out_filename, "w");
-
-	for (uint_t i = 0; i < nc; i++) {
-		fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i], betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
-	}
-
-	fclose(out_file);
-
-
-	free(Ls_c);
-	free(betas_c);
-	free(vals_c1);
-	free(errors_c1);
-	free(vals_c2);
-	free(errors_c2);
-	free(vals_c3);
-	free(errors_c3);
-	free(vals_c4);
-	free(errors_c4);
-
-	return 0;
+  uint_t nc = argc - 1;
+  uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
+  double *betas_c = (double *)malloc(nc * sizeof(double));
+  long double *vals_c1 = (long double *)malloc(nc * sizeof(long double));
+  long double *errors_c1 = (long double *)malloc(nc * sizeof(long double));
+  long double *vals_c2 = (long double *)malloc(nc * sizeof(long double));
+  long double *errors_c2 = (long double *)malloc(nc * sizeof(long double));
+  long double *vals_c3 = (long double *)malloc(nc * sizeof(long double));
+  long double *errors_c3 = (long double *)malloc(nc * sizeof(long double));
+  long double *vals_c4 = (long double *)malloc(nc * sizeof(long double));
+  long double *errors_c4 = (long double *)malloc(nc * sizeof(long double));
+
+  char *out_filename = (char *)malloc(100 * sizeof(char));
+  uint_t out_filename_len = 0;
+
+  for (uint_t i = 0; i < nc; i++) {
+    char *fn = argv[1 + i];
+    char *buff = (char *)malloc(20 * sizeof(char));
+    uint_t pos = 0;
+    uint_t c = 0;
+    while (c < 5) {
+      if (fn[pos] == '_') {
+        c++;
+      }
+      if (i == 0) {
+        out_filename[pos] = fn[pos];
+        out_filename_len++;
+      }
+      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);
+
+    struct stat info;
+    stat(fn, &info);
+    uint_t num_bytes = info.st_size;
+    uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double);
+
+    long double *dist = malloc(dist_len * sizeof(long double));
+    FILE *dist_file = fopen(fn, "rb");
+    fread(dist, sizeof(long double), dist_len, dist_file);
+    fclose(dist_file);
+    {
+      long double mom1[2];
+      get_mean(dist_len, dist, mom1);
+      vals_c1[i] = mom1[0];
+      errors_c1[i] = mom1[1];
+
+      long double mom2[2];
+      get_mom(dist_len, 2, dist, mom1, mom2);
+      vals_c2[i] = mom2[0];
+      errors_c2[i] = mom2[1];
+
+      long double mom3[2];
+      get_mom(dist_len, 3, dist, mom1, mom3);
+      vals_c3[i] = mom3[0];
+      errors_c3[i] = mom3[1];
+
+      long double mom4[2];
+      get_mom(dist_len, 4, dist, mom1, mom4);
+      vals_c4[i] = mom4[0];
+      errors_c4[i] = mom4[1];
+    }
+    free(dist);
+  }
+
+  out_filename[out_filename_len - 1] = '.';
+  out_filename[out_filename_len] = 't';
+  out_filename[out_filename_len + 1] = 'x';
+  out_filename[out_filename_len + 2] = 't';
+  out_filename[out_filename_len + 3] = '\0';
+
+  FILE *out_file = fopen(out_filename, "w");
+
+  for (uint_t i = 0; i < nc; i++) {
+    fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i],
+            betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i],
+            vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]);
+  }
+
+  fclose(out_file);
+
+  free(Ls_c);
+  free(betas_c);
+  free(vals_c1);
+  free(errors_c1);
+  free(vals_c2);
+  free(errors_c2);
+  free(vals_c3);
+  free(errors_c3);
+  free(vals_c4);
+  free(errors_c4);
+
+  return 0;
 }
-