Work on Stokes lines.

1 files changed, 15 insertions, 5 deletions

diff --git a/stokes.hpp b/stokes.hpp
index 3814cf4..6cc5c4a 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -1,4 +1,5 @@
 #include "p-spin.hpp"
+#include <iostream>
 
 class ropeRelaxationStallException: public std::exception {
   virtual const char* what() const throw() {
@@ -130,6 +131,11 @@ class Rope {
         δz[i] = variation(z[i], dz(i), ddz(i), dH, ddH);
       }
 
+      for (unsigned i = 1; i < z.size() - 1; i++) {
+        δz[i] = δz[i].conjugate() - (δz[i].dot(z[i]) / z[i].squaredNorm()) * z[i].conjugate();
+//        δz[i] = δz[i] - ((δz[i].conjugate().dot(dz(i))) / dz(i).squaredNorm()) * dz(i).conjugate();
+      }
+
       // We return a δz with average norm of one.
       double mag = 0;
       for (unsigned i = 1; i < z.size() - 1; i++) {
@@ -152,14 +158,18 @@ class Rope {
       // We rescale the step size by the distance between points.
       double Δl = length() / (n() + 1);
 
-      while (rNew.cost(J) >= cost(J)) {
+      while (true) {
         for (unsigned i = 1; i < z.size() - 1; i++) {
-          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i].conjugate());
+          rNew.z[i] = normalize(z[i] - (δ * Δl) * Δz[i]);
         }
 
-        δ /= 2;
+        if (rNew.cost(J) < cost(J)) {
+          break;
+        } else {
+          δ /= 2;
+        }
 
-        if (δ < 1e-30) {
+        if (δ < 1e-15) {
           throw ropeRelaxationStallException();
         }
       }
@@ -206,7 +216,7 @@ class Rope {
       double δ = δ0;
       try {
         for (unsigned i = 0; i < N; i++) {
-          δ = 2 * step(J, δ);
+          δ = 1.1 * step(J, δ);
         }
       } catch (std::exception& e) {
       }