import numpy as np
from einsteinpy.geodesic import Timelike
from einsteinpy.plotting.geodesic import StaticGeodesicPlotter
Note that, we are working in M-Units ($G = c = M = 1$). Also, since the Schwarzschild spacetime has spherical symmetry, the values of the angular components do not affect the end result (We can always rotate our coordinate system to bring the geodesic in the equatorial plane). Hence, we set $\theta = \pi / 2$ (equatorial plane), with initial $p_\theta = 0$, which implies, that the geodesic should stay in the equatorial plane.
position = [40., np.pi / 2, 0.]
momentum = [0., 0., 3.83405]
a = 0.
steps = 5500
delta = 1.
geod = Timelike(
metric="Schwarzschild",
metric_params=(a,),
position=position,
momentum=momentum,
steps=steps,
delta=delta,
return_cartesian=True
)
sgpl = StaticGeodesicPlotter()
sgpl.plot(geod, color="green", title="3D Geodesic Plot", aspect="auto")
sgpl.show()
sgpl.clear() # In Interactive mode, `clear()` must be called before drawing another plot, to avoid overlap
sgpl.plot(geod, color="green", title="3D Geodesic Plot", aspect="equal")
sgpl.show()
Apsidal Precession is easily observed in the plots above, and as expected, the geodesic is confined to the equatorial plane. We can visualize this better, with a few 2D plots.
sgpl.clear()
sgpl.plot2D(geod, coordinates=(1, 2), color="blue", title="2D Geodesic Plot - XY") # Plot X & Y
sgpl.show()
sgpl.clear()
sgpl.plot2D(geod, coordinates=(1, 3), color="blue", title="2D Geodesic Plot - XZ") # Plot X & Z
sgpl.show()
sgpl.clear()
sgpl.plot2D(geod, coordinates=(2, 3), color="blue", title="2D Geodesic Plot - YZ") # Plot Y & Z
sgpl.show()
Let's see, how the coordinates change, with $\lambda$ (Affine Parameter).
sgpl.clear()
sgpl.parametric_plot(geod, colors=("cyan", "magenta", "lime"), title="2D Geodesic Plot - XYZ vs Lambda") # Plot X, Y, Z vs Lambda (Affine Parameter)
sgpl.show()
The lag between $X1$ ($x$-component) and $X2$ ($y$-component), in the parametric plot, reaffirms the results from above. Also, $X3 \approx 0$ ($z$-component) in this plot, which is expected, as the geodesic never leaves the equatorial plane.