system size now can take on non-integer values

8 files changed, 166 insertions, 103 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index d746ddd..c1ad4f0 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -36,6 +36,6 @@ class graph {
     std::vector<edge> dual_edges;
 
     graph(unsigned int Nx, unsigned int Ny);
-    graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double = 0.25);
+    graph(double Lx, double Ly, std::mt19937& rng, double relax = 0.01, double step = 1.9);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 10f62b7..3d5647f 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -8,6 +8,14 @@
 #define JCV_FLT_MAX 1.7976931348623157E+308
 #include <jc_voronoi.h>
 
+double mod(double a, double b) {
+  if (a >= 0) {
+    return fmod(a, b);
+  } else {
+    return fmod(a + b * ceil(-a / b), b);
+  }
+}
+
 graph::graph(unsigned int Nx, unsigned int Ny) {
   L = {(double)Nx, (double)Ny};
 
@@ -84,45 +92,97 @@ class eulerException: public std::exception
   }
 } eulerex;
 
+graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
+  L = {Lx, Ly};
 
-graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
-  L = {(double)Nx, (double)Ny};
+  std::uniform_real_distribution<double> d(0.0, 1.0);
+  unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5);
 
-  unsigned int nv = Nx * Ny / 2;
+  unsigned int nv = N;
 
   vertices.resize(nv);
-  
-  std::normal_distribution<double> d(0.0, spread);
 
   // the coordinates of the lattice, from which a delaunay triangulation
   // and corresponding voronoi tessellation will be built
   for (unsigned int i = 0; i < nv; i++) {
-    double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
-    double ry = (double)(i / (Nx / 2)) + d(rng);
-    vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
-  }
-
-  // to make the resulting tessellation periodic, we tile eight (!) copies of
-  // the original points for a total of nine. note that the index of each
-  // point quotient 9 is equal to the original index (we'll use this to
-  // translate back)
-  std::vector<jcv_point> points;
-  points.reserve(9 * nv);
-  for (const vertex &v : vertices) {
-    points.push_back({v.r.x, v.r.y});
-    points.push_back({v.r.x + L.x, v.r.y + 0.0});
-    points.push_back({v.r.x - L.x, v.r.y + 0.0});
-    points.push_back({v.r.x + 0.0, v.r.y + L.y});
-    points.push_back({v.r.x + 0.0, v.r.y - L.y});
-    points.push_back({v.r.x + L.x, v.r.y + L.y});
-    points.push_back({v.r.x - L.x, v.r.y + L.y});
-    points.push_back({v.r.x + L.x, v.r.y - L.y});
-    points.push_back({v.r.x - L.x, v.r.y - L.y});
+    vertices[i] = {{L.x * d(rng), L.y * d(rng)}};
   }
 
+  double max_difference = std::numeric_limits<double>::max();
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
+  std::vector<jcv_point> points(9 * nv);
+
+  double rstep = sqrt(step);
+
+  while (max_difference > pow(relax, 2) * 2 * N) {
+    double cur_diff = 0;
+    // to make the resulting tessellation periodic, we tile eight (!) copies of
+    // the original points for a total of nine. note that the index of each
+    // point quotient 9 is equal to the original index (we'll use this to
+    // translate back)
+    for (unsigned int i = 0; i < nv; i++) {
+      const vertex& v = vertices[i];
+      points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
+      points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
+      points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
+      points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
+      points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
+      points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+      points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
+      points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
+      points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+    }
+
+    jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
+    const jcv_site* sites = jcv_diagram_get_sites(&diagram);
+    
+    for (int i = 0; i < diagram.numsites; i++) {
+      const jcv_site* site = &sites[i];
+      unsigned int ind = site->index;
+      if (ind % 9 == 0) {
+        double Cx = 0.0;
+        double Cy = 0.0;
+        unsigned int ne = 0;
+
+        const jcv_graphedge* e = site->edges;
+        while (e) {
+          ne++;
+          Cx += e->pos[0].x;
+          Cy += e->pos[0].y;
+          e = e->next;
+        }
 
+        Cx /= ne;
+        Cy /= ne;
+
+        double dx = Cx - vertices[ind / 9].r.x;
+        double dy = Cy - vertices[ind / 9].r.y;
+
+        double dist = pow(dx, 2) + pow(dy, 2);
+        if (dist > cur_diff) {
+          cur_diff = dist;
+        }
+        vertices[ind / 9] = {{mod(vertices[ind / 9].r.x + rstep * dx, L.x), mod(vertices[ind / 9].r.y + rstep * dy, L.y)}};
+      }
+    }
+    max_difference = cur_diff;
+    jcv_diagram_free(&diagram);
+    memset(&diagram, 0, sizeof(jcv_diagram));
+  }
+
+  for (unsigned int i = 0; i < nv; i++) {
+    const vertex& v = vertices[i];
+    points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
+    points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
+    points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
+    points[9 * i + 3] = {v.r.x + 0.0, v.r.y + L.y};
+    points[9 * i + 4] = {v.r.x + 0.0, v.r.y - L.y};
+    points[9 * i + 5] = {v.r.x + L.x, v.r.y + L.y};
+    points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
+    points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
+    points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
+  }
   jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
