#!/usr/bin/env python
# found this somewhere pasted but can't recall.
# bare-minimum requirement to do geo stuff with python when you can't benefit from
# the numerous and powerful libs that normally do this.

import math

def haversine_distance(p1, p2):
    # points are lat, lon
    # and degrees

    sin = math.sin
    cos = math.cos

    lat1, lon1 = p1
    lat2, lon2 = p2

    lat1 = lat1 / 180 * math.pi
    lat2 = lat2 / 180 * math.pi
    lon1 = lon1 / 180 * math.pi
    lon2 = lon2 / 180 * math.pi

    dlat = lat2 - lat1
    dlon = lon2 - lon1

    sin_dlat = sin(dlat / 2)
    sin_dlon = sin(dlon / 2)

    return 2.0 * math.asin(math.sqrt(sin_dlat * sin_dlat + cos(lat1) * cos(lat2) * sin_dlon * sin_dlon))

def sphere_cos_distance(p1, p2):
    # points are lat, lon
    # and degrees

    sin = math.sin
    cos = math.cos

    lat1, lon1 = p1
    lat2, lon2 = p2

    lat1 = lat1 / 180 * math.pi
    lat2 = lat2 / 180 * math.pi
    lon1 = lon1 / 180 * math.pi
    lon2 = lon2 / 180 * math.pi

    dlon = lon2 - lon1

    return math.acos(sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(dlon))

def linear_distance(p1, p2):
    sin = math.sin
    cos = math.cos

    lat1, lon1 = p1
    lat2, lon2 = p2

    lat1 = lat1 / 180 * math.pi
    lat2 = lat2 / 180 * math.pi
    lon1 = lon1 / 180 * math.pi
    lon2 = lon2 / 180 * math.pi

    average_lat = (lat1 + lat2) / 2

    dlat = lat2 - lat1
    dlon = lon2 - lon1

    scaled_dlon = dlon * cos(average_lat)

    return math.sqrt(dlat * dlat + scaled_dlon * scaled_dlon)

def test_points(p1, p2):
    ha_dist = haversine_distance(p1, p2)
#    cos_dist = sphere_cos_distance(p1, p2)
    cos_dist = linear_distance(p1, p2)

    print "%s %s %s %s" % (ha_dist, cos_dist, ha_dist / (ha_dist - cos_dist), ha_dist * 6375000)

def main():
    p1 = (37.5, 122.5)
    p2 = (37.5, 122.5 + 1e-5)
    p3 = (37.5, 122.5 + 1e-4)
    p4 = (37.5 + 1e-5, 122.5 - 1e-5)
    test_points(p1, p2)
    test_points(p1, p3)
    test_points(p1, p4)

if __name__ == "__main__":
    main()