Merge pull request #36 from bendudson/bugfix-core-mask

Change psi_bndry test in critical.core_mask

freegs/critical.py

@@ -265,7 +265,7 @@ def core_mask(R, Z, psi, opoint, xpoint=[], psi_bndry=None):
 
     # Start and end points
     Ro, Zo, psi_axis = opoint[0]
-    if not psi_bndry:
+    if psi_bndry is None:
         _, _, psi_bndry = xpoint[0]
 
     # Normalise psi

freegs/test_critical.py

@@ -68,3 +68,35 @@ def psi_func(R,Z):
     assert len(opoints) == 2
     assert np.isclose(xpoints[0][0], r0, atol = 1./nx)
     assert np.isclose(xpoints[0][1], z0, atol = 1./ny)
+
+
+def test_mask_zero_psi_bndry():
+    nx = 65
+    ny = 65
+
+    r1d = np.linspace(1.0, 2.0, nx)
+    z1d = np.linspace(-1.0, 1.0, nx)
+    r2d, z2d = np.meshgrid(r1d, z1d, indexing='ij')
+
+    r0 = 1.5
+    z0 = 0.0
+
+    # This has one O-point at (r0,z0) and no x-points
+    # Range from around -0.5 to +0.5
+    def psi_func(R,Z):
+        return np.exp(-((R - r0)**2 + (Z - z0)**2)/0.3**2) - 0.5
+
+    psi = psi_func(r2d, z2d)
+
+    opoints, xpoints = critical.find_critical(r2d, z2d, psi)
+
+    assert len(xpoints) == 0
+    assert len(opoints) == 1
+    assert np.isclose(opoints[0][0], r0, atol = 1./nx)
+    assert np.isclose(opoints[0][1], z0, atol = 1./ny)
+
+    mask = critical.core_mask(r2d, z2d, psi, opoints, xpoints, psi_bndry = 0.0)
+
+    # Some of the mask must equal 1, some 0
+    assert np.any(np.equal(mask, 1))
+    assert np.any(np.equal(mask, 0))