Python вычислить точку пересечения двух больших кругов

Question

Я пытаюсь вычислить точку пересечения (lat и lon в градусах) двух больших кругов, каждая из которых определяется двумя точками на круге. Я пытался следовать методу, указанному here. Но ответ, который я получаю, неверен, мой код ниже, кто-нибудь видит, где я ошибся?

################################################ 
#### Intersection of two great circles. 
# Points on great circle 1. 
glat1 = 54.8639587 
glon1 = -8.177818 

glat2 = 52.65297082 
glon2 = -10.78064876 

# Points on great circle 2. 
cglat1 = 51.5641564 
cglon1 = -9.2754284 

cglat2 = 53.35422063 
cglon2 = -12.5767799 

# 1. Put in polar coords. 

x1 = cos(glat1) * sin(glon1) 
y1 = cos(glat1) * cos(glon1) 
z1 = sin(glat1) 

x2 = cos(glat2) * sin(glon2) 
y2 = cos(glat2) * cos(glon2) 
z2 = sin(glat2) 


cx1 = cos(cglat1) * sin(cglon1) 
cy1 = cos(cglat1) * cos(cglon1) 
cz1 = sin(cglat1) 

cx2 = cos(cglat2) * sin(cglon2) 
cy2 = cos(cglat2) * cos(cglon2) 
cz2 = sin(cglat2) 


# 2. Get normal to planes containing great circles. 
# It's the cross product of vector to each point from the origin. 

N1 = cross([x1, y1, z1], [x2, y2, z2]) 
N2 = cross([cx1, cy1, cz1], [cx2, cy2, cz2]) 


# 3. Find line of intersection between two planes. 
# It is normal to the poles of each plane. 

L = cross(N1, N2) 


# 4. Find intersection points. 

X1 = L/abs(L) 
X2 = -X1 


ilat = asin(X1[2]) * 180./np.pi 
ilon = atan2(X1[1], X1[0]) * 180./np.pi 

Следует также упомянуть об этом на поверхности Земли (предполагая сферу).

Answer 1

Решение от DSM в комментариях выше, ваши углы находятся в градусах, а sin и cos ожидают радианы. Также линия

X1 = L/abs(L) 

должно быть,

X1 = L/np.sqrt(L[0]**2 + L[1]**2 + L[2]**2) 

Answer 2

еще одно исправление, которое должно быть сделано, чтобы изменить соз/грех перед «долготу» в х и у размеры:

x = cos(lat) * cos(lon) 
y = cos(lat) * sin(lon) 
z = sin(lat) 

Это связано с тем, что оригинальное преобразование из угла в сферическую систему выполняется с использованием полярных/азимутальных углов сферы и они не совпадают с углами лат/лон (wiki it https://en.wikipedia.org/wiki/Spherical_coordinate_system).