import numpy as np
import matplotlib
import matplotlib.pyplot as plt

r = 6371000.
r = 100.
h = 2.
y0 = r + h

x = np.linspace(0, 10, 1000)


a = np.sqrt( (1+x/r)**2 - 1.)

xi = ( a * y0 ) / (1 + a**2)
yi = y0 + xi * a

phi = np.arctan( yi / xi )
y = r*( np.pi/2 - phi)

plt.fill_between(x, 0, y, color='g',alpha=0.3)
plt.plot( x, np.sqrt( (yi-y0)**2 + xi**2), 'r-', label='calculation')
plt.plot( x, np.sqrt(2*r*x + x**2), 'b--', label='wikipedia')
plt.xlabel("height [m]")
plt.ylabel("horizon distance [m]")
plt.title("Horizon; distance? (radius = %.3f)" % r)
plt.legend(loc=2)
plt.show()