from kerrgeodesic_gw import KerrBH
a = var('a')
M = KerrBH(a)
M.boyer_lindquist_coordinates()
init_point = M((0, 6, pi/2, 0), name='p0')
geod = M.geodesic([0, 300], init_point, mu=1, E=0.883, \
L=1.982, Q=0.467, a_num=0.998)
geod.integrate(step=0.005)
graph = geod.plot(label_axes=True)
graph._extra_kwds['aspect_ratio'] = 1
sphinx_plot(graph)