Lambert-Andoyer

get_dist.py

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

import sys
import math

DEG2RAD = math.pi / 180.0

# WGS84
dLA_A = 6378137.0 # 赤道半径
dLA_F = 1.0 / 298.257222101 # 扁平率
# 扁平率 F = (A - B) / A
dLA_B = dLA_A * (1.0 - dLA_F) # 極半径
dLA_Coef1 = math.atan(dLA_B / dLA_A)

def distance_lambert(dLat1, dLon1, dLat2, dLon2):
    dLat1 = dLat1 * DEG2RAD
    dLon1 = dLon1 * DEG2RAD
    dLat2 = dLat2 * DEG2RAD
    dLon2 = dLon2 * DEG2RAD
    dLA_P1 = dLA_Coef1 * math.tan(dLat1)
    dLA_P2 = dLA_Coef1 * math.tan(dLat2)
    # Spherical Distance
    dSD = math.acos(math.sin(dLA_P1) * math.sin(dLA_P2) \
                        + math.cos(dLA_P1) * math.cos(dLA_P2) * math.cos(dLon1 - dLon2))
    # Lambert-Andoyer Correction
    dCos_SD = math.cos(dSD / 2.0)
    dSin_SD = math.sin(dSD / 2.0)
    dSin_SD2 = math.sin(dSD)
    dLA_P1 = math.sin(dLA_P1)
    dLA_P2 = math.sin(dLA_P2)
    dTemp = dLA_P1 * dLA_P1 + dLA_P2 * dLA_P2
    dTemp2 = 2.0 * dLA_P1 * dLA_P2
    dCC = (dSin_SD2 - dSD) * (dTemp + dTemp2) / (dCos_SD * dCos_SD)
    dSS = (dSin_SD2 + dSD) * (dTemp - dTemp2) / (dSin_SD * dSin_SD)
    dDelta = dLA_F / 8.0 * (dCC - dSS)
    # Geodetic Distance
    return dLA_A * (dSD + dDelta) / 1000.0

while True:
    line = raw_input()
    line = line.strip()
    tokens = line.split(' ')
    if len(tokens) >= 4:
        lat1 = float(tokens[0])
        lon1 = float(tokens[1])
        lat2 = float(tokens[2])
        lon2 = float(tokens[3])
        print distance_lambert(lat1, lon1, lat2, lon2)

http://www.gadgety.net/shin/web/php/distance.html
http://petit-noise.net/blog/%E7%B7%AF%E5%BA%A6%E7%B5%8C%E5%BA%A6%E3%81%8B%E3%82%892%E7%82%B9%E9%96%93%E3%81%AE%E8%B7%9D%E9%9B%A2%E3%82%92%E6%B1%82%E3%82%81%E3%82%8B/