こんにちは。
回転楕円体面上の二点間距離の計算は測地線の逆問題ですが¹、短距離条件で近似計算²を行ってみました。下記の近似精度評価例では、12.457km の距離に対する計算誤差は 1.3mm です。

$ ./geodesic_inverse_approx.py 60.0 60.1 0.1
input: lat1= 60.000000  lat2= 60.100000  lam12=     0.100000
output: az1= 26.526      az2= 26.612       s12= 12456.777
error:  az1= 1.2e-05     az2= 1.2e-05      s12=   1.3e-03

geodesic_inverse_approx.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division
from __future__ import print_function
import sys
import numpy as np
from easydict import EasyDict as edict
from geographiclib.geodesic import Geodesic

RE = 6378137.0 # IS-GPS
FE = (1/298.257223563) # IS-GPS
E2 = FE * (2 - FE)
DEGREE = np.pi/180
geod = Geodesic.WGS84

# "Algorithms for geodesics" Charles F. F. Karney (2013)
#  5. Starting point for Newton’s method
#  solving for the great circle on the auxiliary sphere
def geodesic_inverse_approx(beta1, beta2, lam12):
    cosBeta = (beta1.cos + beta2.cos)/2
    wm = np.sqrt(1-E2*cosBeta**2)
    omg12 = angle_approx(lam12/wm)
    az1, r = angle(beta1.cos*beta2.sin - beta1.sin*beta2.cos*omg12.cos, beta2.cos*omg12.sin)
    az2, _ = angle(-beta1.sin*beta2.cos + beta1.cos*beta2.sin*omg12.cos, beta1.cos*omg12.sin)
    sig12 = arg_approx(beta1.sin*beta2.sin + beta1.cos*beta2.cos*omg12.cos, r)
    s12 = sig12 * wm
    return az1, az2, s12

def reducedLatitude(lat):
    chi = np.tan(lat/2)
    beta, _ = angle(1-chi**2, 2*chi*(1-FE))
    return beta

def arg(*args):
    if len(args) == 1:
        c, s = args[0].cos, args[0].sin
    else:
        c, s = args[0], args[1]
    return np.arctan2(s, c)

def arg_approx(c, s):
    x = s/c
    return x*(1-x**2/3)

def angle(*args):
    if len(args) == 1:
        x = np.tan(args[0]/2)
        return angle(1-x**2, 2*x)
    r = np.hypot(args[0], args[1])
    c, s = args[0]/r, args[1]/r
    return edict({'cos': c, 'sin': s}), r

def angle_approx(x):
    return edict({'cos': 1-x**2/2, 'sin': x})

def print_geodesic_inverse(lat1, lat2, lam12):
    print("input: lat1=%10.6f  lat2=%10.6f  lam12=%13.6f" % (lat1, lat2, lam12))
    beta1 = reducedLatitude(lat1*DEGREE)
    beta2 = reducedLatitude(lat2*DEGREE)
    az1, az2, s12 = geodesic_inverse_approx(beta1, beta2, lam12*DEGREE)
    az1, az2, s12 = arg(az1)/DEGREE, arg(az2)/DEGREE, s12*RE
    print("output: az1=%7.3f      az2=%7.3f       s12=%10.3f" % (az1, az2, s12))
    g = edict(geod.Inverse(lat1, 0, lat2, lam12))
    print("error:  az1=%8.1e     az2=%8.1e      s12=%10.1e" % (az1-g.azi1, az2-g.azi2, s12-g.s12))
    return

def main():
    args = sys.argv[1:]
    if len(args) == 3:
        print_geodesic_inverse(*map(lambda x:np.float(x), args))
    return

if __name__ == '__main__':
    main()
    exit(0)

これは GeographicLib を使って計算できます（Inverse()）。 ↩
これは、GeographicLib の逆問題計算の中の Newton 法向け初期値計算にも使われている方法です。 ↩