def test_normal():
    rng = np.random.default_rng(577)
    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.normal(x, y)
            r = np.hypot(x, y)
            rat = r/R
            dzdr = rat/np.sqrt(1-rat*rat)
            nz = 1/np.sqrt(1+dzdr*dzdr)
            normal = np.array([-x/r*dzdr*nz, -y/r*dzdr*nz, nz])
            np.testing.assert_allclose(result, normal)
        # Check 0,0
        np.testing.assert_equal(sphere.normal(0, 0), np.array([0, 0, 1]))
        # 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))
        r = np.hypot(x, y)
        rat = r/R
        dzdr = rat/np.sqrt(1-rat*rat)
        nz = 1/np.sqrt(1+dzdr*dzdr)
        normal = np.dstack([-x/r*dzdr*nz, -y/r*dzdr*nz, nz])
        np.testing.assert_allclose(
            sphere.normal(x, y),
            normal
        )
        # Make sure non-unit stride arrays also work
        np.testing.assert_allclose(
            sphere.normal(x[::5,::2], y[::5,::2]),
            normal[::5, ::2]
        )