How to interpolate circle arc with lines on Earth surface?

I have the coordinates (lat, lon) of 3 points (P0, P1, P2) on the Earth surface. The distance between P1 and P0 is the same as the distance between P2 and P0. Therefore, one can make a circle arc from P1 to P2 with P0 as a centre. My problem is that I need to interpolate this circle arc by N lines to be able to work with it in the GeoJSON format (which does not support circles/circle arcs). Each line segment should ideally have the same legth.

Is there any library suitable for this problem? I have found some code for transformation of a circle to a polygon (such as polycircles), but none of them seemed to support circle arcs.

A simple code is here:

import geopy.distance

# DME OKL: 50 05 45.12 N 14 15 56.19 E
lon0 = 14.265608
lat0 = 50.095867
coords0 = (lat0, lon0)

# 1st circle point: 50 03 10,23 N  014 28 30,47 E
lon1 = 14.475131
lat1 = 50.052842
coords1 = (lat1, lon1)

# 2nd circle point: 49 59 33,96 N  014 24 58,76 E
lon2 = 14.416322
lat2 = 49.992767
coords2 = (lat2, lon2)

# verification of dinstances:
print(geopy.distance.distance(coords1, coords0).nautical)
print(geopy.distance.distance(coords2, coords0).nautical)

# How to interpolate circle arc between points 1 and 2?

I'm not aware of any libraries to do it, so I did it geometrically using vector geometry.

The routine arc in the code below will return two numpy arrays, lon[] and lat[] which contain the longitudes and latitudes of N+1 points along the required circular arc. (The first and last points are P1 and P2).

Method (with C the centre of the sphere): convert to cartesians; find the angle between planes C.P1.P0 and C.P2.P0 (NOT between lines P0.P1 and P0.P2 as I mistakenly first tried); rotate C.P1 about C.P0 by successive increments of angle/N; convert back to longitude and latitude.

Your original points are too close together on a sphere to be interesting, so I've made up some more points to try as well.

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

PI = 4 * math.atan( 1.0 )

def lonlat_to_cartesian( R, lon, lat ):
    lon_rad = lon * PI / 180
    lat_rad = lat * PI / 180
    x = R * math.cos( lat_rad ) * math.cos( lon_rad )
    y = R * math.cos( lat_rad ) * math.sin( lon_rad )
    z = R * math.sin( lat_rad )
    return x, y, z

def cartesian_to_lonlat( x, y, z ):
    R = math.sqrt( x ** 2 + y ** 2 + z ** 2 )
    lon = math.atan2( y, x ) * 180 / PI
    lat = math.asin( z / R ) * 180 / PI
    return R, lon, lat

def rotationMatrix( a, theta ):        # matrix for rotation theta about axis a
    R = np.zeros( ( 3, 3 ) )
    n = a / np.linalg.norm( a )
    C = math.cos( theta )
    S = math.sin( theta )
    R[0,0] = C + n[0] * n[0] * ( 1 - C )
    R[0,1] =     n[0] * n[1] * ( 1 - C ) - n[2] * S
    R[0,2] =     n[0] * n[2] * ( 1 - C ) + n[1] * S
    R[1,1] = C + n[1] * n[1] * ( 1 - C )
    R[1,2] =     n[1] * n[2] * ( 1 - C ) - n[0] * S
    R[1,0] =     n[1] * n[0] * ( 1 - C ) + n[2] * S
    R[2,2] = C + n[2] * n[2] * ( 1 - C )
    R[2,0] =     n[2] * n[0] * ( 1 - C ) - n[1] * S
    R[2,1] =     n[2] * n[1] * ( 1 - C ) + n[0] * S
    return R

def plotSphere( radius ):
    u, v = np.mgrid[0:2*math.pi:50j, 0:math.pi:50j]
    x = radius * np.cos(u) * np.sin(v)
    y = radius * np.sin(u) * np.sin(v)
    z = radius * np.cos(v)
    ax.plot_surface( x, y, z, alpha=0.2 )

def arc( p0_lon, p0_lat, p1_lon, p1_lat, p2_lon, p2_lat, N ):
    p0 = np.array( lonlat_to_cartesian( 1.0, p0_lon, p0_lat ) )      # cartesian coordinates of points
    p1 = np.array( lonlat_to_cartesian( 1.0, p1_lon, p1_lat ) )
    p2 = np.array( lonlat_to_cartesian( 1.0, p2_lon, p2_lat ) )

    n1 = np.cross( p1, p0 );   n1 = n1 / np.linalg.norm( n1 )        # unit normal to plane C.P1.P0
    n2 = np.cross( p2, p0 );   n2 = n2 / np.linalg.norm( n2 )        # unit normal to plane C.P2.P0

    alpha = math.acos( np.dot( n1, n2 ) )                            # angle between planes
    if np.dot( np.cross( p1 - p0, p2 - p0 ), p0 ) < 0: alpha = -alpha

    dalpha = alpha / N                                                                   # angular increment
    derror = ( np.linalg.norm( p2 - p0 ) - np.linalg.norm( p1 - p0 ) ) / N               # error in distance (should be 0)
    lon, lat = np.zeros( N + 1 ), np.zeros( N + 1 )
    for i in range( N + 1 ):
        q = p0 + ( 1 + i * derror) * rotationMatrix( p0, i * dalpha ) @ ( p1 - p0 )      # rotate p1 about p0
        R, lon[i], lat[i] = cartesian_to_lonlat( q[0], q[1], q[2] )
    return lon, lat

########################################################################

p0_lon, p0_lat = 14.265608, 50.095867 #
p1_lon, p1_lat = 14.475131, 50.052842 #   Original points - too close to be interesting
p2_lon, p2_lat = 14.416322, 49.992767 #

p0_lon, p0_lat =  10      ,  0        #
p1_lon, p1_lat = -25      ,  -30      #   Alternatives
p2_lon, p2_lat =  40      ,  +30      #

N = 10
lon, lat = arc( p0_lon, p0_lat, p1_lon, p1_lat, p2_lon, p2_lat, N )

########################################################################

pts = np.zeros( (3,3) )
pts[:,0] = lonlat_to_cartesian( 1.0, p0_lon, p0_lat )
pts[:,1] = lonlat_to_cartesian( 1.0, p1_lon, p1_lat )
pts[:,2] = lonlat_to_cartesian( 1.0, p2_lon, p2_lat )
edge = np.zeros( (3,N+1) )
for i in range( N + 1 ): 
    edge[:,i] = lonlat_to_cartesian( 1.0, lon[i], lat[i] )

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')
plotSphere( 1 )
ax.scatter( pts[0], pts[1], pts[2], 'b' )
ax.plot( edge[0,:], edge[1,:], edge[2,:], 'r-' )
ax.set_aspect('equal')
plt.show()

Original points:

More interesting points: