Found two bugs in current calc

src/gvec_to_python/GVEC_functions.py

@@ -449,9 +449,6 @@ class GVEC:
                 Second entry are the Cartesian coordinates of evaluation points.
         '''
         j = self.jv(s, a1, a2)
-        print(np.max(np.abs(j[0])))
-        print(np.max(np.abs(j[1])))
-        print(np.max(np.abs(j[2])))
         return self.transform(self.jv(s, a1, a2), s, a1, a2, kind='push_vector'), self.f(s, a1, a2)
 
     # test equilibirum condition grad p = J x B / mu0
@@ -474,9 +471,9 @@ class GVEC:
         print('b1:', np.max(np.abs(b[0])))
         print('b2:', np.max(np.abs(b[1])))
         print('b3:', np.max(np.abs(b[2])))
-        print('jxb, 1st comp, rel error:', np.max(np.abs(p_s - (j[1] * b[2] - j[2] * b[1])) / (abs_j * abs_b)))
-        print('jxb, 2nd comp, rel error:', np.max(np.abs(j[2] * b[0] - j[0] * b[2]) / (abs_j * abs_b)))
-        print('jxb, 3rd comp, rel error:', np.max(np.abs(j[0] * b[1] - j[1] * b[0]) / (abs_j * abs_b)))
+        print('p_s = (j2 x bv)_1, rel error:', np.max(np.abs(p_s - (j[1] * b[2] - j[2] * b[1])) / (abs_j * abs_b)))
+        print('0   = (j2 x bv)_2, rel error:', np.max(np.abs(j[2] * b[0] - j[0] * b[2]) / (abs_j * abs_b)))
+        print('0   = (j2 x bv)_3, rel error:', np.max(np.abs(j[0] * b[1] - j[1] * b[0]) / (abs_j * abs_b)))
         #assert np.allclose(p_s, j[1] * b[2] - j[2] * b[1])
         #assert np.allclose(0*b[0], j[2] * b[0] - j[0] * b[2])
         #assert np.allclose(0*b[0], j[0] * b[1] - j[1] * b[0])
@@ -531,16 +528,16 @@ class GVEC:
         -------
             A 3-tuple of either floats (single-point evaluation) or 3d numpy arrays (meshgrid evaluation).'''
 
-        mu0 = 1.#1.2566370621219e-6
+        mu0 = 1.2566370621219e-6
         b = np.ascontiguousarray(self._bv_gvec(s, th, ze))
 
         grad_bt = self._grad_b_theta(s, th, ze)
         grad_bz = self._grad_b_zeta(s, th, ze)
-        db_ds = np.ascontiguousarray((0*grad_bt[0], grad_bt[0], grad_bz[0]))
-        db_dt = np.ascontiguousarray((0*grad_bt[1], grad_bt[1], grad_bz[1]))
-        db_dz = np.ascontiguousarray((0*grad_bt[2], grad_bt[2], grad_bz[2]))
+        db_ds = np.ascontiguousarray((0.*grad_bt[0], grad_bt[0], grad_bz[0]))
+        db_dt = np.ascontiguousarray((0.*grad_bt[1], grad_bt[1], grad_bz[1]))
+        db_dz = np.ascontiguousarray((0.*grad_bt[2], grad_bt[2], grad_bz[2]))
 
-        g = self.domain.dg_gvec(s, th, ze, der=None)
+        g = self.domain.g_gvec(s, th, ze)
         dg_ds = self.domain.dg_gvec(s, th, ze, der='s')
         dg_dt = self.domain.dg_gvec(s, th, ze, der='th')
         dg_dz = self.domain.dg_gvec(s, th, ze, der='ze')

src/gvec_to_python/geometry/domain.py

@@ -320,7 +320,7 @@ class GVEC_domain:
 
         q1, q2, ze2 = self.Xmap(s, th, ze)
 
-        return q1 * dq1_di * self.det_dX(s, th, ze) + self.det_dh(q1, q2, ze2) * self.dJ_X(s, th, ze)[comp]
+        return dq1_di * self.det_dX(s, th, ze) + self.det_dh(q1, q2, ze2) * self.dJ_X(s, th, ze)[comp]
 
     # gvec_pest coordinates
     def f_pest(self, s2, th2, ze2):
@@ -857,11 +857,11 @@ class GVEC_domain:
             # meshgrid evaluation
             if dX.ndim == 5:
                 tmp += [ddX[0][:, :, :, 0, n] * dX2_dt + ddX[1][:, :, :, 1, n] * dX1_ds
-                    - ddX[0][:, :, :, 1, n] * dX1_dt + ddX[1][:, :, :, 0, n] * dX2_ds]
+                    - ddX[0][:, :, :, 1, n] * dX1_dt - ddX[1][:, :, :, 0, n] * dX2_ds]
             # single point evaluation
             else:
                 tmp += [ddX[0][0, n] * dX2_dt + ddX[1][1, n] * dX1_ds
-                    - ddX[0][1, n] * dX1_dt + ddX[1][0, n] * dX2_ds]
+                    - ddX[0][1, n] * dX1_dt - ddX[1][0, n] * dX2_ds]
 
         return tmp[0], tmp[1], tmp[2]