# coding:utf-8

import math

R = 6371000

def calc_greate_circle_distance(lat1, lon1, lat2, lon2):

    lat1_rad = lat1 * math.pi / 180

    lon1_rad = lon1 * math.pi / 180

    lat2_rad = lat2 * math.pi / 180

    lon2_rad = lon2 * math.pi / 180

    d = math.sqrt(math.pow(math.cos((lat1_rad + lat2_rad) / 2) * R * (lon1_rad - lon2_rad), 2) +

                  math.pow(R * (lat1_rad - lat2_rad), 2))

    if d < 10000:

        return float(d)

    else:

        delta_lat = lat2_rad - lat1_rad

        delta_lon = lon2_rad - lon1_rad

        a = math.sin(delta_lat / 2) * math.sin(delta_lat / 2) + \

            math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(delta_lon / 2) * math.sin(delta_lon / 2)

        c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))

        return float(R * c)

def calc_azimuth(lat1, lon1, lat2, lon2):

    lat1_rad = lat1 * math.pi / 180

    lon1_rad = lon1 * math.pi / 180

    lat2_rad = lat2 * math.pi / 180

    lon2_rad = lon2 * math.pi / 180

    y = math.sin(lon2_rad - lon1_rad) * math.cos(lat2_rad)

    x = math.cos(lat1_rad) * math.sin(lat2_rad) - \

        math.sin(lat1_rad) * math.cos(lat2_rad) * math.cos(lon2_rad - lon1_rad)

    brng = math.atan2(y, x) * 180 / math.pi

    return float((brng + 360.0) % 360.0)

所谓方位角是与正北方向、顺时针之间的夹角