def test_intersect():
    rng = np.random.default_rng(5772)
    size = 10_000
    for i in range(100):
        R = 1./rng.normal(0.0, 0.3)
        sphereCoordSys = batoid.CoordSys(origin=[0, 0, -1])
        sphere = batoid.Sphere(R)
        x = rng.uniform(-0.3*abs(R), 0.3*abs(R), size=size)
        y = rng.uniform(-0.3*abs(R), 0.3*abs(R), size=size)
        z = np.full_like(x, -2*abs(R))
        # If we shoot rays straight up, then it's easy to predict the intersection
        vx = np.zeros_like(x)
        vy = np.zeros_like(x)
        vz = np.ones_like(x)
        rv = batoid.RayVector(x, y, z, vx, vy, vz)
        np.testing.assert_allclose(rv.z, -2*abs(R))
        rv2 = batoid.intersect(sphere, rv.copy(), sphereCoordSys)
        assert rv2.coordSys == sphereCoordSys

        rv2 = rv2.toCoordSys(batoid.CoordSys())
        np.testing.assert_allclose(rv2.x, x)
        np.testing.assert_allclose(rv2.y, y)
        np.testing.assert_allclose(rv2.z, sphere.sag(x, y)-1, rtol=0, atol=1e-9)

        # Check default intersect coordTransform
        rv2 = rv.copy().toCoordSys(sphereCoordSys)
        batoid.intersect(sphere, rv2)
        assert rv2.coordSys == sphereCoordSys
        rv2 = rv2.toCoordSys(batoid.CoordSys())
        np.testing.assert_allclose(rv2.x, x)
        np.testing.assert_allclose(rv2.y, y)
        np.testing.assert_allclose(rv2.z, sphere.sag(x, y)-1, rtol=0, atol=1e-9)