@@ -167,7 +227,7 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
 
         if (it1 == known_vertices.end()) {
           vi1 = dual_vertices.size();
-          dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
+          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{{site->p.x, site->p.y}}});
           known_vertices[t1] = vi1;
         } else {
           vi1 = it1->second;
@@ -176,8 +236,8 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
         
         if (i1 < i2) {
           edges.push_back({{i1, i2},
-              {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x),
-               fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)}
+              {mod((site->p.x + neighbor->p.x) / 2, L.x),
+               mod((site->p.y + neighbor->p.y) / 2, L.y)}
                });
 
           jcv_graphedge *en;
@@ -197,15 +257,15 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread)
 
           if (it2 == known_vertices.end()) {
             vi2 = dual_vertices.size();
-            dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}});
+            dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}});
             known_vertices[t2] = vi2;
           } else {
             vi2 = it2->second;
           }
 
           dual_edges.push_back({{vi1, vi2},
-              {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x),
-               fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)}
+              {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
+               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)}
                });
         }
 
diff --git a/src/animate.cpp b/src/animate.cpp
index 8c2d12f..26f4996 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,15 +1,15 @@
 
 #include "animate.hpp"
 
-animate::animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]) : Lx(Lx), Ly(Ly), G(Lx * Ly) {
+animate::animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]) : G(2 * (unsigned int)ceil(Lx * Ly / 2)) {
   glutInit(&argc, argv);
   glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
-  glutInitWindowSize(window_size, (unsigned int )(window_size * ((double)Ly / (double)Lx)));
+  glutInitWindowSize((unsigned int)(Lx / Ly * window_size), window_size);
   glutCreateWindow("wolff");
   glClearColor(0.0,0.0,0.0,0.0);
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
-  gluOrtho2D(-1.0, Lx + 3.0, -1.0 , Ly + 3.0);
+  gluOrtho2D(-1, Lx + 3, -1 , Ly + 3);
 }
 
 void animate::pre_fracture(const network &n) {
@@ -25,18 +25,18 @@ void animate::pre_fracture(const network &n) {
     graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
     graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-    if (fabs(r1.x - r2.x) > Lx / 2) {
-      if (r1.x < Lx / 2) {
-        r1.x += Lx;
+    if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+      if (r1.x < n.G.L.x / 2) {
+        r1.x += n.G.L.x;
       } else {
-        r2.x += Lx;
+        r2.x += n.G.L.x;
       }
     }
-    if (fabs(r1.y - r2.y) > Ly / 2) {
-      if (r1.y < Ly / 2) {
-        r1.y += Ly;
+    if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+      if (r1.y < n.G.L.y / 2) {
+        r1.y += n.G.L.y;
       } else {
-        r2.y += Ly;
+        r2.y += n.G.L.y;
       }
     }
 
@@ -47,8 +47,8 @@ void animate::pre_fracture(const network &n) {
   glFlush();
 }
 
-void animate::bond_broken(const network& net, const current_info& cur, unsigned int i) {
-  long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i];
+void animate::bond_broken(const network& n, const current_info& cur, unsigned int i) {
+  long double c = logl(cur.conductivity / fabs(cur.currents[i])) + n.thresholds[i];
   if (c > lv && avalanches.back().size() > 0) {
     lv = c;
     avalanches.push_back({i});
@@ -56,25 +56,25 @@ void animate::bond_broken(const network& net, const current_info& cur, unsigned
     avalanches.back().push_back(i);
   }
 
-  boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+  boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G);
 
    glBegin(GL_LINES);          // Each set of 3 vertices form a triangle
     glColor3f(1.0f, 1.0f, 1.0f); // Blue
-    graph::coordinate r1 = net.G.vertices[net.G.edges[i].v[0]].r;
-    graph::coordinate r2 = net.G.vertices[net.G.edges[i].v[1]].r;
+    graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
+    graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-    if (fabs(r1.x - r2.x) > Lx / 2) {
-      if (r1.x < Lx / 2) {
-        r2.x -= Lx;
+    if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+      if (r1.x < n.G.L.x / 2) {
+        r2.x -= n.G.L.x;
       } else {
-        r2.x += Lx;
+        r2.x += n.G.L.x;
       }
     }
-    if (fabs(r1.y - r2.y) > Ly / 2) {
-      if (r1.y < Ly / 2) {
-        r2.y -= Ly;
+    if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+      if (r1.y < n.G.L.y / 2) {
+        r2.y -= n.G.L.y;
       } else {
-        r2.y += Ly;
+        r2.y += n.G.L.y;
       }
     }
 
@@ -111,18 +111,18 @@ void animate::post_fracture(network &n) {
         graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
         graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
 
-        if (fabs(r1.x - r2.x) > Lx / 2) {
-          if (r1.x < Lx / 2) {
-            r1.x += Lx;
+        if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+          if (r1.x < n.G.L.x / 2) {
+            r1.x += n.G.L.x;
           } else {
-            r2.x += Lx;
+            r2.x += n.G.L.x;
           }
         }
-        if (fabs(r1.y - r2.y) > Ly / 2) {
-          if (r1.y < Ly / 2) {
-            r1.y += Ly;
+        if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+          if (r1.y < n.G.L.y / 2) {
+            r1.y += n.G.L.y;
           } else {
-            r2.y += Ly;
+            r2.y += n.G.L.y;
           }
         }
 
diff --git a/src/animate.hpp b/src/animate.hpp
index 793771d..423eba7 100644
--- a/src/animate.hpp
+++ b/src/animate.hpp
@@ -9,14 +9,12 @@
 
 class animate : public hooks {
   private:
-    unsigned int Lx;
-    unsigned int Ly;
     Graph G;
   public:
     long double lv;
     std::list<std::list<unsigned int>> avalanches;
 
-    animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]);
+    animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]);
 
     void pre_fracture(const network &) override;
     void bond_broken(const network& net, const current_info& cur, unsigned int i) override;
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 10b63d0..62c5b1e 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -31,24 +31,28 @@ int main(int argc, char* argv[]) {
   int opt;
 
   unsigned int N = 1;
-  unsigned int Lx = 16;
-  unsigned int Ly = 16;
+  double Lx = 16.0;
+  double Ly = 16.0;
   double  beta = 0.5;
+  double w = 0.01;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+  while ((opt = getopt(argc, argv, "X:Y:N:b:w:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned int)atof(optarg);
         break;
       case 'X':
-        Lx = atoi(optarg);
+        Lx = atof(optarg);
         break;
       case 'Y':
-        Ly = atoi(optarg);
+        Ly = atof(optarg);
         break;
       case 'b':
         beta = atof(optarg);
         break;
+      case 'w':
+        w = atof(optarg);
+        break;
       default:
         exit(1);
     }
@@ -63,7 +67,7 @@ int main(int argc, char* argv[]) {
   std::mt19937 rng{seeds};
 
   for (unsigned int trial = 0; trial < N; trial++) {
-    graph G(Lx, Ly, rng, 0.4);
+    graph G(Lx, Ly, rng, w);
     network network(G, &c);
     network.set_thresholds(beta, rng);
     network.fracture(meas);
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 4b6fb5b..742a3f9 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -31,8 +31,8 @@ int main(int argc, char* argv[]) {
   int opt;
 
   unsigned int N = 1;
-  unsigned int Lx = 16;
-  unsigned int Ly = 16;
+  double Lx = 16;
+  double Ly = 16;
   double  beta = 0.5;
 
   while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
@@ -41,10 +41,10 @@ int main(int argc, char* argv[]) {
         N = (unsigned int)atof(optarg);
         break;
       case 'X':
-        Lx = atoi(optarg);
+        Lx = atof(optarg);
         break;
       case 'Y':
-        Ly = atoi(optarg);
+        Ly = atof(optarg);
         break;
       case 'b':
         beta = atof(optarg);
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  ma meas(Lx, Ly, Lx, Ly, beta);
+  ma meas(Lx, Ly, ceil(Lx), ceil(Ly), beta);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
@@ -65,7 +65,7 @@ int main(int argc, char* argv[]) {
   for (unsigned int trial = 0; trial < N; trial++) {
     while (true) {
       try {
-        graph G(Lx, Ly, rng, 0.4);
+        graph G(Lx, Ly, rng);
         network network(G, &c);
         network.set_thresholds(beta, rng);
         network.fracture(meas);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 43483be..27f239d 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -2,7 +2,7 @@
 #include "measurements.hpp"
 #include <iostream>
 
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, double Lx, double Ly, double beta) {
   std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
   std::ifstream file(filename);
 
@@ -31,7 +31,7 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
 }
 
 template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta, unsigned int Mx, unsigned int My) {
+void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, double Lx, double Ly, double beta, unsigned int Mx, unsigned int My) {
   std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
   std::ifstream file(filename);
 
@@ -57,16 +57,16 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
 }
 
 template <class T>
-std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
-  for (unsigned int i = 0; i < Lx * Ly; i++) {
+std::vector<fftw_complex> data_transform(unsigned int Mx, unsigned int My, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
+  for (unsigned int i = 0; i < Mx * My; i++) {
     fftw_forward_in[i] = (double)data[i];
   }
 
   fftw_execute(forward_plan);
 
-  std::vector<fftw_complex> output(Lx * (Ly / 2 + 1));
+  std::vector<fftw_complex> output(Mx * (My / 2 + 1));
 
-  for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
+  for (unsigned int i = 0; i < Mx * (My / 2 + 1); i++) {
     output[i][0] = fftw_forward_out[i][0];
     output[i][1] = fftw_forward_out[i][1];
   }
@@ -75,43 +75,44 @@ std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const
 }
 
 template <class T>
-void correlation(unsigned int Lx, unsigned int Ly, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
-  for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
+void correlation(unsigned int Mx, unsigned int My, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+  for (unsigned int i = 0; i < Mx * (My / 2 + 1); i++) {
     fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
     fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
   }
 
   fftw_execute(reverse_plan);
 
-  for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
-    data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
+  for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
+    data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My));
   }
 }
 
 template <class T>
-void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+void autocorrelation(unsigned int Mx, unsigned int My, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
   fftw_execute(forward_plan);
 
-  for (unsigned int i = 0; i < Ly * (Lx / 2 + 1); i++) {
+  for (unsigned int i = 0; i < My * (Mx / 2 + 1); i++) {
     fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
     fftw_reverse_in[i][1] = 0.0;
   }
 
   fftw_execute(reverse_plan);
 
-  for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
-    out_data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
+  for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
+    out_data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My));
   }
 }
 
-unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) {
+unsigned int edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned int Mx, unsigned int My) {
   return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
 }
 
-ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly),
-  sc(Lx * Ly, 0),
-  sa(Lx * Ly, 0),
-  sd(Lx * Ly, 0),
+ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
+  Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)),
+  sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sd(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sb(log2(Mx < My ? Mx : My) + 1, 0),
   Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
   Css((Mx / 2 + 1) * (My / 2 + 1), 0),
@@ -217,7 +218,7 @@ void ma::post_fracture(network &n) {
   // non-spanning clusters
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
   for (unsigned int i = 0; i < num; i++) {
-    if (i != crack_component) {
+    if (i != crack_component && components[i].size() > 0) {
       for (auto it = components[i].begin(); it != components[i].end(); it++) {
         fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
       }
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 685806b..bd93cfc 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -26,8 +26,8 @@ class ma : public hooks {
     fftw_plan forward_plan;
     fftw_plan reverse_plan;
     unsigned int N;
-    unsigned int Lx;
-    unsigned int Ly;
+    double Lx;
+    double Ly;
     unsigned int Mx;
     unsigned int My;
     double beta;
@@ -54,7 +54,7 @@ class ma : public hooks {
 
     std::list<std::list<unsigned int>> avalanches;
 
-    ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta);
+    ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta);
     ~ma();
 
     void pre_fracture(const network &) override;