Fixed sphere calculation

ogn/ognutils.py

@@ -1,6 +1,7 @@
 import requests
 import csv
 from io import StringIO
+from math import sin, cos, asin, atan2, sqrt, pi
 
 from .model import Flarm, AddressOrigin
 
@@ -8,6 +9,9 @@ from geopy.geocoders import Nominatim
 
 DDB_URL = "http://ddb.glidernet.org/download"
 
+deg2rad = pi/180
+rad2deg = 1/deg2rad
+
 
 def get_ddb(csvfile=None):
     if csvfile is None:
@@ -47,3 +51,14 @@ def get_country_code(latitude, longitude):
     except KeyError:
         country_code = None
     return country_code
+
+
+def wgs84_to_sphere(lat1, lat2, lon1, lon2, alt1, alt2):
+    lat1 = lat1*deg2rad
+    lat2 = lat2*deg2rad
+    lon1 = lon1*deg2rad
+    lon2 = lon2*deg2rad
+    radius = 6366000*2*asin(sqrt((sin((lat1-lat2)/2))**2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))**2))
+    theta = atan2(alt2-alt1, radius)*rad2deg
+    phi = atan2(sin(lon2-lon1)*cos(lat2), cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon2-lon1))*rad2deg
+    return radius, theta, phi