diff options
| -rw-r--r-- | lib/orthogonal.h | 22 | ||||
| -rw-r--r-- | src/wolff_On.cpp | 15 |
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"; |