Oct-04-2016, 10:02 PM
This is a small part of a bigger project I'm working on to Pythonize (a real word I just made up) the formulas in Jean Meeus' book "Astronomical Algorithms" copywrite 1991, pages 80 and 81 and calculates the geodesic distance between two points on the earths surface. The two co-ordinates referenced at the end are, per Mr. Meeus for the Observatoire de Paris (France) and the US Naval Observatory (Washington DC) with an accuracy of 50 meters. Comments and improvements are welcome.
#! /usr/bin/python3 # -*- coding: utf-8 -*- from math import atan, cos, pi, sin, sqrt def calc_geo_dist(lon1, lat1, lon2, lat2): a = 6378.137 # in km, semi-major axis # b = 6356.752 # polar radius in km, rounded approx 6356752.314140347 meters per GRS 80 rad = pi/180 flat = 1/298.257 # ok with nasa actual 0.003352810681183637418 e = ((lon1 - lon2) / 2) * rad f = ((lat1+lat2) / 2) * rad g = ((lat1-lat2) / 2) * rad s = sin(g) * sin(g) * cos(e) * cos(e) + cos(f) * cos(f) * sin(e) * sin(e) c = cos(g) * cos(g) * cos(e) * cos(e) + sin(f) * sin(f) * sin(e) * sin(e) w = atan(sqrt(s/c)) r = sqrt(s * c) / w d = 2 * w * a h1 = (3 * r - 1) / (2 * c) h2 = (3 * r + 1) / (2 * s) s1 = d * (1 + flat * h1 * sin(f) * sin(f) * cos(g) * cos(g) - flat * h2 * cos(f) * cos(f) * sin(g) * sin(g)) return s1 while True: try: my_long1 = float(input('\nEnter longitude in decimal degrees: ')) if -180 <= my_long1 <= 180: break else: print('\n WARNING: Value must be between -180 and 180 ') except ValueError as err: print('\n ERROR: Enter value between -180 and 180 {}'.format(err)) while True: try: my_lat1 = float(input('Enter latitude in decimal degrees: ')) if -90 <= my_lat1 <= 90: break else: print('\n WARNING: Value must be between -90 and 90 ') except ValueError: print('\n ERROR: Enter value between -90 and 90 ') distance = input('\nDo you require the distance between two locations (y or n)? ') if distance == 'y': my_long2 = float(input('For distance, enter second longitude in decimal degrees: ')) my_lat2 = float(input('Enter second latitude in decimal degrees: ')) if distance == 'y': dbp = calc_geo_dist(my_long1, my_lat1, my_long2, my_lat2) print('The distance between the two points = {} km'.format(str(dbp))) # lon1 = -2.3372222, lat1 = 48.8363889; lon2 = 77.0655556, lat 2 = 38.9213889. Should be 6181.63 km
If it ain't broke, I just haven't gotten to it yet.
OS: Windows 10, openSuse 42.3, freeBSD 11, Raspian "Stretch"
Python 3.6.5, IDE: PyCharm 2018 Community Edition
OS: Windows 10, openSuse 42.3, freeBSD 11, Raspian "Stretch"
Python 3.6.5, IDE: PyCharm 2018 Community Edition