added many correlation types

1 files changed, 291 insertions, 53 deletions

diff --git a/src/fracture.c b/src/fracture.c
index 49fb039..5b44238 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -175,14 +175,25 @@ int main(int argc, char *argv[]) {
   }
 
   double *dd_correlations; // damage-damage
-  double *dp_correlations; // damage-precursor
+  double *dc_correlations; // damage-crack
+  double *db_correlations; // damage-backbone
   double *ds_correlations; // damage-stress
-  double *pp_correlations; // precursor-precursor
-  double *ps_correlations; // precursor-stress
+  double *dA_correlations; // damage-avalanche
+  double *cc_correlations; // crack-crack
+  double *cb_correlations; // crack-backbone
+  double *cs_correlations; // crack-stress
+  double *cA_correlations; // crack-avalanche
+  double *bb_correlations; // backbone-backbone
+  double *bs_correlations; // backbone-stress
+  double *bA_correlations; // backbone-avalanche
   double *ss_correlations; // stress-stress
+  double *AA_correlations; // avalanche-avalanche
   double *DD_correlations; // after-crack damage-damage
-  double *DP_correlations; // after-crack damage-precursor
-  double *PP_correlations; // after-crack precursor-precursor
+  double *DC_correlations; // after-crack damage-crack
+  double *DB_correlations; // after-crack damage-backbone
+  double *CC_correlations; // after-crack crack-crack
+  double *CB_correlations; // after-crack crack-backbone
+  double *BB_correlations; // after-crack backbone-backbone
   double *fftw_forward_in;
   fftw_complex *fftw_forward_out;
   fftw_complex *fftw_reverse_in;
@@ -191,22 +202,36 @@ int main(int argc, char *argv[]) {
   fftw_plan reverse_plan;
   uint64_t N_averaged = 0;
   double mean_D = 0;
-  double mean_P = 0;
+  double mean_A = 0;
+  double mean_B = 0;
+  double mean_C = 0;
   double mean_d = 0;
-  double mean_p = 0;
+  double mean_c = 0;
+  double mean_b = 0;
   double mean_s = 0;
   char *correlations_filename;
   if (save_correlations) {
     assert(lattice == DIAGONAL_LATTICE);
     dd_correlations = (double *)calloc(pow(L, 2), sizeof(double));
-    dp_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    dc_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    db_correlations = (double *)calloc(pow(L, 2), sizeof(double));
     ds_correlations = (double *)calloc(pow(L, 2), sizeof(double));
-    pp_correlations = (double *)calloc(pow(L, 2), sizeof(double));
-    ps_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    dA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    cc_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    cb_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    cs_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    cA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    bb_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    bs_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    bA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
     ss_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    AA_correlations = (double *)calloc(pow(L, 2), sizeof(double));
     DD_correlations = (double *)calloc(pow(L, 2), sizeof(double));
-    DP_correlations = (double *)calloc(pow(L, 2), sizeof(double));
-    PP_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    DC_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    DB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    CC_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    CB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
+    BB_correlations = (double *)calloc(pow(L, 2), sizeof(double));
     fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double));
     fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
     fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex));
