New measurements for moduli jumps.

2 files changed, 55 insertions, 2 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 5132358..0a5471e 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -162,7 +162,7 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u
 
 ma::ma(unsigned n, double a, double beta, double weight)
     : sc(2 * n), sC(2 * n), sn(2 * n), ss(2 * n), sm(2 * n), sl(2 * n), sb(n + 1), sd(3 * n), sD(3 * n), sf(3 * n), sa(3 * n),
-      sA(3 * n), se(3 * n), sE(3 * n), si(10000), sI(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cC(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
+      sA(3 * n), se(3 * n), sE(3 * n), si(10000), sI(10000), sy(10000), sY(10000), sz(10000), sZ(10000), cc(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cC(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cn(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cs(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cm(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cl(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cb(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ca(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cA(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), ce(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), cE(pow(1+2*(unsigned)ceil(sqrt(n)), 2)),
       cq(pow(1+2*(unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n), last_clusters_v(2 * n) {
@@ -191,7 +191,7 @@ ma::ma(unsigned n, double a, double beta, double weight)
 
 ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
     : sc(Lx * Ly / 2), sC(Lx * Ly / 2), sn(Lx * Ly / 2), ss(Lx * Ly / 2), sm(Lx * Ly / 2), sl(Lx * Ly / 2),
-      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sD(Lx * Ly), sf(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), se(Lx * Ly), sE(Lx * Ly), si(10000), sI(10000),
+      sb(Lx * Ly / 2 + 1), sd(Lx * Ly), sD(Lx * Ly), sf(Lx * Ly), sa(Lx * Ly), sA(Lx * Ly), se(Lx * Ly), sE(Lx * Ly), si(10000), sI(10000), sy(10000), sY(10000), sz(10000), sZ(10000),
       cc((Lx / 2 + 1) * (Ly / 2 + 1)), cC((Lx / 2 + 1) * (Ly / 2 + 1)), cn((Lx / 2 + 1) * (Ly / 2 + 1)),
       cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)),
       cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)),
@@ -238,6 +238,10 @@ ma::~ma() {
   update_distribution_file("sE", sE, model_string);
   update_distribution_file("si", si, model_string);
   update_distribution_file("sI", sI, model_string);
+  update_distribution_file("sy", sy, model_string);
+  update_distribution_file("sY", sY, model_string);
+  update_distribution_file("sz", sz, model_string);
+  update_distribution_file("sZ", sZ, model_string);
   update_field_file("cc", Nc, cc, model_string);
   update_field_file("cC", NC, cC, model_string);
   update_field_file("cl", Nl, cl, model_string);
@@ -255,6 +259,8 @@ ma::~ma() {
 void ma::pre_fracture(const network&) {
   lc = std::numeric_limits<long double>::lowest();
   lv = std::numeric_limits<long double>::lowest();
+  moduliA = {};
+  moduliE = {};
   avalanches = {};
   evalanches = {};
   num = 0;
@@ -266,6 +272,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   if (c > lc) {
     last_clusters = net.C;
     lc = c;
+    moduliA.push_back(cur.conductivity[0]);
     avalanches.push_back({i});
   } else {
     avalanches.back().push_back(i);
@@ -273,6 +280,7 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   if (v > lv) {
     last_clusters_v = net.C;
     lv = v;
+    moduliE.push_back(cur.conductivity[0]);
     evalanches.push_back({i});
   } else {
     evalanches.back().push_back(i);
@@ -394,6 +402,44 @@ void ma::post_fracture(network& n) {
       av_it++;
     }
 
+    double mAFinal = log(moduliA.back());
+    if (mAFinal < 2 && mAFinal > -18) {
+      sZ[(unsigned)(10000 * (mAFinal + 18) / 20)]++;
+    }
+
+    auto am_it_curr = moduliA.rbegin();
+    auto am_it_next = std::next(am_it_curr);
+
+    while (am_it_next != moduliA.rend()) {
+      double mm = log(*am_it_next - *am_it_curr);
+
+      if (mm < 2 && mm > -18) {
+        sz[(unsigned)(10000 * (mm + 18) / 20)]++;
+      }
+
+      am_it_curr++;
+      am_it_next++;
+    }
+
+    double mEFinal = log(moduliE.back());
+    if (mEFinal < 2 && mEFinal > -18) {
+      sY[(unsigned)(10000 * (mEFinal + 18) / 20)]++;
+    }
+
+    auto em_it_curr = moduliE.rbegin();
+    auto em_it_next = std::next(em_it_curr);
+
+    while (em_it_next != moduliE.rend()) {
+      double mm = log(*em_it_next - *em_it_curr);
+
+      if (mm < 2 && mm > -18) {
+        sy[(unsigned)(10000 * (mm + 18) / 20)]++;
+      }
+
+      em_it_curr++;
+      em_it_next++;
+    }
+
     auto ev_it = evalanches.rbegin();
     ev_it++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 14e6580..a4ed5d7 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -51,6 +51,10 @@ class ma : public hooks {
 
     std::vector<uint64_t> si;
     std::vector<uint64_t> sI;
+    std::vector<uint64_t> sy;
+    std::vector<uint64_t> sY;
+    std::vector<uint64_t> sz;
+    std::vector<uint64_t> sZ;
 
     std::vector<uint64_t> cc;
     std::vector<uint64_t> cC;
@@ -71,7 +75,10 @@ class ma : public hooks {
   public:
     long double lc;
     long double lv;
+    long double lm;
 
+    std::list<double> moduliA;
+    std::list<double> moduliE;
     std::list<std::list<unsigned>> avalanches;
     std::list<std::list<unsigned>> evalanches;
     std::list<unsigned> last_avalanche;