Calculate distance between two coordinates on a globe
Question: I get two coordinate pairs in the form `90°0′0″N 0°0′0″E` as string and want
to calculate the distance between those points on a sphere with radius
R=6371km.

I found two formulas on the internet [here](http://www.movable-
type.co.uk/scripts/latlong.html), the "haversine" and the "spherical law of
cosines", but they don't seem to work. For a 90° angle which should return
`2*pi*R / 4`, the haversine operates correct but the cosines fail and return
0. A different point with more random coordinates returns false values with
both algorithms: the haversine is too high and the cosines are too low.

Is my implementation wrong or did I chose an incorrect algorithm?

How should I make these calculations (coordinate pairs to distance on globe
surface) instead?

(And yes, I know that I'm not checking for N/S and E/W yet, but the tested
coordinates are all in the north-eastern hemisphere.)

Here's my Python 3 code:

    
    
    import math, re
    R = 6371
    PAT = r'(\d+)°(\d+)′(\d+)″([NSEW])'
    
    def distance(first, second):
        def coords_to_rads(s):  
            return [math.radians(int(d) +int(m)/60 +int(s)/3600) \
                    for d, m, s, nswe in re.findall(PAT, s)]
    
        y1, x1 = coords_to_rads(first)
        y2, x2 = coords_to_rads(second)  
        dx = x1 - x2  
        dy = y1 - y2  
    
        print("coord string:", first, "|", second)
        print("coord radians:", y1, x1, "|", y2, x2)
        print("x/y-distances:", dy, dx)
    
        a = math.sin(dx/2)**2 + math.cos(x1) * math.cos(x2) * math.sin(dy/2)**2  
        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))  
        haversine =  R * c  
    
        law_of_cosines = math.acos( math.sin(x1) * math.sin(x2) + \
                                    math.cos(x1) * math.cos(x2) ) * R
    
        print("HS:", round(haversine, 2), "LOC:", round(law_of_cosines, 2))
    
        return haversine
        #return law_of_cosines
    
    if __name__ == '__main__':
        def test(result, correct):
            print("result: ", result)
            print("correct:", correct)
    
        test(distance("90°0′0″N 0°0′0″E", "0°0′0″N, 0°0′0″E"), 10007.5)
        test(distance("51°28′48″N 0°0′0″E", "46°12′0″N, 6°9′0″E"), 739.2)
        test(distance("90°0′0″N 0°0′0″E", "90°0′0″S, 0°0′0″W"), 20015.1)
        test(distance("33°51′31″S, 151°12′51″E", "40°46′22″N 73°59′3″W"), 15990.2)
    

Here is some output:

    
    
    coord string: 90°0′0″N 0°0′0″E | 0°0′0″N, 0°0′0″E
    coord radians: 1.5707963267948966 0.0 | 0.0 0.0
    x/y-distances: 1.5707963267948966 0.0
    HS: 10007.54 LOC: 0.0
    result: 10007.543398010286
    correct: 10007.5
    
    coord string: 51°28′48″N 0°0′0″E | 46°12′0″N, 6°9′0″E
    coord radians: 0.8984954989266809 0.0 | 0.8063421144213803 0.10733774899765128
    x/y-distances: 0.09215338450530064 -0.10733774899765128
    HS: 900.57 LOC: 683.85
    result: 900.5669567853056
    correct: 739.2
    


Answer: It looks like you mixed up `x` and `y` in your calculation of `a`. You're
supposed to take the cosine of latitude (`y`), not longitude (`x`).

I discovered this by changing your `distance` to `angular_distance` (i.e.
don't multiply by `R`) and adding some additional tests:

    
    
    test(angular_distance("90°0′0″N 0°0′0″E", "89°0′0″N, 0°0′0″E"), math.radians(1))
    test(angular_distance("90°0′0″N 0°0′0″E", "80°0′0″N, 0°0′0″E"), math.radians(10))
    test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 0°0′0″E"), math.radians(40))
    test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 20°0′0″E"), math.radians(40))
    

