diff --git a/src/simsopt/geo/surface.py b/src/simsopt/geo/surface.py
index 5cd9cf32a..2382013a5 100644
--- a/src/simsopt/geo/surface.py
+++ b/src/simsopt/geo/surface.py
@@ -438,7 +438,8 @@ def is_self_intersecting(self, angle=0., thetas=None):
         r"""
         This function computes a cross section of self at the input cylindrical angle.  Then,
         approximating the cross section as a piecewise linear polyline, the Bentley-Ottmann algorithm
-        is used to check for self-intersecting segments of the cross section.
+        is used to check for self-intersecting segments of the cross section.  NOTE: if this function returns False,
+        this is not a guarantee that the surface is everywhere not self-intersecting.
 
         Args:
             angle: the cylindrical angle at which we would like to check whether the surface is self-intersecting.  Note that a
diff --git a/tests/geo/test_surface.py b/tests/geo/test_surface.py
index a71ad3132..8d95b4d7e 100755
--- a/tests/geo/test_surface.py
+++ b/tests/geo/test_surface.py
@@ -16,7 +16,7 @@
     best_nphi_over_ntheta
 from simsopt.geo.curverzfourier import CurveRZFourier
 from simsopt._core.json import GSONDecoder, GSONEncoder, SIMSON
-from .surface_test_helpers import get_surface
+from .surface_test_helpers import get_surface, get_boozer_surface
 
 TEST_DIR = (Path(__file__).parent / ".." / "test_files").resolve()
 
@@ -388,6 +388,21 @@ def test_is_self_intersecting(self):
         
         s = get_surface('SurfaceRZFourier', True, full=True, nphi=200, ntheta=200, mpol=2, ntor=2)
         assert not s.is_self_intersecting()
+        
+        # make sure it works on an NCSX BoozerSurface
+        bs, boozer_surf = get_boozer_surface()
+        s = boozer_surf.surface
+        assert not s.is_self_intersecting(angle=0.123*np.pi)
+        assert not s.is_self_intersecting(angle=0.123*np.pi, thetas=200)
+        assert not s.is_self_intersecting(thetas=231)
+        
+        # make sure it works on a perturbed NCSX BoozerSurface
+        dofs = s.x.copy()
+        dofs[14]+=0.2
+        s.x = dofs
+        assert s.is_self_intersecting(angle=0.123*np.pi)
+        assert s.is_self_intersecting(angle=0.123*np.pi, thetas=200)
+        assert s.is_self_intersecting(thetas=202) 
 
 
 class UtilTests(unittest.TestCase):