From 10510b7e2a09c102c307165e0cd716e17b17ef84 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 5 Sep 2018 22:45:39 -0400
Subject: fixed a bug in perturbative generation and set up epsilon to be
 automatically set by T and H

---
 lib/orthogonal.h |  2 +-
 src/wolff_On.cpp | 16 ++++++++++++++++
 2 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 829df10..4ae3ca6 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -182,7 +182,7 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const vecto
 
   v2 = 0;
 
-  double factor = epsilon * gsl_ran_ugaussian(r);
+  double factor = gsl_ran_gaussian(r, epsilon);
 
   for (q_t i = 0; i < q; i++) {
     m[0][i] = m[0][i] + factor * v[i]; // perturb orthogonal vector in original direction
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 6cd8777..736eb1b 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -142,6 +142,22 @@ int main(int argc, char *argv[]) {
   std::function <orthogonal_R_t(gsl_rng *, vector_R_t)> gen_R;
 
   if (use_pert) {
+    double Hish;
+    if (modulated_field) {
+      Hish = fabs(H_vec[0]);
+    } else {
+      double H2 = 0;
+      for (q_t i = 0; i < N_COMP; i++) {
+        H2 += pow(H_vec[i], 2);
+      }
+      Hish = sqrt(H2);
+    }
+
+    if (Hish > 1.0) {
+      epsilon = sqrt(T / Hish);
+    } else {
+      epsilon = sqrt(T);
+    }
     gen_R = std::bind(generate_rotation_perturbation <N_COMP>, std::placeholders::_1, std::placeholders::_2, epsilon, order);
     pert_type = "PERTURB";
   } else {
-- 
cgit v1.2.3-70-g09d2