calculate the shortest distance along the surface of an ellipsoid between two points on the surface of earth given longitudes and latitudes https en wikipedia orgwikigeographicaldistancelambert sformulaforlonglines note this algorithm uses geodesyhaversinedistance py to compute central angle sigma representing the earth as an ellipsoid allows us to approximate distances between points on the surface much better than a sphere ellipsoidal formulas treat the earth as an oblate ellipsoid which means accounting for the flattening that happens at the north and south poles lambert s formulae provide accuracy on the order of 10 meteres over thousands of kilometeres other methods can provide millimeterlevel accuracy but this is a simpler method to calculate long range distances without increasing computational intensity args lat1 lon1 latitude and longitude of coordinate 1 lat2 lon2 latitude and longitude of coordinate 2 returns geographical distance between two points in metres from collections import namedtuple point2d namedtuplepoint2d lat lon sanfrancisco point2d37 774856 122 424227 yosemite point2d37 864742 119 537521 newyork point2d40 713019 74 012647 venice point2d45 443012 12 313071 flambertsellipsoidaldistancesanfrancisco yosemite 0 0f meters 254 351 meters flambertsellipsoidaldistancesanfrancisco newyork 0 0f meters 4 138 992 meters flambertsellipsoidaldistancesanfrancisco venice 0 0f meters 9 737 326 meters constants per wgs84 https en wikipedia orgwikiworldgeodeticsystem distance in metresm equation parameters https en wikipedia orgwikigeographicaldistancelambert sformulaforlonglines parametric latitudes https en wikipedia orgwikilatitudeparametricorreducedlatitude compute central angle between two points using haversine theta sigma haversinedistance equatorial radius intermediate p and q values intermediate x value x sigma sinsigma sin2pcos2q cos2sigma2 intermediate y value y sigma sinsigma cos2psin2q sin2sigma2 calculate the shortest distance along the surface of an ellipsoid between two points on the surface of earth given longitudes and latitudes https en wikipedia org wiki geographical_distance lambert s_formula_for_long_lines note this algorithm uses geodesy haversine_distance py to compute central angle sigma representing the earth as an ellipsoid allows us to approximate distances between points on the surface much better than a sphere ellipsoidal formulas treat the earth as an oblate ellipsoid which means accounting for the flattening that happens at the north and south poles lambert s formulae provide accuracy on the order of 10 meteres over thousands of kilometeres other methods can provide millimeter level accuracy but this is a simpler method to calculate long range distances without increasing computational intensity args lat1 lon1 latitude and longitude of coordinate 1 lat2 lon2 latitude and longitude of coordinate 2 returns geographical distance between two points in metres from collections import namedtuple point_2d namedtuple point_2d lat lon san_francisco point_2d 37 774856 122 424227 yosemite point_2d 37 864742 119 537521 new_york point_2d 40 713019 74 012647 venice point_2d 45 443012 12 313071 f lamberts_ellipsoidal_distance san_francisco yosemite 0 0f meters 254 351 meters f lamberts_ellipsoidal_distance san_francisco new_york 0 0f meters 4 138 992 meters f lamberts_ellipsoidal_distance san_francisco venice 0 0f meters 9 737 326 meters constants per wgs84 https en wikipedia org wiki world_geodetic_system distance in metres m equation parameters https en wikipedia org wiki geographical_distance lambert s_formula_for_long_lines parametric latitudes https en wikipedia org wiki latitude parametric_ or_reduced _latitude compute central angle between two points using haversine theta sigma haversine_distance equatorial radius intermediate p and q values intermediate x value x sigma sin sigma sin 2pcos 2q cos 2 sigma 2 intermediate y value y sigma sin sigma cos 2psin 2q sin 2 sigma 2
from math import atan, cos, radians, sin, tan

from .haversine_distance import haversine_distance

AXIS_A = 6378137.0
AXIS_B = 6356752.314245
EQUATORIAL_RADIUS = 6378137


def lamberts_ellipsoidal_distance(
    lat1: float, lon1: float, lat2: float, lon2: float
) -> float:
    

    
    
    
    
    flattening = (AXIS_A - AXIS_B) / AXIS_A
    
    
    b_lat1 = atan((1 - flattening) * tan(radians(lat1)))
    b_lat2 = atan((1 - flattening) * tan(radians(lat2)))

    
    
    sigma = haversine_distance(lat1, lon1, lat2, lon2) / EQUATORIAL_RADIUS

    
    p_value = (b_lat1 + b_lat2) / 2
    q_value = (b_lat2 - b_lat1) / 2

    
    
    x_numerator = (sin(p_value) ** 2) * (cos(q_value) ** 2)
    x_demonimator = cos(sigma / 2) ** 2
    x_value = (sigma - sin(sigma)) * (x_numerator / x_demonimator)

    
    
    y_numerator = (cos(p_value) ** 2) * (sin(q_value) ** 2)
    y_denominator = sin(sigma / 2) ** 2
    y_value = (sigma + sin(sigma)) * (y_numerator / y_denominator)

    return EQUATORIAL_RADIUS * (sigma - ((flattening / 2) * (x_value + y_value)))


if __name__ == "__main__":
    import doctest

    doctest.testmod()