# Angle between 2 arc of circle on a sphere [closed]

I have 4 points with (lat, long) coordinates. Lets say

A = (lat_A, lon_A)

B = (lat_B, lon_B)

C = (lat_C, lon_C)

D = (lat_D, lon_D)

A and B make the first arc, C and D the second. I would like to compute the angle between these 2 arcs.

I can see that these arcs make a great circle on my sphere and calculate the angle between these arcs would be equivalent to calculate the angle between the 2 circles.

I also tried using median point and azimuth but the result doesnt seem to be correct.

Here is my code :

``````import math as m

def midpoint(lat1, lon1, lat2, lon2):
# median point between 2 points on a sphere
mid_lat = m.radians((lat1 + lat2) / 2)
mid_lon = m.radians((lon1 + lon2) / 2)
return mid_lat, mid_lon

def azimuth(lat1, lon1, lat2, lon2):
# compute azimut between 2 points on the sphere
delta_lon = lon2 - lon1
azimuth = m.atan2(m.sin(delta_lon) * m.cos(lat2),
m.cos(lat1)*m.sin(lat2)-m.sin(lat1)*m.cos(lat2)*m.cos(delta_lon))
return azimuth

def angle_between_arcs(A, B, C, D):
# median point
center_AB = midpoint(A[0], A[1], B[0], B[1])
center_CD = midpoint(C[0], C[1], D[0], D[1])

# azimuth
azimuth_A = azimuth(center_AB[0], center_AB[1], A[0], A[1])
azimuth_B = azimuth(center_AB[0], center_AB[1], B[0], B[1])

azimuth_C = azimuth(center_CD[0], center_CD[1], C[0], C[1])
azimuth_D = azimuth(center_CD[0], center_CD[1], D[0], D[1])

# difference
angle_diff = m.degrees(abs(azimuth_pg1 - azimuth_po1) - abs(azimuth_pg2 - azimuth_po2))

return angle_diff
``````

EDIT

I finally found my answer and i'm pretty sure it will be usefull to somebody one day.

First you have to understand that and arc of circle given by two points on the globe define a great circle with the center of the earth as center.

Then you have two great circles and you want to get the angle between these circle. These two circles intersect each other in two points on the globe, and you want to get the angle of the closest point.

Here is the python code :

``````def sphericalToCartesian(longitude, latitude):

# to cartesian
return x, y, z

def cartesianToSpherical(x, y, z):
lat = math.asin(z) * 180./np.pi
lon = math.atan2(y, x) * 180./np.pi
return lon, lat

def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance in kilometers between two
points
on the earth (specified in decimal degrees)
"""
lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2,
lat2])

# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) *
math.sin(dlon/2)**2
c = 2 * math.asin(math.sqrt(a))
r = 6371
return c * r

def getIntersectionPoints(point1_cercle1, point2_cercle1,
point1_cercle2, point2_cercle2):

# each pair allows to create a great circle around the globe
# to cartesian first
p1c1x, p1c1y, p1c1z = sphericalToCartesian(*point1_cercle1)
p2c1x, p2c1y, p2c1z = sphericalToCartesian(*point2_cercle1)

p1c2x, p1c2y, p1c2z = sphericalToCartesian(*point1_cercle2)
p2c2x, p2c2y, p2c2z = sphericalToCartesian(*point2_cercle2)

# normal of each plan corresponding to a great circle
norm1 = np.cross([p1c1x, p1c1y, p1c1z], [p2c1x, p2c1y, p2c1z])
norm2 = np.cross([p1c2x, p1c2y, p1c2z], [p2c2x, p2c2y, p2c2z])

# intersection of normals
line = np.cross(norm1, norm2)

# two points opposed on earth
p0 = line / np.sqrt(line[0]**2 + line[1]**2 + line[2]**2)
p1 = -p0

# to (lon, lat)
i_long0, i_lat0 = cartesianToSpherical(*p0)
i_long1, i_lat1 = cartesianToSpherical(*p1)

# lets keep the closest point
dist_i0, dist_i1 = haversine(i_long0, i_lat0, *point1_cercle1),
haversine(i_long1, i_lat1, *point1_cercle1)
if dist_i0 < dist_i1:
i_point = (i_long0, i_lat0)
else:
i_point = (i_long1, i_lat1)
return i_point

def getAngleInSphericalTriangle(point_0, point_1, point_2):
"""
point_0 is a point of your first great circle or arc of
circle
point_1 is the intersection point
point_2 is a point of your 2nd great circle or arc of
circle
"""
# to cartesian
cart_p0 = sphericalToCartesian(*point_0)
cart_p1 = sphericalToCartesian(*point_1)
cart_p2 = sphericalToCartesian(*point_2)

# angle between points and earth center
p1p2 = np.arccos(np.dot(cart_p1, cart_p2))
p0p1 = np.arccos(np.dot(cart_p0, cart_p1))
p2p0 = np.arccos(np.dot(cart_p2, cart_p0))

#  third fundamental law of spherical trigonometry