def test_sag():
    rng = np.random.default_rng(57)
    for i in range(100):
        R = 1./rng.normal(0.0, 0.3)
        sphere = batoid.Sphere(R)
        for j in range(10):
            x = rng.uniform(-0.7*abs(R), 0.7*abs(R))
            y = rng.uniform(-0.7*abs(R), 0.7*abs(R))
            result = sphere.sag(x, y)
            np.testing.assert_allclose(
                result,
                R*(1-np.sqrt(1.0-(x*x + y*y)/R/R))
            )
            # Check that it returned a scalar float and not an array
            assert isinstance(result, float)
        # Check 0,0
        np.testing.assert_allclose(sphere.sag(0, 0), 0.0, rtol=0, atol=1e-17)
        # Check vectorization
        x = rng.uniform(-0.7*abs(R), 0.7*abs(R), size=(10, 10))
        y = rng.uniform(-0.7*abs(R), 0.7*abs(R), size=(10, 10))
        np.testing.assert_allclose(
            sphere.sag(x, y),
            R*(1-np.sqrt(1.0-(x*x + y*y)/R/R))
        )
        # Make sure non-unit stride arrays also work
        np.testing.assert_allclose(
            sphere.sag(x[::5,::2], y[::5,::2]),
            R*(1-np.sqrt(1.0-(x*x + y*y)/R/R))[::5,::2]
        )
        do_pickle(sphere)