fixed recording of the spin stiffness

1 files changed, 6 insertions, 6 deletions

diff --git a/src/new_model.cpp b/src/new_model.cpp
index 31e7b5a..4393a3d 100644
--- a/src/new_model.cpp
+++ b/src/new_model.cpp
@@ -52,6 +52,7 @@ typedef struct _xy_bond_props {
 typedef ising::graph<ising_site_props> graph_i;
 typedef xy::graph<std::tuple<>, xy_bond_props> graph_x;
 
+// the sin of the angle between two vectors, signed: vsin(v1, v2) = -vsin(v2, v1)
 double vsin(const vector_t<2, double>& v1, const vector_t<2, double>& v2) {
   return v1[0] * v2[1] - v1[1] * v2[0];
 }
@@ -105,14 +106,15 @@ class x_measurement : public xy::wolff::measurement<orthogonal_t<2, double>, vec
     }
 
     void plain_bond_visited(const xy::wolff::system<orthogonal_t<2, double>, vector_t<2, double>, graph_x>& S, const graph_x::halfedge& e, const vector_t<2, double>& s1_new, double dE) override {
-      cos_xy -= dE;
 
       vector_t<2, double> dsin = (vsin(s1_new, S.s[e.neighbor.ind]) - vsin(S.s[e.self.ind], S.s[e.neighbor.ind])) * e.prop.dNeighborSelf;
 
       if (e.prop.distance == 0) {
+        cos_xy -= dE;
         sin_xy += K1y * (1 + (*si0) * a * *(e.prop.spin)) * dsin;
       } else if (e.prop.distance == 1) {
-        sin_xy += K1yn * dsin;
+        cos_xy -= 2 * dE;
+        sin_xy += K1y * K1yn * dsin;
       }
     }
 };
@@ -329,7 +331,7 @@ int main (int argc, char *argv[]) {
     si.s[pow(L, D) + i].x = true;
   }
 
-  cos_xy = K1 * 2 * pow(L, 2) + K2 * 2 * pow(L, 2);
+  cos_xy = K1 * 2 * pow(L, 2) + K1 * K2 * 4 * pow(L, 2);
   sin_xy.fill(0);
 
   i_measurement mi(si, α, K1);
@@ -343,8 +345,6 @@ int main (int argc, char *argv[]) {
     sxy.run_wolff(1, xy::wolff::generate_rotation_uniform<2, graph_x>, mxy, rng);
     total_cos += cos_xy;
     total_sin2 += sin_xy * sin_xy;
-
-//    std::cout << cos_xy / sxy.ne << " " << sin_xy / pow(sxy.ne, 2) << " " << E << "\n";
   }
 
   std::ifstream checkfile("out.dat");
@@ -358,7 +358,7 @@ int main (int argc, char *argv[]) {
   if (!already_exists) {
     outfile << "N L T J K1 K2 \\[Alpha] M \\[Rho] \\[Sigma]\n";
   }
-  outfile << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2) / N / sxy.ne << " " <<(total_cos + total_sin2) / N / sxy.ne << "\n";
+  outfile << N << " " << L << " " << T << " " << J << " " << K1 << " " << K2 << " " << α << " " << mi.totalA2 / N / pow(si.nv, 2) << " " << (total_cos - total_sin2 / T) / N / sxy.ne << " " <<(total_cos + total_sin2 / T) / N / sxy.ne << "\n";
   outfile.close();
 }