Merge branch 'master' of 10.9.0.2:wolff

2 files changed, 32 insertions, 5 deletions

diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 85d0e11..e1bf33a 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -163,22 +163,36 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const vecto
   double m_dot_v = 0;
 
   for (q_t i = 0; i < q; i++) {
-    m[0][i] = gsl_ran_ugaussian(r);
+    m[0][i] = gsl_ran_ugaussian(r); // create a random vector
     m_dot_v += m[0][i] * v[i];
   }
 
   double v2 = 0;
-  double factor = epsilon * gsl_ran_ugaussian(r);
 
   for (q_t i = 0; i < q; i++) {
-    m[0][i] = m[0][i] - m_dot_v * v[i] + factor * v[i];
+    m[0][i] = m[0][i] - m_dot_v * v[i]; // find the component orthogonal to v
     v2 += pow(m[0][i], 2);
   }
 
   double mag_v = sqrt(v2);
 
   for (q_t i = 0; i < q; i++) {
-    m[0][i] /= mag_v;
+    m[0][i] /= mag_v; // normalize
+  }
+
+  v2 = 0;
+
+  double factor = gsl_ran_gaussian(r, epsilon);
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] += factor * v[i]; // perturb orthogonal vector in original direction
+    v2 += pow(m[0][i], 2);
+  }
+
+  mag_v = sqrt(v2);
+
+  for (q_t i = 0; i < q; i++) {
+    m[0][i] /= mag_v; // normalize
   }
 
   return m;
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 6cd8777..2356875 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -142,8 +142,21 @@ 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);
+    }
+
+    epsilon = sqrt((N_COMP - 1) * T / (D + Hish / 2)) / 2;
+
     gen_R = std::bind(generate_rotation_perturbation <N_COMP>, std::placeholders::_1, std::placeholders::_2, epsilon, order);
-    pert_type = "PERTURB";
+    pert_type = "PERTURB5";
   } else {
     gen_R = generate_rotation_uniform <N_COMP>;
     pert_type = "UNIFORM";