hubeny

ilh.py

"""
Copyright(c) 2015 ryo hashiguchi
This software is released under the MIT License.
http://opensource.org/licenses/mit-license.php
"""
import math
import sys


def convert(data):
    data = list(map(lambda x: float(x), data.split(" ")))
    return (data[0] + data[1] / 60.0 + data[2] / 3600.0) * math.pi / 180


def hubeny(pos1, pos2, display=False):
    lat_ave = (pos1[0] + pos2[0]) / 2
    d_lat = abs(pos1[0] - pos2[0])
    d_lon = abs(pos1[1] - pos2[1])
    M = 6334832.10663254 / math.sqrt((1 - 0.006674 * math.sin(lat_ave) * math.sin(lat_ave)) ** 3)
    N = 6377397.155 / math.sqrt(1 - 0.006674 * math.sin(lat_ave) * math.sin(lat_ave))
    if display:
        print("P: (%.10f + %.10f) / 2 = %.10f" % (pos1[0], pos2[0], lat_ave))
        print("dp: |%.10f - %.10f| = %.10f" % (pos1[0], pos2[0], d_lat))
        print("dr: |%.10f - %.10f| = %.10f" % (pos1[1], pos2[1], d_lon))
        print("M: " + str(M))
        print("N: " + str(N))
    return math.sqrt((M * d_lat) * (M * d_lat) + (N * math.cos(lat_ave) * d_lon) * (N * math.cos(lat_ave) * d_lon))


if __name__ == "__main__":
    while True:
        text = input('')
        if len(text) == 0:
            break
        lat1, lon1, lat2, lon2 = text.split(" ")
        print()
        detail = len(sys.argv) > 1 and sys.argv[1] == "true"
        d = hubeny([convert("31 43 " + lat1), convert("130 43 " + lon1)], [convert("31 43 " + lat2), convert("130 43 " + lon2)], detail)
        print("D: " + str(d if detail else round(d, 3)))