Moved some repeated code into the pSpinModel class.

1 files changed, 7 insertions, 16 deletions

diff --git a/langevin.cpp b/langevin.cpp
index 8ecd69c..0e5a616 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -81,17 +81,10 @@ int main(int argc, char* argv[]) {
     }
   }
 
-  Complex κ(Rκ, Iκ);
-  Real σp = sqrt(factorial(p) / ((Real)2 * pow(N, p - 1)));
-  Real σ2 = sqrt(factorial(2) / ((Real)2 * pow(N, 2 - 1)));
-  Real σ4 = sqrt(factorial(4) / ((Real)2 * pow(N, 4 - 1)));
-
   Rng r;
 
-  pSpinModel<Real, 2, 4> pSpin;
-
-  std::get<0>(pSpin.Js) = generateCouplings<Real, 2>(N, std::normal_distribution<Real>(0, σ2), r.engine());
-  std::get<1>(pSpin.Js) = generateCouplings<Real, 4>(N, std::normal_distribution<Real>(0, σ4 / 10), r.engine());
+//  pSpinModel<Real, 2, 4> pSpin({J2, J4});
+  pSpinModel<Real, 2, 4>pSpin(N, r.engine(), 1, 0.1);
 
   std::normal_distribution<Real> Red(0, 1);
 
@@ -100,18 +93,17 @@ int main(int argc, char* argv[]) {
   RealVector dHr;
   RealMatrix ddHr;
   std::tie(Hr, dHr, ddHr, std::ignore) = pSpin.hamGradHess(zMin);
-  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + 2.0 * zMin * zMin.transpose());
-  std::cout << eigenS.eigenvalues().transpose() << std::endl;
+  Eigen::EigenSolver<Matrix<Real>> eigenS(ddHr - (dHr * zMin.transpose() + zMin.dot(dHr) * Matrix<Real>::Identity(zMin.size(), zMin.size()) + (ddHr * zMin) * zMin.transpose()) / (Real)zMin.size() + zMin * zMin.transpose() / (Real)N);
   for (unsigned i = 0; i < N; i++) {
-    RealVector zNew = normalize(zMin + 0.01 * eigenS.eigenvectors().col(i).real());
-    std::cout << pSpin.getHamiltonian(zNew) - Hr << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << std::endl;
+    RealVector zNew = normalize(zMin + 1e-3 * eigenS.eigenvectors().col(i).real());
+    std::cout << eigenS.eigenvalues()(i) << " " << 2.0 * (pSpin.getHamiltonian(zNew) - Hr) / 1e-6 << " " << real(eigenS.eigenvectors().col(i).dot(zMin)) << " " << eigenS.eigenvectors().col(0).dot(eigenS.eigenvectors().col(i)) << std::endl;
   }
   std::cout << std::endl;
   getchar();
 
   complex_normal_distribution<Real> d(0, 1, 0);
 
-  pSpinModel<Complex, 2, 4> complexPSpin = pSpin.cast<Complex>();;
+  pSpinModel complexPSpin = pSpin.cast<Complex>();;
   ComplexVector zSaddle = zMin.cast<Complex>();
 
   ComplexVector zSaddleNext;
@@ -134,8 +126,7 @@ int main(int argc, char* argv[]) {
   Real φ = atan2( H2.imag() - H1.imag(), H1.real() - H2.real());
   std::cerr << (zSaddle - zSaddleNext).norm() / (Real)N << " " << φ << " " << H1 * exp(Complex(0, φ)) << " " << H2 * exp(Complex(0, φ)) << std::endl;
 
-  std::get<0>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<0>(complexPSpin.Js);
-  std::get<1>(complexPSpin.Js) = exp(Complex(0, φ)) * std::get<1>(complexPSpin.Js);
+  complexPSpin *= exp(Complex(0, φ));
 
   Cord test(complexPSpin, zSaddle, zSaddleNext, 3);
   test.relaxNewton(10, 1, 1e4);