bugfix in Pmap for multiple s or zeta positions.

bugfix in Pmap (solve th+lamdba-th2=0),th2 was not correctly set if multiple s or zeta positions were evaluated. Extended the domain tests to check all combinations of 1d arrays of s,th,ze with multiple / single element

src/gvec_to_python/geometry/domain.py

@@ -1025,9 +1025,9 @@ class GVEC_domain:
             for i in range(s.shape[0]):
                 for j in range(th2.shape[1]):
                     for k in range(ze.shape[2]):
-                        func = lambda th: th + self.LA_base.eval_stz(s[i, 0, 0], th, ze[0, 0, k], self.LA_coef) - th2[0, j, 0]
+                        func = lambda th: th + self.LA_base.eval_stz(s[i, 0, 0], th, ze[0, 0, k], self.LA_coef) - th2[i, j, k]
                         fprime = lambda th: 1 + self.LA_base.eval_stz(s[i, 0, 0], th, ze[0, 0, k], self.LA_coef, der='th')
-                        th[i, j, k] = newton(func, th2[0, j, 0], fprime)
+                        th[i, j, k] = newton(func, th2[i, j, k], fprime)
 
         return s + 0*th, th, ze + 0*th