Got the variation calculation fixed and working.

1 files changed, 21 insertions, 16 deletions

diff --git a/stokes.hpp b/stokes.hpp
index 43a69e3..6a00637 100644
--- a/stokes.hpp
+++ b/stokes.hpp
@@ -21,32 +21,37 @@ Vector<Scalar> variation(const Vector<Scalar>& z, const Vector<Scalar>& z´, con
   Vector<Scalar> ż = zDot(z, dH);
   double ż² = ż.squaredNorm();
 
-  double Reż·z´² = real(pow(ż.dot(z´), 2));
+  double Reż·z´ = real(ż.dot(z´));
 
   Matrix<Scalar> dż = (dH.conjugate() - (dH.dot(z) / z²) * z.conjugate()) * z.adjoint() / z²;
   Matrix<Scalar> dżc = -ddH + (ddH * z.conjugate()) * z.transpose() / z²
-    + (z.dot(dH) / z²) * (Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²);
+    + (z.dot(dH) / z²) * (
+        Matrix<Scalar>::Identity(ddH.rows(), ddH.cols()) - z.conjugate() * z.transpose() / z²
+      );
 
-  Vector<Scalar> dLdz = - 0.5 * (
-        ż.dot(z´) * (dżc * z´) + z´.dot(ż) * (dż * z´.conjugate())
-        - (dż * ż.conjugate() + dżc * ż) * Reż·z´² / ż²
+  Vector<Scalar> dLdz = - Reż·z´ * (
+        dżc * z´ + dż * z´.conjugate() - (dż * ż.conjugate() + dżc * ż) * Reż·z´ / ż²
       ) / (ż² * z´²);
 
   Vector<Scalar> ż´ = -(ddH * z´).conjugate() + ((ddH * z´).dot(z) / z²) * z.conjugate() + (
-        dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²) * z.conjugate()
+        dH.dot(z) * z´.conjugate() + dH.dot(z´) * z.conjugate() - (
+          dH.dot(z) * (z´.dot(z) + z.dot(z´)) / z²
+        ) * z.conjugate()
       ) / z²;
 
-  Vector<Scalar> ddtdLdz´ = -0.5 * (
-      ż´.dot(z´) * ż.conjugate() + ż.dot(z´´) * ż.conjugate() + ż.dot(z´) * ż´.conjugate()
-      - (
-        Reż·z´² * z´´.conjugate() + real(2.0 * ż.dot(z´) * (ż.dot(z´´) + ż´.dot(z´))) * z´.conjugate()
-        - Reż·z´² * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() / z´²
-        ) / z´²
-      - (
-        (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
-        ) * (
-        ż.dot(z´) * ż.conjugate() - Reż·z´² / z´² * z´.conjugate()
+  double dReż·z´ = real(ż.dot(z´´) + ż´.dot(z´));
+
+  Vector<Scalar> ddtdLdz´ = - (
+      Reż·z´ * (
+        ż´.conjugate() - (
+          Reż·z´ * z´´.conjugate() + dReż·z´ * z´.conjugate()
+          - (Reż·z´ / z´²) * (z´´.dot(z´) + z´.dot(z´´)) * z´.conjugate() 
+          ) / z´²
         )
+      + dReż·z´ * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate())
+      - Reż·z´ * (
+        (ż.dot(ż´) + ż´.dot(ż)) / ż² + (z´´.dot(z´) + z´.dot(z´´)) / z´²
+        ) * (ż.conjugate() - Reż·z´ / z´² * z´.conjugate())
       ) / (ż² * z´²);
 
   return dLdz - ddtdLdz´;