Skip to content

Instantly share code, notes, and snippets.

@lingxz lingxz/schwarzschild.py
Last active Jul 6, 2019

Embed
What would you like to do?
# code for generating the graphs in https://theconfused.me/blog/numerical-integration-of-light-paths-in-a-schwarzschild-metric/
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as spi
def schw_null_geodesics(w, t, p):
r, rdot, phi = w
M, L = p
phidot = L / r**2
return [
rdot,
L**2 * (r - 3*M) / r**4,
phidot
]
def expected_schw_light_bending(r, M):
return 4. * M / r
def solve(t, initial, p):
sol = spi.odeint(schw_null_geodesics, initial, t, args=(p,),)
r = sol[:,0]
phi = sol[:,2]
return r, phi
def calculate_deflection(x, y, t):
num_entries = len(t)
gradient = (y[num_entries*4//5] - y[num_entries-1]) / (x[num_entries*4//5] - x[num_entries-1])
return np.arctan(-gradient)
def calculate_initial_conditions(b, initial_x):
# calculate initial conditions for each impact parameter
initial_r = np.sqrt(b**2 + initial_x**2)
initial_phi = np.arccos(initial_x / initial_r)
initial_rdot = np.cos(initial_phi)
initial_phidot = -np.sqrt((1 - initial_rdot**2) / initial_r**2)
L = initial_r**2 * initial_phidot
return [initial_r, initial_rdot, initial_phi], L
def plot_light_rays():
plt.figure()
M = 1. # natural units
initial_x = -50
max_t = 10000.
# impact parameters to plot
bs = np.array([4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40])
# time points to evaluate
t = np.arange(0, max_t, 0.01)
for b in bs:
initial, L = calculate_initial_conditions(b, initial_x)
r, phi = solve(t, initial, [M, L])
x = r * np.cos(phi)
y = r * np.sin(phi)
plt.plot(x, y)
axes = plt.gca()
lim = -initial_x
axes.set_xlim([-lim, lim])
axes.set_ylim([-lim, lim])
# plot the black hole, schwarzschild radius is 2M
circle = plt.Circle((0., 0.), 2*M, color='black', fill=True, zorder=10)
plt.legend(bs)
axes.add_artist(circle)
axes.set_aspect('equal', adjustable='box')
def plot_deflection_angles():
plt.figure()
M = 1. # natural units
initial_x = -50
max_t = 10000.
# impact parameters to plot
bs = np.linspace(20, 100, 40)
# time points to evaluate
t = np.arange(0, max_t, 0.01)
deflections = []
for b in bs:
initial, L = calculate_initial_conditions(b, initial_x)
r, phi = solve(t, initial, [M, L])
x = r * np.cos(phi)
y = r * np.sin(phi)
deflections.append(calculate_deflection(x, y, t))
# plot deflections angles against analytical predicted deflections angles
expected = expected_schw_light_bending(bs, M)
plt.plot(bs, expected, 'ro')
plt.plot(bs, deflections, 'b+')
plt.ylabel('Deflection angle')
plt.xlabel('Distance of closest approach')
plt.legend([r'Theoretical deflection ($\frac{4M}{R}$)', 'Numerical result'])
if __name__ == '__main__':
plt.style.use('ggplot')
plot_light_rays()
plot_deflection_angles()
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.