flatness of middle earth

"""
In S. Weinberg's book "Gravitation and Cosmology", page 7 there's an exercise where the curvature of Middle earth is estimated by the distances between every pair of 4 points: Umbar, Hobbiton, Erebor and Dagorlad (which we denote as points A,B,C & D). There are overall 6 distinct pairs of points. Euclidean geometry together with the distances of 5 pairs can be used to estimate the 6th. The relative difference between the estimated distance and the actual distance should be the same order of magnitude as the ratio between the actual distance and the radius of curvature. If the estimated distance is greater than the actual distance, then the curvature is positive, otherwise it's negative.
"""

def cosine_angle(os, ns1, ns2):
    """
    Returns the cosine of the angle in a triangle, accoring to the cosine theorem
    os - Length of opposing side
    ns1 - Length of neighboring side
    ns2 - Length of the other neighboring side

    >>> abs(cosine_angle(2.**0.5,1.,1.))<1e-14
    True
    >>> import math; abs(cosine_angle(1.,1.,1.)-math.cos(math.pi/3))<1e-14
    True
    """

    return (ns1**2+ns2**2-os**2)/(2*ns1*ns2)

def main():

    import math

    # Distances between points, in miles
    AB = 1112
    AC = 1498
    AD = 780
    BC = 813
    BD = 960
    CD = 735

    cos_ADB = cosine_angle(AB,AD,BD)
    cos_BDC = cosine_angle(BC,BD,CD)
    sin_ADB = math.sqrt(1-cos_ADB**2)
    sin_BDC = math.sqrt(1-cos_BDC**2)
    cos_ADC = cos_ADB*cos_BDC-sin_ADB*sin_BDC
    estimated_AC = math.sqrt(AD**2 + CD**2 - 2*AD*CD*cos_ADC)

    print('Estimated length of AC[miles]: '+str(estimated_AC))
    print('Actual length of AC[miles]: '+str(AC))
    print('Difference between them [miles]: '+str(AC-estimated_AC))
    print('Relative difference: '+str(100*(AC-estimated_AC)/AC)[:5]+'%')

if __name__ == '__main__':
    import doctest
    doctest.testmod()

    main()