From 58e62d50243bc45dbd4ae466642e801a1b46dd69 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 9 Nov 2021 17:59:53 +0100
Subject: Progress towards a systematized approach.

---
 collectStokesData.hpp | 52 +++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)
 create mode 100644 collectStokesData.hpp

(limited to 'collectStokesData.hpp')

diff --git a/collectStokesData.hpp b/collectStokesData.hpp
new file mode 100644
index 0000000..50aad24
--- /dev/null
+++ b/collectStokesData.hpp
@@ -0,0 +1,52 @@
+#include "stokes.hpp"
+#include <iostream>
+
+using Complex = std::complex<Real>;
+
+template<int ...ps, class Generator, typename... T>
+void collectStokesData(unsigned N, Generator& r, double ε, T... μs) {
+  pSpinModel<Real, ps...> ReM(N, r, μs...);
+
+  std::normal_distribution<Real> Red(0, 1);
+
+  Vector<Real> xMin = randomMinimum(ReM, Red, r, ε);
+
+  Real Hx;
+  Vector<Real> dHx;
+  Matrix<Real> ddHx;
+  std::tie(Hx, dHx, ddHx, std::ignore) = ReM.hamGradHess(xMin);
+  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHx - xMin.dot(dHx) * Matrix<Real>::Identity(xMin.size(), xMin.size()) / Real(N));
+
+  complex_normal_distribution<Real> d(0, 1, 0);
+
+  pSpinModel M = ReM.template cast<Complex>();;
+  Vector<Complex> zMin = xMin.cast<Complex>();
+
+  Vector<Complex> zSaddle = zMin;
+
+  while ((zSaddle - zMin).norm() < 1e-3 * N || abs(imag(M.getHamiltonian(zSaddle))) < 1e-10) {
+    Vector<Complex> z0 = normalize(zSaddle + 0.5 * randomVector<Complex>(N, d, r));
+    zSaddle = findSaddle(M, z0, ε);
+  }
+
+  Complex H1 = M.getHamiltonian(zMin);
+  Complex H2 = M.getHamiltonian(zSaddle);
+
+  Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
+
+  M *= exp(Complex(0, φ));
+
+  Cord c(M, zMin, zSaddle, 3);
+  c.relaxNewton(10, 1, 1e4);
+
+  ((std::cout << ps), ...) << std::endl;;
+  ((std::cout << μs), ...) << std::endl;;
+
+  std::apply([](const Tensor<Real, ps>&... Js) -> void {
+        ((std::cout << Js << std::endl), ...);
+      } , ReM.Js);
+
+  std::cout << xMin.transpose() << std::endl;
+  std::cout << Hx << std::endl;
+  std::cout << eigenS.eigenvalues().real().transpose() << std::endl;
+}
-- 
cgit v1.2.3-54-g00ecf