added modulated field to wolff_planar

2 files changed, 55 insertions, 20 deletions

diff --git a/src/wolff_heisenberg.cpp b/src/wolff_heisenberg.cpp
index ad3f3c6..e05453f 100644
--- a/src/wolff_heisenberg.cpp
+++ b/src/wolff_heisenberg.cpp
@@ -4,7 +4,7 @@
 #include <correlation.h>
 #include <wolff.h>
 
-typedef state_t <orthogonal_t <3, double>, vector_t <3, double>> sim_t;
+typedef state_t <orthogonal_t <3, double>, vector_t <3, double>> heisenberg_t;
 
 int main(int argc, char *argv[]) {
 
@@ -63,7 +63,7 @@ int main(int argc, char *argv[]) {
     timestamp = spec.tv_sec*1000000000LL + spec.tv_nsec;
   }
 
-  std::function <orthogonal_t <3, double>(gsl_rng *, const sim_t *)> gen_R;
+  std::function <orthogonal_t <3, double>(gsl_rng *, const heisenberg_t *)> gen_R;
 
   const char *pert_type;
 
@@ -110,20 +110,20 @@ int main(int argc, char *argv[]) {
   free(filename_S);
   free(filename_X);
 
-  std::function <void(const sim_t *)> *measurements = (std::function <void(const sim_t *)> *)malloc(4 * sizeof(std::function <void(const sim_t *)>));
+  std::function <void(const heisenberg_t *)> *measurements = (std::function <void(const heisenberg_t *)> *)malloc(4 * sizeof(std::function <void(const heisenberg_t *)>));
 
-  measurements[0] = [&](const sim_t *s) {
+  measurements[0] = [&](const heisenberg_t *s) {
     float smaller_E = (float)s->E;
     fwrite(&smaller_E, sizeof(float), 1, outfile_E);
   };
-  measurements[1] = [&](const sim_t *s) {
+  measurements[1] = [&](const heisenberg_t *s) {
     float smaller_X = (float)correlation_length(s);
     fwrite(&smaller_X, sizeof(float), 1, outfile_X);
   };
-  measurements[2] = [&](const sim_t *s) {
+  measurements[2] = [&](const heisenberg_t *s) {
     write_magnetization(s->M, outfile_M);
   };
-  measurements[3] = [&](const sim_t *s) {
+  measurements[3] = [&](const heisenberg_t *s) {
     fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S);
   };
 
diff --git a/src/wolff_planar.cpp b/src/wolff_planar.cpp
index bfb382f..4b9b5f0 100644
--- a/src/wolff_planar.cpp
+++ b/src/wolff_planar.cpp
@@ -4,7 +4,27 @@
 #include <wolff.h>
 #include <correlation.h>
 
-typedef state_t <orthogonal_t <2, double>, vector_t <2, double>> sim_t;
+typedef state_t <orthogonal_t <2, double>, vector_t <2, double>> planar_t;
+
+// angle from the x-axis of a two-vector
+double theta(vector_t <2, double> v) {
+  double x = v.x[0];
+  double y = v.x[1];
+
+  double val = atan(y / x);
+
+  if (x < 0.0 && y > 0.0) {
+    return M_PI + val;
+  } else if ( x < 0.0 && y < 0.0 ) {
+    return - M_PI + val;
+  } else {
+    return val;
+  }
+}
+
+double H_modulated(vector_t <2, double> v, int order, double mag) {
+  return mag * cos(order * theta(v));
+}
 
 int main(int argc, char *argv[]) {
 
@@ -13,17 +33,20 @@ int main(int argc, char *argv[]) {
   D_t D = 2;
   L_t L = 128;
   double T = 2.26918531421;
-  double *H = (double *)calloc(MAX_Q, sizeof(double));
+  double *H_vec = (double *)calloc(MAX_Q, sizeof(double));
 
   bool silent = false;
   bool use_pert = false;
 
+  bool modulated_field = false;
+  int order = 2;
+
   int opt;
   q_t J_ind = 0;
   q_t H_ind = 0;
   double epsilon = 1;
 
-  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:q:D:L:T:J:H:spe:mo:")) != -1) {
     switch (opt) {
     case 'N': // number of steps
       N = (count_t)atof(optarg);
@@ -38,7 +61,7 @@ int main(int argc, char *argv[]) {
       T = atof(optarg);
       break;
     case 'H': // external field. nth call couples to state n
-      H[H_ind] = atof(optarg);
+      H_vec[H_ind] = atof(optarg);
       H_ind++;
       break;
     case 's': // don't print anything during simulation. speeds up slightly
@@ -50,6 +73,12 @@ int main(int argc, char *argv[]) {
     case 'e':
       epsilon = atof(optarg);
       break;
+    case 'm':
+      modulated_field = true;
+      break;
+    case 'o':
+      order = atoi(optarg);
+      break;
     default:
       exit(EXIT_FAILURE);
     }
@@ -65,7 +94,7 @@ int main(int argc, char *argv[]) {
 
   const char *pert_type;
 
-  std::function <orthogonal_t <2, double>(gsl_rng *, const sim_t *)> gen_R;
+  std::function <orthogonal_t <2, double>(gsl_rng *, const planar_t *)> gen_R;
 
   if (use_pert) {
     gen_R = std::bind(generate_rotation_perturbation <2>, std::placeholders::_1, std::placeholders::_2, epsilon);
@@ -81,7 +110,7 @@ int main(int argc, char *argv[]) {
   fprintf(outfile_info, "<| \"ID\" -> %lu, \"MODEL\" -> \"PLANAR\", \"q\" -> 2, \"D\" -> %" PRID ", \"L\" -> %" PRIL ", \"NV\" -> %" PRIv ", \"NE\" -> %" PRIv ", \"T\" -> %.15f, \"H\" -> {", timestamp, D, L, L * L, D * L * L, T);
 
   for (q_t i = 0; i < 2; i++) {
-    fprintf(outfile_info, "%.15f", H[i]);
+    fprintf(outfile_info, "%.15f", H_vec[i]);
     if (i < 2 - 1) {
       fprintf(outfile_info, ", ");
     }
@@ -111,26 +140,32 @@ int main(int argc, char *argv[]) {
   free(filename_S);
   free(filename_X);
 
-  std::function <void(const sim_t *)> *measurements = (std::function <void(const sim_t *)> *)calloc(4, sizeof(std::function <void(const sim_t *)>));
+  std::function <void(const planar_t *)> *measurements = (std::function <void(const planar_t *)> *)calloc(4, sizeof(std::function <void(const planar_t *)>));
 
-  measurements[0] = (std::function <void(const sim_t *)>)[&](const sim_t *s) {
+  measurements[0] = (std::function <void(const planar_t *)>)[&](const planar_t *s) {
     float smaller_E = (float)s->E;
     fwrite(&smaller_E, sizeof(float), 1, outfile_E);
   };
-  measurements[1] = [&](const sim_t *s) {
+  measurements[1] = [&](const planar_t *s) {
     float smaller_X = (float)correlation_length(s);
     fwrite(&smaller_X, sizeof(float), 1, outfile_X);
   };
-  measurements[2] = [&](const sim_t *s) {
+  measurements[2] = [&](const planar_t *s) {
     write_magnetization(s->M, outfile_M);
   };
-  measurements[3] = [&](const sim_t *s) {
+  measurements[3] = [&](const planar_t *s) {
     fwrite(&(s->last_cluster_size), sizeof(uint32_t), 1, outfile_S);
   };
 
+  std::function <double(vector_t <2, double>)> H;
 
+  if (modulated_field) {
+    H = std::bind(H_modulated, std::placeholders::_1, order, H_vec[0]);
+  } else {
+    H = std::bind(H_vector <2, double>, std::placeholders::_1, H_vec);
+  }
 
-  wolff <orthogonal_t <2, double>, vector_t <2, double>> (N, D, L, T, dot <2, double>, std::bind(H_vector <2, double>, std::placeholders::_1, H), gen_R, 4, measurements, silent);
+  wolff <orthogonal_t <2, double>, vector_t <2, double>> (N, D, L, T, dot <2, double>, H, gen_R, 4, measurements, silent);
 
   free(measurements);
 
@@ -139,7 +174,7 @@ int main(int argc, char *argv[]) {
   fclose(outfile_S);
   fclose(outfile_X);
 
-  free(H);
+  free(H_vec);
 
   fftw_cleanup();