summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--animation.hpp111
-rw-r--r--dimers.cpp122
-rw-r--r--euclidean.hpp52
-rw-r--r--ising.hpp18
-rw-r--r--ising_animate.cpp50
-rw-r--r--spheres.hpp112
-rw-r--r--spheres_infinite.cpp250
-rw-r--r--spin.hpp2
8 files changed, 400 insertions, 317 deletions
diff --git a/animation.hpp b/animation.hpp
new file mode 100644
index 0000000..857354e
--- /dev/null
+++ b/animation.hpp
@@ -0,0 +1,111 @@
+
+#include "space_wolff.hpp"
+#include <GL/glut.h>
+
+template <typename T, int D>
+void drawSphere(const Vector<T, D>& x, Radius r, std::array<double, 3> color = {0, 0, 0},
+ unsigned nPoints = 20) {
+ glColor3d(color[0], color[1], color[2]);
+ glBegin(GL_POLYGON);
+ for (unsigned i = 0; i < nPoints; i++) {
+ glVertex2d(x(0) + r * cos(2 * i * M_PI / nPoints), x(1) + r * sin(2 * i * M_PI / nPoints));
+ }
+ glEnd();
+}
+
+template <typename T, int D>
+void draw(const Spin<T, D, Radius>& s, std::array<double, 3> color = {0, 0, 0}) {
+ drawSphere<T, D>(s.x, s.s, color);
+}
+
+template <typename T, int D>
+void draw(const Spin<T, D, Dimer<T, D>>& s, std::array<double, 3> color = {0, 0, 0}) {
+ drawSphere<T, D>(s.x + s.s.relativePosition, s.s.radius, color);
+ drawSphere<T, D>(s.x - s.s.relativePosition, s.s.radius, color);
+}
+
+template <int D>
+void draw(const Spin<signed, D, IsingSpin>& s, std::array<double, 3> color = {0, 0, 0}) {
+ if (s.s == 1) {
+ glColor3d(color[0], color[1], color[2]);
+ } else if (s.s == -1) {
+ glColor3d(1 - color[0], 1 - color[1], 1 - color[2]);
+ }
+ glRecti(s.x(0), s.x(1), s.x(0) + 1, s.x(1) + 1);
+}
+
+template <typename T, int D, class R, class S> class Animation : public measurement<T, D, R, S> {
+private:
+ unsigned n;
+ unsigned wait;
+ double L;
+ R s₀⁻¹;
+
+public:
+ Animation(T L, unsigned w, int argc, char* argv[], unsigned wait, bool periodic = false)
+ : s₀⁻¹(0), wait(wait), L(1000 * L) {
+ n = 0;
+
+ glutInit(&argc, argv);
+ glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
+ glutInitWindowSize(w, w);
+ glutCreateWindow("wolffWindow");
+ glClearColor(0.0, 0.0, 0.0, 0.0);
+ glMatrixMode(GL_PROJECTION);
+ glLoadIdentity();
+ if (periodic) {
+ gluOrtho2D(0, L, 0, L);
+ } else {
+ gluOrtho2D(-L, L, -L, L);
+ }
+ }
+
+ void pre_cluster(const Model<T, D, R, S>& m, unsigned,
+ const Transformation<T, D, R, S>* t) override {
+ s₀⁻¹ = m.s0.inverse();
+ glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
+ glClear(GL_COLOR_BUFFER_BIT);
+ for (const Spin<T, D, S>* s : m.s) {
+ draw(s₀⁻¹.act(*s));
+ }
+ if (n > wait) {
+ glColor3d(1.0, 0.0, 0.0);
+ glLineWidth(3);
+ glBegin(GL_LINES);
+ Vector<double, D> r = (t->r.t).template cast<double>() / 2.0;
+ double θ = atan2(r(1), r(0));
+ Vector<double, D> v1 = s₀⁻¹.template act<double>({r(0) + L * sin(θ), r(1) - L * cos(θ)});
+ Vector<double, D> v2 = s₀⁻¹.template act<double>({r(0) - L * sin(θ), r(1) + L * cos(θ)});
+ glVertex2d(v1(0), v1(1));
+ glVertex2d(v2(0), v2(1));
+ glEnd();
+ }
+ glFlush();
+ }
+
+ void plain_site_transformed(const Model<T, D, R, S>& m,
+ const Transformation<T, D, R, S>& t) override {
+ if (n > wait) {
+ std::array<double, 3> color;
+ if (t.current().size() > 1) {
+ color = {0, 0, 1};
+ } else {
+ color = {0, 1, 0};
+ }
+ for (const Spin<T, D, S>* s : t.current()) {
+ if (s != NULL) {
+ draw(s₀⁻¹.act(*s), color);
+ } else {
+ }
+ }
+ }
+ }
+
+ void post_cluster(const Model<T, D, R, S>& m) override {
+ if (n > wait) {
+ glFlush();
+ sleep(2);
+ }
+ n++;
+ }
+};
diff --git a/dimers.cpp b/dimers.cpp
new file mode 100644
index 0000000..4416333
--- /dev/null
+++ b/dimers.cpp
@@ -0,0 +1,122 @@
+
+#include <chrono>
+#include <fstream>
+#include <iostream>
+
+#include "animation.hpp"
+#include "space_wolff.hpp"
+#include "spheres.hpp"
+
+const unsigned D = 2;
+typedef Model<double, D, Euclidean<double, D>, Dimer<double, 2>> model;
+
+int main(int argc, char* argv[]) {
+ const unsigned D = 2;
+
+ double L = 32;
+ unsigned N = 1000;
+ double T = 2.0 / log(1.0 + sqrt(2.0));
+ double H = 1.0;
+ unsigned n = 25;
+ unsigned wait = 1000;
+
+ double k = 1e2;
+ double a = 0.1;
+
+ int opt;
+
+ while ((opt = getopt(argc, argv, "n:N:L:T:H:a:k:w:")) != -1) {
+ switch (opt) {
+ case 'n':
+ n = (unsigned)atof(optarg);
+ break;
+ case 'N':
+ N = (unsigned)atof(optarg);
+ break;
+ case 'L':
+ L = atof(optarg);
+ break;
+ case 'T':
+ T = atof(optarg);
+ break;
+ case 'H':
+ H = atof(optarg);
+ break;
+ case 'a':
+ a = atof(optarg);
+ break;
+ case 'k':
+ k = atof(optarg);
+ break;
+ case 'w':
+ wait = atoi(optarg);
+ break;
+ default:
+ exit(1);
+ }
+ }
+
+ auto zSingle = zSpheres<D>(a, k);
+
+ std::function<double(const Spin<double, D, Dimer<double, D>>&,
+ const Spin<double, D, Dimer<double, D>>&)>
+ Z = [zSingle, a, k](const Spin<double, D, Dimer<double, D>>& s1,
+ const Spin<double, D, Dimer<double, D>>& s2) -> double {
+ Spin<double, D, Radius> s11 = {.x = s1.x + s1.s.relativePosition, .s = s1.s.radius};
+ Spin<double, D, Radius> s12 = {.x = s1.x - s1.s.relativePosition, .s = s1.s.radius};
+ Spin<double, D, Radius> s21 = {.x = s2.x + s2.s.relativePosition, .s = s2.s.radius};
+ Spin<double, D, Radius> s22 = {.x = s2.x - s2.s.relativePosition, .s = s2.s.radius};
+
+ return zSingle(s11, s21) + zSingle(s12, s21) + zSingle(s11, s22) + zSingle(s12, s22);
+ };
+
+ auto g1 = nudgeGen<D, Dimer<double, D>>(1);
+ auto g2 = swapGen<D, Dimer<double, D>>(0.1);
+ auto g3 = accrossGen<D, Dimer<double, D>>(0.1);
+ auto g4 = centerGen<D, Dimer<double, D>>(0);
+
+ auto tag = std::chrono::high_resolution_clock::now();
+
+ std::string filename = "flips_" + std::to_string(n) + "_" + std::to_string(T) + "_" +
+ std::to_string(H) + "_" + std::to_string(a) + "_" + std::to_string(k) +
+ "_" + std::to_string(tag.time_since_epoch().count()) + ".dat";
+
+ std::ofstream file;
+ file.open(filename);
+
+ Animation<double, D, Euclidean<double, D>, Dimer<double, D>> A(L, 750, argc, argv, wait);
+ model sphere(1.0, Z, bCenter<D, Dimer<double, D>>(H));
+
+ Rng rng;
+
+ sphere.s.resize(n);
+
+ unsigned nx = floor(sqrt(n));
+ for (unsigned i = 0; i < sphere.s.size(); i++) {
+ Spin<double, 2, Dimer<double, D>>* ss = new Spin<double, 2, Dimer<double, D>>();
+ ss->x = {(i / nx) * L / nx - L / 2, (i % nx) * L / nx - L / 2};
+ ss->s = Dimer<double, D>{.relativePosition = {0.2, 0}, .radius = 0.25};
+ sphere.s[i] = ss;
+ sphere.dict.insert(ss);
+ }
+
+ sphere.wolff(T, {g1, g2, g3, g4}, A, N);
+
+ file.close();
+
+ /*
+ std::ofstream snapfile;
+ snapfile.open("sphere_snap.dat");
+
+ for (Spin<double, D, double>* s : sphere.s) {
+ Spin<double, D, double> rs = sphere.s0.inverse().act(*s);
+ snapfile << rs.s << " " << rs.x.transpose() << "\n";
+ delete s;
+ }
+
+ snapfile.close();
+ */
+
+ return 0;
+}
+
diff --git a/euclidean.hpp b/euclidean.hpp
index b1ce2c6..fb2dbcd 100644
--- a/euclidean.hpp
+++ b/euclidean.hpp
@@ -5,6 +5,15 @@
#include "spin.hpp"
#include "vector.hpp"
+typedef double Radius;
+typedef signed IsingSpin;
+
+template <class T, int D> class Dimer {
+public:
+ double radius;
+ Vector<T, D> relativePosition;
+};
+
template <class T, int D> class Euclidean {
public:
Vector<T, D> t;
@@ -25,36 +34,27 @@ public:
r = r0;
}
- Vector<T, D> act(const Vector<T, D>& x) const {
- return t + r * x;
- }
-
- template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
- Spin<T, D, S> s_new;
-
- s_new.x = this->act(s.x);
- s_new.s = s.s;
+ Vector<T, D> act(const Vector<T, D>& x) const { return t + r * x; }
- return s_new;
+ template <typename U> Vector<U, D> act(const Vector<U, D>& x) const {
+ return t.template cast<U>() + r.template cast<U>() * x;
}
- Euclidean act(const Euclidean& x) const {
- Vector<T, D> tnew = r * x.t + t;
- Matrix<T, D> rnew = r * x.r;
+ Radius act(Radius r) const { return r; }
- Euclidean pnew(tnew, rnew);
+ IsingSpin act(IsingSpin s) const { return s; }
- return pnew;
+ Dimer<T, D> act(const Dimer<T, D>& d) const {
+ return {.radius = d.radius, .relativePosition = r * d.relativePosition};
}
- Euclidean inverse() const {
- Vector<T, D> tnew = -r.transpose() * t;
- Matrix<T, D> rnew = r.transpose();
+ template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
+ return {.x = act(s.x), .s = act(s.s)};
+ }
- Euclidean pnew(tnew, rnew);
+ Euclidean act(const Euclidean& x) const { return Euclidean(r * x.t + t, r * x.r); }
- return pnew;
- }
+ Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); }
};
template <class T, int D> class TorusGroup {
@@ -105,6 +105,16 @@ public:
return s_new;
}
+ template <typename U> Vector<U, D> act(const Vector<U, D>& s) const {
+ Vector<U, D> s_new = t.template cast<U>() + r.template cast<U>() * s;
+
+ for (unsigned i = 0; i < D; i++) {
+ s_new(i) = fmod(L + s_new(i), L);
+ }
+
+ return s_new;
+ }
+
TorusGroup act(const TorusGroup& x) const {
Vector<T, D> tnew = r * x.t + t;
Matrix<T, D> rnew = r * x.r;
diff --git a/ising.hpp b/ising.hpp
index 6f140be..a97da26 100644
--- a/ising.hpp
+++ b/ising.hpp
@@ -91,15 +91,15 @@ void isingPopulate(isingModel& m, unsigned L, double pop, double mag) {
for (unsigned i = 0; i < L; i++) {
for (unsigned j = 0; j < L; j++) {
if (rng.uniform<double>(0, 1) < pop) {
- Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
- ss->x = {i, j};
- if (rng.uniform<double>(0, 1) < mag) {
- ss->s = 1;
- } else {
- ss->s = -1;
- }
- m.s.push_back(ss);
- m.dict.insert(ss);
+ Spin<signed, D, signed>* ss = new Spin<signed, D, signed>();
+ ss->x = {i, j};
+ if (rng.uniform<double>(0, 1) < mag) {
+ ss->s = 1;
+ } else {
+ ss->s = -1;
+ }
+ m.s.push_back(ss);
+ m.dict.insert(ss);
}
}
}
diff --git a/ising_animate.cpp b/ising_animate.cpp
index cd65b6b..f3ab314 100644
--- a/ising_animate.cpp
+++ b/ising_animate.cpp
@@ -1,45 +1,5 @@
+#include "animation.hpp"
#include "ising.hpp"
-#include <GL/glut.h>
-
-class Animation : public measurement<signed, D, TorusGroup<signed, D>, signed> {
-private:
- bool color;
-
-public:
- Animation(double L, unsigned w, bool tcolor, int argc, char* argv[]) {
- color = tcolor;
-
- glutInit(&argc, argv);
- glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
- glutInitWindowSize(w, w);
- glutCreateWindow("wolffWindow");
- glClearColor(0.0, 0.0, 0.0, 0.0);
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- gluOrtho2D(0, L, 0, L);
- }
-
- void post_cluster(const isingModel& m) override {
- glClearColor(1.0, 1.0, 1.0, 1.0);
- glClear(GL_COLOR_BUFFER_BIT);
- for (const Spin<signed, 2, signed>* s : m.s) {
- if (s->s == 1) {
- if (color)
- glColor3f(1.0, 0.0, 0.0);
- else
- glColor3f(1.0, 1.0, 1.0);
- } else if (s->s == -1) {
- if (color)
- glColor3f(0.0, 0.0, 1.0);
- else
- glColor3f(0.0, 0.0, 0.0);
- }
- Vector<signed, 2> xx = m.s0.inverse().act(*s).x;
- glRecti(xx(0), xx(1), xx(0) + 1, xx(1) + 1);
- }
- glFlush();
- }
-};
int main(int argc, char* argv[]) {
unsigned L = 32;
@@ -50,10 +10,11 @@ int main(int argc, char* argv[]) {
double T = 2.0 / log(1.0 + sqrt(2.0));
double H = 1.0;
bool color = false;
+ unsigned wait = N;
int opt;
- while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:c")) != -1) {
+ while ((opt = getopt(argc, argv, "N:L:T:H:m:r:p:cw:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned)atof(optarg);
@@ -79,6 +40,9 @@ int main(int argc, char* argv[]) {
case 'c':
color = true;
break;
+ case 'w':
+ wait = atoi(optarg);
+ break;
default:
exit(1);
}
@@ -97,7 +61,7 @@ int main(int argc, char* argv[]) {
auto g = isingGen(L);
- Animation A(L, 750, color, argc, argv);
+ Animation<signed, D, TorusGroup<signed, D>, IsingSpin> A(L, 750, argc, argv, wait, true);
ising.wolff(T, {g}, A, N);
diff --git a/spheres.hpp b/spheres.hpp
new file mode 100644
index 0000000..040313c
--- /dev/null
+++ b/spheres.hpp
@@ -0,0 +1,112 @@
+
+#include "space_wolff.hpp"
+
+template <int D>
+std::function<double(const Spin<double, D, Radius>&, const Spin<double, D, Radius>&)>
+zSpheres(double a, double k) {
+ return [a, k](const Spin<double, D, Radius>& s1, const Spin<double, D, Radius>& s2) -> double {
+ Vector<double, D> d = s1.x - s2.x;
+
+ double σ = s1.s + s2.s;
+ double δ = σ - sqrt(d.transpose() * d);
+
+ if (δ > -a * σ) {
+ return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
+ } else if (δ > -2 * a * σ) {
+ return 0.5 * k * pow(δ + 2 * a * σ, 2);
+ } else {
+ return 0;
+ }
+ };
+}
+
+template <int D, class S> std::function<double(Spin<double, D, S>)> bCenter(double H) {
+ return [H](Spin<double, D, S> s) -> double { return H * s.x.norm(); };
+}
+
+template <int D> Euclidean<double, D> flipAround(double θ, Vector<double, D>& t₀) {
+ Matrix<double, D> m;
+
+ m(0, 0) = -cos(2 * θ);
+ m(1, 1) = cos(2 * θ);
+ m(0, 1) = -2 * cos(θ) * sin(θ);
+ m(1, 0) = -2 * cos(θ) * sin(θ);
+
+ return Euclidean<double, D>(t₀ - m * t₀, m);
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> nudgeGen(double ε) {
+ return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+ Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+ double θ = rng.uniform((double)0.0, 2 * M_PI);
+
+ Spin<double, D, S>* s = rng.pick(M.s);
+ Vector<double, D> t = s->x;
+ for (unsigned j = 0; j < D; j++) {
+ t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
+ }
+
+ return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s);
+ };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> swapGen(double ε) {
+ return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+ Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+ Spin<double, D, S>* s1 = rng.pick(M.s);
+ Spin<double, D, S>* s2 = rng.pick(M.s);
+
+ while (s1 == s2) {
+ s2 = rng.pick(M.s);
+ }
+
+ Vector<double, D> t1 = s1->x;
+ Vector<double, D> t2 = s2->x;
+ Vector<double, D> t12 = t1 - t2;
+ Vector<double, D> t = (t1 + t2) / 2;
+
+ double θ =
+ atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
+
+ return new PairFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s1, s2);
+ };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> accrossGen(double ε) {
+ return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+ Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+ Spin<double, D, S>* s1 = rng.pick(M.s);
+ Spin<double, D, S>* s2 = rng.pick(M.s);
+
+ while (s1 == s2) {
+ s2 = rng.pick(M.s);
+ }
+
+ Vector<double, D> t1 = s1->x;
+ Vector<double, D> t2 = s2->x;
+ Vector<double, D> t12 = t1 - t2;
+ Vector<double, D> t = t2;
+
+ double θ =
+ atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
+
+ return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s1);
+ };
+}
+
+template <int D, class S> Gen<double, D, Euclidean<double, D>, S> centerGen(double ε) {
+ return [ε](Model<double, D, Euclidean<double, D>, S>& M,
+ Rng& rng) -> Transformation<double, D, Euclidean<double, D>, S>* {
+ double θ = rng.uniform((double)0.0, 2 * M_PI);
+
+ Vector<double, D> t = M.s0.t;
+ for (unsigned j = 0; j < D; j++) {
+ t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
+ }
+
+ Spin<double, D, S>* s = rng.pick(M.s);
+
+ return new SpinFlip<double, D, Euclidean<double, D>, S>(M, flipAround(θ, t), s);
+ };
+}
+
diff --git a/spheres_infinite.cpp b/spheres_infinite.cpp
index 56abe77..af91427 100644
--- a/spheres_infinite.cpp
+++ b/spheres_infinite.cpp
@@ -4,226 +4,12 @@
#include <chrono>
#include "space_wolff.hpp"
-#include <GL/glut.h>
+#include "spheres.hpp"
+#include "animation.hpp"
const unsigned D = 2;
typedef Model<double, D, Euclidean<double, D>, double> model;
-class animation : public measurement<double, 2, Euclidean<double, D>, double> {
-private:
- uint64_t t1;
- uint64_t t2;
- unsigned n;
- unsigned tmp;
- unsigned wait;
- double L;
- Euclidean<double, D> s0_tmp;
- std::ofstream& file;
-
-public:
- animation(double L, unsigned w, int argc, char* argv[], unsigned wait, std::ofstream& file) : s0_tmp(0), wait(wait), L(50 * L), file(file) {
- t1 = 0;
- t2 = 0;
- n = 0;
-
- glutInit(&argc, argv);
- glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
- glutInitWindowSize(w, w);
- glutCreateWindow("wolffWindow");
- glClearColor(0.0, 0.0, 0.0, 0.0);
- glMatrixMode(GL_PROJECTION);
- glLoadIdentity();
- gluOrtho2D(-L, L, -L, L);
- }
-
- void pre_cluster(const model& m, unsigned, const Transformation<double, D, Euclidean<double, D>, double>* t) override {
- s0_tmp = m.s0.inverse();
- tmp = 0;
- glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
- glClear(GL_COLOR_BUFFER_BIT);
- glColor3d(1.0, 0.0, 0.0);
- glBegin(GL_LINES);
- Vector<double, 2> r = t->r.t / 2;
- double θ = atan2(r(1), r(0));
- Vector<double, 2> v1 = s0_tmp.act({r(0) + L * sin(θ), r(1) - L * cos(θ)});
- Vector<double, 2> v2 = s0_tmp.act({r(0) - L * sin(θ), r(1) + L * cos(θ)});
- glVertex2d(v1(0), v1(1));
- glVertex2d(v2(0), v2(1));
- if (n > wait)
- file << v1.transpose() << " " << v2.transpose() << std::endl;
- glEnd();
- for (const Spin<double, 2, double>* s : m.s) {
- glBegin(GL_POLYGON);
- unsigned n_points = 50;
- glColor3f(0.0f, 0.0f, 0.0f);
- if (n > wait)
- file << s << " " << s->s << " " << s0_tmp.act(*s).x.transpose() << " ";
- for (unsigned i = 0; i < n_points; i++) {
- glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
- s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
- }
- glEnd();
- }
- if (n > wait)
- file << std::endl;
- glLineWidth(3);
- glFlush();
- }
-
- void plain_site_transformed(const model& m, const Transformation<double, D, Euclidean<double, D>, double>& t) override {
- if (n > wait) {
- if (t.current().size() > 1) {
- glColor3f(0.0f, 0.0f, 1.0f);
- } else {
- glColor3f(0.0f, 1.0f, 0.0f);
- }
- for (const Spin<double, 2, double>* s : t.current()) {
- if (s != NULL) {
- file << s << " " << s0_tmp.act(t.r.act(*s)).x.transpose() << " ";
- glBegin(GL_POLYGON);
- unsigned n_points = 50;
- for (unsigned i = 0; i < n_points; i++) {
- glVertex2d(s0_tmp.act(*s).x(0) + s->s * cos(2 * i * M_PI / n_points),
- s0_tmp.act(*s).x(1) + s->s * sin(2 * i * M_PI / n_points));
- }
- glEnd();
- } else {
- file << s << " " << s0_tmp.act({0, 0}).transpose() << " ";
- }
- }
- }
- tmp++;
- }
-
- void post_cluster(const model& m) override {
- t1 += tmp;
- t2 += tmp * tmp;
- if (n > wait) {
- file << std::endl;
- glFlush();
- sleep(2);
- }
- n++;
- }
-
- void clear() {
- t1 = 0;
- t2 = 0;
- n = 0;
- }
-
- double var() { return (t2 - t1 * t1 / (double)n) / (double)n; }
-};
-
-Gen<double, D, Euclidean<double, D>, double> eGen(double ε) {
- return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
- Vector<double, D> t;
- Matrix<double, D> m;
-
- double θ = rng.uniform((double)0.0, 2 * M_PI);
- m(0, 0) = -cos(2 * θ);
- m(1, 1) = cos(2 * θ);
- m(0, 1) = -2 * cos(θ) * sin(θ);
- m(1, 0) = -2 * cos(θ) * sin(θ);
-
- Spin<double, D, double>* s = rng.pick(M.s);
- t = s->x;
- for (unsigned j = 0; j < D; j++) {
- t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
- }
-
- Euclidean<double, D> g(t - m * t, m);
- return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
- };
-}
-
-Gen<double, D, Euclidean<double, D>, double> mGen(double ε) {
- return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
- Matrix<double, D> m;
-
- Spin<double, D, double>* s1 = rng.pick(M.s);
- Spin<double, D, double>* s2 = rng.pick(M.s);
-
- while (s1 == s2) {
- s2 = rng.pick(M.s);
- }
-
- Vector<double, D> t1 = s1->x;
- Vector<double, D> t2 = s2->x;
- Vector<double, D> t12 = t1 - t2;
- Vector<double, D> t = (t1 + t2) / 2;
-
- double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
-
- m(0, 0) = -cos(2 * θ);
- m(1, 1) = cos(2 * θ);
- m(0, 1) = -2 * cos(θ) * sin(θ);
- m(1, 0) = -2 * cos(θ) * sin(θ);
-
- Vector<double, D> t3 = t - m * t;
-
- if (t3(0) != t3(0)) {
- std::cout << t3 << "\n" << t << "\n" << m << "\n" << t12 << "\n" ;
- getchar();
- }
-
- Euclidean<double, D> g(t3, m);
- return new PairFlip<double, D, Euclidean<double, D>, double>(M, g, s1, s2);
- };
-}
-
-Gen<double, D, Euclidean<double, D>, double> tGen(double ε) {
- return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
- Matrix<double, D> m;
-
- Spin<double, D, double>* s1 = rng.pick(M.s);
- Spin<double, D, double>* s2 = rng.pick(M.s);
-
- while (s1 == s2) {
- s2 = rng.pick(M.s);
- }
-
- Vector<double, D> t1 = s1->x;
- Vector<double, D> t2 = s2->x;
- Vector<double, D> t12 = t1 - t2;
- Vector<double, D> t = t2;
-
- double θ = atan2(t12(1), t12(0)) + rng.variate<double, std::normal_distribution>(0.0, ε) / t12.norm();
-
- m(0, 0) = -cos(2 * θ);
- m(1, 1) = cos(2 * θ);
- m(0, 1) = -2 * cos(θ) * sin(θ);
- m(1, 0) = -2 * cos(θ) * sin(θ);
-
- Vector<double, D> t3 = t - m * t;
-
- Euclidean<double, D> g(t3, m);
- return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s1);
- };
-}
-
-Gen<double, D, Euclidean<double, D>, double> rGen(double ε) {
- return [ε](model& M, Rng& rng) -> Transformation<double, D, Euclidean<double, D>, double>* {
- Vector<double, D> t;
- Matrix<double, D> m;
-
- double θ = rng.uniform((double)0.0, 2 * M_PI);
- m(0, 0) = -cos(2 * θ);
- m(1, 1) = cos(2 * θ);
- m(0, 1) = -2 * cos(θ) * sin(θ);
- m(1, 0) = -2 * cos(θ) * sin(θ);
-
- Spin<double, D, double>* s = rng.pick(M.s);
- t = M.s0.t;
- for (unsigned j = 0; j < D; j++) {
- t(j) += rng.variate<double, std::normal_distribution>(0.0, ε);
- }
-
- Euclidean<double, D> g(t - m * t, m);
- return new SpinFlip<double, D, Euclidean<double, D>, double>(M, g, s);
- };
-}
-
int main(int argc, char* argv[]) {
const unsigned D = 2;
@@ -270,26 +56,6 @@ int main(int argc, char* argv[]) {
}
}
- std::function<double(const Spin<double, D, double>&, const Spin<double, D, double>&)> Z =
- [L, a, k](const Spin<double, D, double>& s1, const Spin<double, D, double>& s2) -> double {
- Vector<double, D> d = s1.x - s2.x;
-
- double σ = s1.s + s2.s;
- double δ = σ - sqrt(d.transpose() * d);
-
- if (δ > -a * σ) {
- return 0.5 * k * (2 * pow(a * σ, 2) - pow(δ, 2));
- } else if (δ > -2 * a * σ) {
- return 0.5 * k * pow(δ + 2 * a * σ, 2);
- } else {
- return 0;
- }
- };
-
- std::function<double(Spin<double, D, double>)> B = [L, H](Spin<double, D, double> s) -> double {
- return H * s.x.norm();
- };
-
std::function<double(Spin<double, D, double>)> B_hard = [L, H](Spin<double, D, double> s) -> double {
if (fabs(s.x(0)) < 3 * L / 4 && fabs(s.x(1)) < 3 * L / 4) {
return 0;
@@ -298,10 +64,10 @@ int main(int argc, char* argv[]) {
}
};
- auto g1 = eGen(1);
- auto g2 = mGen(0.1);
- auto g3 = tGen(0.1);
- auto g4 = rGen(0);
+ auto g1 = nudgeGen<D, Radius>(1);
+ auto g2 = swapGen<D, Radius>(0.1);
+ auto g3 = accrossGen<D, Radius>(0.1);
+ auto g4 = centerGen<D, Radius>(0);
auto tag = std::chrono::high_resolution_clock::now();
@@ -311,8 +77,8 @@ int main(int argc, char* argv[]) {
std::ofstream file;
file.open(filename);
- animation A(L, 750, argc, argv, wait, file);
- model sphere(1.0, Z, B);
+ Animation<double, D, Euclidean<double, D>, Radius> A(L, 750, argc, argv, wait);
+ model sphere(1.0, zSpheres<D>(a, k), bCenter<D, Radius>(H));
Rng rng;
diff --git a/spin.hpp b/spin.hpp
index bb8369b..d4d39cc 100644
--- a/spin.hpp
+++ b/spin.hpp
@@ -7,7 +7,5 @@ template <class T, int D, class S> class Spin {
public:
Vector<T, D> x;
S s;
-
- Spin() { x = Vector<T, D>::Zero(); }
};