@@ -223,19 +248,33 @@ int main(int argc, char *argv[]) {
     if (correlations_out != NULL) {
       fread(&N_averaged, sizeof(uint64_t), 1, correlations_out);
       fread(&mean_d, sizeof(double), 1, correlations_out);
-      fread(&mean_p, sizeof(double), 1, correlations_out);
+      fread(&mean_c, sizeof(double), 1, correlations_out);
+      fread(&mean_b, sizeof(double), 1, correlations_out);
       fread(&mean_s, sizeof(double), 1, correlations_out);
+      fread(&mean_A, sizeof(double), 1, correlations_out);
       fread(&mean_D, sizeof(double), 1, correlations_out);
-      fread(&mean_P, sizeof(double), 1, correlations_out);
+      fread(&mean_C, sizeof(double), 1, correlations_out);
+      fread(&mean_B, sizeof(double), 1, correlations_out);
       fread(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
-      fread(dp_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(dc_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(db_correlations, sizeof(double), pow(L, 2), correlations_out);
       fread(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
-      fread(pp_correlations, sizeof(double), pow(L, 2), correlations_out);
-      fread(ps_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(dA_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(cc_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(cb_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(cs_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(cA_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(bb_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(bs_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(bA_correlations, sizeof(double), pow(L, 2), correlations_out);
       fread(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(AA_correlations, sizeof(double), pow(L, 2), correlations_out);
       fread(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
-      fread(DP_correlations, sizeof(double), pow(L, 2), correlations_out);
-      fread(PP_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(DC_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(DB_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(CC_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(CB_correlations, sizeof(double), pow(L, 2), correlations_out);
+      fread(BB_correlations, sizeof(double), pow(L, 2), correlations_out);
       fclose(correlations_out);
     }
   }
@@ -428,8 +467,37 @@ int main(int argc, char *argv[]) {
       }
 
       fftw_execute(forward_plan);
-      fftw_complex *P_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
-      memcpy(P_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+      fftw_complex *C_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
+      memcpy(C_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        fftw_forward_in[j] = 0.0;
+      }
+
+      uint32_t in_backbone1 = 0;
+      dll_t *tmp_cycle = cycle;
+      while (tmp_cycle != NULL) {
+        fftw_forward_in[tmp_cycle->e] = 1.0;
+        in_backbone1++;
+        tmp_cycle = tmp_cycle->right;
+      }
+
+      fftw_execute(forward_plan);
+      fftw_complex *B_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
+      memcpy(B_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+      uint32_t in_avalanche = 0;
+      for (uint32_t j = 0; j < g->ne; j++) {
+        if (tmp_net->fuses[j] && !(net->fuses[j])) {
+          fftw_forward_in[j] = 1.0;
+        } else {
+          fftw_forward_in[j] = 0.0;
+        }
+      }
+
+      fftw_execute(forward_plan);
+      fftw_complex *A_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
+      memcpy(A_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
 
       uint32_t damage2 = 0;
       for (uint32_t j = 0; j < g->ne; j++) {
@@ -456,11 +524,28 @@ int main(int argc, char *argv[]) {
       }
 
       fftw_execute(forward_plan);
-      fftw_complex *p_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex));
-      memcpy(p_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+      fftw_complex *c_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
+      memcpy(c_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        fftw_forward_in[j] = 0.0;
+      }
+
+      uint32_t in_backbone2 = 0;
+      tmp_cycle = cycle;
+      while (tmp_cycle != NULL) {
+        if (net->fuses[tmp_cycle->e]) {
+          fftw_forward_in[tmp_cycle->e] = 1.0;
+          in_backbone2++;
+        }
+        tmp_cycle = tmp_cycle->right;
+      }
+
+      fftw_execute(forward_plan);
+      fftw_complex *b_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex));
+      memcpy(b_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
 
       graph_components_free(comp);
-      int test = 0;
       while (cycle != NULL) {
         dll_t *old = cycle;
         cycle = cycle->right;
@@ -484,9 +569,12 @@ int main(int argc, char *argv[]) {
       memcpy(s_transform, fftw_forward_out, g->ne * sizeof(fftw_complex));
 
       mean_D = (double)damage1 / (1.0 + N_averaged) + (double)N_averaged * mean_D / (N_averaged + 1.0);
-      mean_P = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_P / (N_averaged + 1.0);
+      mean_C = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_C / (N_averaged + 1.0);
+      mean_B = (double)in_backbone1 / (1.0 + N_averaged) + (double)N_averaged * mean_B / (N_averaged + 1.0);
+      mean_A = (double)in_avalanche / (1.0 + N_averaged) + (double)N_averaged * mean_A / (N_averaged + 1.0);
       mean_d = (double)damage2 / (1.0 + N_averaged) + (double)N_averaged * mean_d / (N_averaged + 1.0);
-      mean_p = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_p / (N_averaged + 1.0);
+      mean_c = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_c / (N_averaged + 1.0);
+      mean_b = (double)in_backbone2 / (1.0 + N_averaged) + (double)N_averaged * mean_b / (N_averaged + 1.0);
       mean_s = (double)t_stress / (1.0 + N_averaged) + (double)N_averaged * mean_s / (N_averaged + 1.0);
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
@@ -501,25 +589,69 @@ int main(int argc, char *argv[]) {
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
-        fftw_reverse_in[j][0] = D_transform[j][0] * P_transform[j][0] + D_transform[j][1] * P_transform[j][1];
-        fftw_reverse_in[j][1] = D_transform[j][0] * P_transform[j][1] - D_transform[j][1] * P_transform[j][0];
+        fftw_reverse_in[j][0] = D_transform[j][0] * C_transform[j][0] + D_transform[j][1] * C_transform[j][1];
+        fftw_reverse_in[j][1] = D_transform[j][0] * C_transform[j][1] - D_transform[j][1] * C_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        DC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DC_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = D_transform[j][0] * B_transform[j][0] + D_transform[j][1] * B_transform[j][1];
+        fftw_reverse_in[j][1] = D_transform[j][0] * B_transform[j][1] - D_transform[j][1] * B_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        DB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DB_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = C_transform[j][0] * C_transform[j][0] + C_transform[j][1] * C_transform[j][1];
+        fftw_reverse_in[j][1] = 0.0;
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        CC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CC_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = C_transform[j][0] * B_transform[j][0] + C_transform[j][1] * B_transform[j][1];
+        fftw_reverse_in[j][1] = C_transform[j][0] * B_transform[j][1] - C_transform[j][1] * B_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        CB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CB_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = B_transform[j][0] * B_transform[j][0] + B_transform[j][1] * B_transform[j][1];
+        fftw_reverse_in[j][1] = 0.0;
       }
 
       fftw_execute(reverse_plan);
 
       for (uint32_t j = 0; j < g->ne; j++) {
-        DP_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DP_correlations[j] / (N_averaged + 1.0);
+        BB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * BB_correlations[j] / (N_averaged + 1.0);
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
-        fftw_reverse_in[j][0] = P_transform[j][0] * P_transform[j][0] + P_transform[j][1] * P_transform[j][1];
+        fftw_reverse_in[j][0] = A_transform[j][0] * A_transform[j][0] + A_transform[j][1] * A_transform[j][1];
         fftw_reverse_in[j][1] = 0.0;
       }
 
       fftw_execute(reverse_plan);
 
       for (uint32_t j = 0; j < g->ne; j++) {
-        PP_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * PP_correlations[j] / (N_averaged + 1.0);
+        AA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * AA_correlations[j] / (N_averaged + 1.0);
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
@@ -534,14 +666,14 @@ int main(int argc, char *argv[]) {
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
-        fftw_reverse_in[j][0] = d_transform[j][0] * p_transform[j][0] + d_transform[j][1] * p_transform[j][1];
-        fftw_reverse_in[j][1] = d_transform[j][0] * p_transform[j][1] - d_transform[j][1] * p_transform[j][0];
+        fftw_reverse_in[j][0] = d_transform[j][0] * c_transform[j][0] + d_transform[j][1] * c_transform[j][1];
+        fftw_reverse_in[j][1] = d_transform[j][0] * c_transform[j][1] - d_transform[j][1] * c_transform[j][0];
       }
 
       fftw_execute(reverse_plan);
 
       for (uint32_t j = 0; j < g->ne; j++) {
-        dp_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dp_correlations[j] / (N_averaged + 1.0);
+        dc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dc_correlations[j] / (N_averaged + 1.0);
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
@@ -556,27 +688,105 @@ int main(int argc, char *argv[]) {
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
-        fftw_reverse_in[j][0] = p_transform[j][0] * p_transform[j][0] + p_transform[j][1] * p_transform[j][1];
+        fftw_reverse_in[j][0] = d_transform[j][0] * b_transform[j][0] + d_transform[j][1] * b_transform[j][1];
+        fftw_reverse_in[j][1] = d_transform[j][0] * b_transform[j][1] - d_transform[j][1] * b_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        db_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * db_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = d_transform[j][0] * A_transform[j][0] + d_transform[j][1] * A_transform[j][1];
+        fftw_reverse_in[j][1] = d_transform[j][0] * A_transform[j][1] - d_transform[j][1] * A_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        dA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dA_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = c_transform[j][0] * c_transform[j][0] + c_transform[j][1] * c_transform[j][1];
+        fftw_reverse_in[j][1] = 0.0;
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        cc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cc_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = c_transform[j][0] * s_transform[j][0] + c_transform[j][1] * s_transform[j][1];
+        fftw_reverse_in[j][1] = c_transform[j][0] * s_transform[j][1] - c_transform[j][1] * s_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        cs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cs_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = c_transform[j][0] * b_transform[j][0] + c_transform[j][1] * b_transform[j][1];
+        fftw_reverse_in[j][1] = c_transform[j][0] * b_transform[j][1] - c_transform[j][1] * b_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        cb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cb_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = c_transform[j][0] * A_transform[j][0] + c_transform[j][1] * A_transform[j][1];
+        fftw_reverse_in[j][1] = c_transform[j][0] * A_transform[j][1] - c_transform[j][1] * A_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        cA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cA_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = b_transform[j][0] * b_transform[j][0] + b_transform[j][1] * b_transform[j][1];
         fftw_reverse_in[j][1] = 0.0;
       }
 
       fftw_execute(reverse_plan);
 
       for (uint32_t j = 0; j < g->ne; j++) {
-        pp_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * pp_correlations[j] / (N_averaged + 1.0);
+        bb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bb_correlations[j] / (N_averaged + 1.0);
+      }
+
+      for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
+        fftw_reverse_in[j][0] = b_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1];
+        fftw_reverse_in[j][1] = b_transform[j][0] * s_transform[j][1] - s_transform[j][1] * s_transform[j][0];
+      }
+
+      fftw_execute(reverse_plan);
+
+      for (uint32_t j = 0; j < g->ne; j++) {
+        bs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bs_correlations[j] / (N_averaged + 1.0);
       }
 
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
-        fftw_reverse_in[j][0] = p_transform[j][0] * s_transform[j][0] + p_transform[j][1] * s_transform[j][1];
-        fftw_reverse_in[j][1] = p_transform[j][0] * s_transform[j][1] - p_transform[j][1] * s_transform[j][0];
+        fftw_reverse_in[j][0] = A_transform[j][0] * b_transform[j][0] + A_transform[j][1] * b_transform[j][1];
+        fftw_reverse_in[j][1] = A_transform[j][0] * b_transform[j][1] - A_transform[j][1] * b_transform[j][0];
       }
 
       fftw_execute(reverse_plan);
 
       for (uint32_t j = 0; j < g->ne; j++) {
-        ps_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ps_correlations[j] / (N_averaged + 1.0);
+        bA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bA_correlations[j] / (N_averaged + 1.0);
       }
 
+
       for (uint32_t j = 0; j < L * (L / 2 + 1); j++) {
         fftw_reverse_in[j][0] = s_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1];
         fftw_reverse_in[j][1] = 0.0;
@@ -589,9 +799,12 @@ int main(int argc, char *argv[]) {
       }
 
       free(D_transform);
-      free(P_transform);
+      free(C_transform);
+      free(A_transform);
+      free(B_transform);
+      free(b_transform);
       free(d_transform);
-      free(p_transform);
+      free(c_transform);
       free(s_transform);
       N_averaged++;
     }
@@ -789,30 +1002,55 @@ int main(int argc, char *argv[]) {
     FILE *correlations_out = fopen(correlations_filename, "wb");
     fwrite(&N_averaged, sizeof(uint64_t), 1, correlations_out);
     fwrite(&mean_d, sizeof(double), 1, correlations_out);
-    fwrite(&mean_p, sizeof(double), 1, correlations_out);
+    fwrite(&mean_c, sizeof(double), 1, correlations_out);
+    fwrite(&mean_b, sizeof(double), 1, correlations_out);
     fwrite(&mean_s, sizeof(double), 1, correlations_out);
+    fwrite(&mean_A, sizeof(double), 1, correlations_out);
     fwrite(&mean_D, sizeof(double), 1, correlations_out);
-    fwrite(&mean_P, sizeof(double), 1, correlations_out);
+    fwrite(&mean_C, sizeof(double), 1, correlations_out);
+    fwrite(&mean_B, sizeof(double), 1, correlations_out);
     fwrite(dd_correlations, sizeof(double), pow(L, 2), correlations_out);
-    fwrite(dp_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(dc_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(db_correlations, sizeof(double), pow(L, 2), correlations_out);
     fwrite(ds_correlations, sizeof(double), pow(L, 2), correlations_out);
-    fwrite(pp_correlations, sizeof(double), pow(L, 2), correlations_out);
-    fwrite(ps_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(dA_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(cc_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(cb_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(cs_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(cA_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(bb_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(bs_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(bA_correlations, sizeof(double), pow(L, 2), correlations_out);
     fwrite(ss_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(AA_correlations, sizeof(double), pow(L, 2), correlations_out);
     fwrite(DD_correlations, sizeof(double), pow(L, 2), correlations_out);
-    fwrite(DP_correlations, sizeof(double), pow(L, 2), correlations_out);
-    fwrite(PP_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(DC_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(DB_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(CC_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(CB_correlations, sizeof(double), pow(L, 2), correlations_out);
+    fwrite(BB_correlations, sizeof(double), pow(L, 2), correlations_out);
     fclose(correlations_out);
 
     free(dd_correlations);
-    free(dp_correlations);
+    free(dc_correlations);
+    free(db_correlations);
     free(ds_correlations);
-    free(pp_correlations);
-    free(ps_correlations);
+    free(dA_correlations);
+    free(cc_correlations);
+    free(cb_correlations);
+    free(cs_correlations);
+    free(cA_correlations);
+    free(bb_correlations);
+    free(bs_correlations);
+    free(bA_correlations);
     free(ss_correlations);
+    free(AA_correlations);
     free(DD_correlations);
-    free(DP_correlations);
-    free(PP_correlations);
+    free(DC_correlations);
+    free(DB_correlations);
+    free(CC_correlations);
+    free(CB_correlations);
+    free(BB_correlations);
     fftw_free(fftw_forward_in);
     fftw_free(fftw_forward_out);
     fftw_free(fftw_reverse_in);