#!/usr/bin/env python2
#
# Copyright © 2010 Elizabeth Jennifer Myers. 
# All rights reserved.
# 
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#    * Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#    * Redistributions in binary form must reproduce the above copyright
#      notice, this list of conditions and the following disclaimer in the
#      documentation and/or other materials provided with the distribution.
#    * Neither the name of the copyright holders nor the names of its
#      contributors may be used to endorse or promote products
#      derived from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER BE LIABLE FOR ANY
# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

# =============================================================================

# Based on the formula for great-circle distance
# (http://en.wikipedia.org/wiki/Great-circle_distance)
#
# Note that this is not 100% accurate as the earth isn't a "true" sphere, but
# for our purposes it's plenty (guaranteed within 0.5%).

import math
                                                                                                                                                             
def GeoDistance(lat1, lon1, lat2, lon2):                                                                                                                     
    """                                                                                                                                                      
    Calculate the distance between two latitude and longitude pairs on earth.                                                                                
                                                                                                                                                             
    Note that negative numbers denote west and south, and positive numbers                                                                                   
    denote east and north.                                                                                                                                   
    """                                                                                                                                                      
                                                                                                                                                             
    phi1 = math.radians(90.0 - lat1)                                                                                                                         
    phi2 = math.radians(90.0 - lat2)                                                                                                                         
    theta1 = math.radians(lon1)                                                                                                                              
    theta2 = math.radians(lon2)                                                                                                                              
                                                                                                                                                             
    arc = math.acos(math.sin(phi1) * math.sin(phi2) * 
            math.cos(theta1 - theta2) + math.cos(phi1) * math.cos(phi2))

    # 6335.439 is the average radius of the earth in km.
    return arc * 6335.439