simd doesn't work but otherwise things are fine

examples/2_Intermediate/boozerQA.py

@@ -114,7 +114,7 @@ def fun(dofs):
     outstr += f", Len=sum([{cl_string}])={sum(J.J().real for J in ls):.1f}"
     outstr += f", ║∇J║={np.linalg.norm(grad.real):.1e}"
     print(outstr)
-    return J.real, grad.real
+    return J, grad
 
 
 print("""
@@ -128,6 +128,16 @@ def fun(dofs):
 h = np.random.uniform(size=dofs.shape)
 J0, dJ0 = f(dofs)
 
+#eps = 1e-10
+#idx = 10
+#dJ_exact = dJ0[idx].real
+#x0 = dofs.copy().astype(complex)
+#x0[idx]+=eps*1j
+#J1, _ = f(x0)
+#print(dJ_exact, np.imag(J1)/eps, (dJ_exact-np.imag(J1)/eps)/dJ_exact)
+#import ipdb;ipdb.set_trace()
+
+
 dJh = sum(dJ0 * h)
 for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8, 1e-9]:
     J1, _ = f(dofs + 2*eps*h)

src/simsopt/geo/boozersurface.py

@@ -775,7 +775,7 @@ def solve_residual_equation_exactly_newton(self, tol=1e-10, maxiter=10, iota=0.,
             else:
                 b = np.concatenate((r[mask], [(label.J()-self.targetlabel), s.gamma()[0, 0, 2]]))
             norm = np.linalg.norm(b.real)
-            if norm <= tol:
+            if norm <= tol and i > 2:
                 break
             if s.stellsym:
                 J = np.vstack((