from math import acos, atan2, cos, degrees, pi, radians, sin, sqrt

R = 6378137.0

def distance(pos1, pos2):
    if pos1 == pos2:
        return 0.0

    lat1 = radians(pos1[0])
    lon1 = radians(pos1[1])
    lat2 = radians(pos2[0])
    lon2 = radians(pos2[1])
    
    return acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * R
    
def intermediate_point(pos1, pos2, f):
    if pos1 == pos2:
        return pos1

    lat1 = radians(pos1[0])
    lon1 = radians(pos1[1])
    lat2 = radians(pos2[0])
    lon2 = radians(pos2[1])
    
    delta = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1))
    
    if delta == 0: # too close
        return pos1 if f < 0.5 else pos2

    a = sin((1 - f) * delta) / delta
    b = sin(f * delta) / delta
    x = a * cos(lat1) * cos(lon1) + b * cos(lat2) * cos(lon2)
    y = a * cos(lat1) * sin(lon1) + b * cos(lat2) * sin(lon2)
    z = a * sin(lat1) + b * sin(lat2)
    
    lat3 = atan2(z, sqrt(x**2 + y**2))
    lon3 = atan2(y, x)
    
    def normalize(pos):
        return ((pos[0]+540)%360)-180, ((pos[1]+540)%360)-180
    
    return normalize((degrees(lat3), degrees(lon3)))

@@ -268,20 +268,12 @@
     
-return acos(
+a = 
 sin(
@@ -326,16 +326,77 @@
 n2-lon1)
+%0A    %0A    if a %3E 1:%0A        return 0.0%0A    %0A    return acos(a
 ) * R%0A  
@@ -604,21 +604,12 @@
     
-delta = acos(
+a = 
 sin(
@@ -662,16 +662,112 @@
 n2-lon1)
+%0A    %0A    if a %3E 1: # too close%0A        return pos1 if f %3C 0.5 else pos2%0A    %0A    delta = acos(a
 )%0A    %0A 

Fix math domain error (#1209)