1 files changed, 36 insertions, 3 deletions

diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index d718440..76831ba 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -40,6 +40,7 @@ int main(int argc, char *argv[]) {
 
   bool silent = false;
   bool use_pert = false;
+  bool N_is_sweeps = false;
 
   bool modulated_field = false;
   int order = 2;
@@ -51,7 +52,7 @@ int main(int argc, char *argv[]) {
 
   unsigned char measurement_flags = measurement_energy | measurement_clusterSize;
 
-  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:M:S")) != -1) {
     switch (opt) {
     case 'N': // number of steps
       N = (count_t)atof(optarg);
@@ -82,11 +83,14 @@ int main(int argc, char *argv[]) {
       modulated_field = true;
       break;
     case 'M':
-      measurement_flags |= 1 << atoi(optarg);
+      measurement_flags ^= 1 << atoi(optarg);
       break;
     case 'o':
       order = atoi(optarg);
       break;
+    case 'S':
+      N_is_sweeps = true;
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -179,6 +183,13 @@ int main(int argc, char *argv[]) {
     n_measurements++;
   }
 
+  meas_t *meas_sweeps;
+  if (N_is_sweeps) {
+    meas_sweeps = (meas_t *)calloc(1, sizeof(meas_t));
+    measurements[n_measurements] = measurement_average_cluster<orthogonal_R_t, vector_R_t> (meas_sweeps);
+    n_measurements++;
+  }
+
   std::function <double(vector_R_t)> H;
 
   if (modulated_field) {
@@ -187,7 +198,24 @@ int main(int argc, char *argv[]) {
     H = std::bind(H_vector <N_COMP, double>, std::placeholders::_1, H_vec);
   }
 
-  wolff <orthogonal_R_t, vector_R_t> (N, D, L, T, dot <N_COMP, double>, H, gen_R, n_measurements, measurements, silent);
+  // initialize random number generator
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, rand_seed());
+
+  state_t <orthogonal_R_t, vector_R_t> s(D, L, T, dot <N_COMP, double>, H);
+
+  if (N_is_sweeps) {
+    count_t N_rounds = 0;
+    printf("\n");
+    while (N_rounds * N * meas_sweeps->x < N * s.nv) {
+      printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+      wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, measurements, r, silent);
+      N_rounds++;
+    }
+    printf("\033[F\033[J\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n\n", (count_t)(N_rounds * N * meas_sweeps->x / s.nv), N, s.E, s.last_cluster_size);
+  } else {
+    wolff <orthogonal_R_t, vector_R_t> (N, &s, gen_R, n_measurements, measurements, r, silent);
+  }
 
   free(measurements);
 
@@ -204,7 +232,12 @@ int main(int argc, char *argv[]) {
     fclose(outfile_F);
   }
 
+  if (N_is_sweeps) {
+    free(meas_sweeps);
+  }
+
   free(H_vec);
+  gsl_rng_free(r);
 
   return 0;
 }