diff --git a/sbapp/sideband/geo.py b/sbapp/sideband/geo.py index 0932607..c4ca0d5 100644 --- a/sbapp/sideband/geo.py +++ b/sbapp/sideband/geo.py @@ -223,42 +223,78 @@ def orthodromic_distance(c1, c2, ellipsoid=True): else: return spherical_distance(c1, c2) -# def tests(): -# import RNS -# import numpy as np -# from geographiclib.geodesic import Geodesic -# geod = Geodesic.WGS84 -# coords = [ -# [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 0.0)], -# [(51.230800, 4.38703, 0.0), (51.230801, 4.38703, 0.0)], -# [(35.3524, 135.0302, 100), (35.3532,135.0305, 500)], -# [(57.758793, 22.605194, 0.0), (43.048838, -9.241343, 0.0)], -# [(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)], -# [(-90.0, 0.0, 0.0), (90.0, 0.0, 0.0)], -# [(-90.0, 0.0, 0.0), (78.0, 0.0, 0.0)], -# [(0.0, 0.0, 0.0), (0.5, 179.5, 0.0)], -# [(0.7, 0.0, 0.0), (0.0, -180.0, 0.0)], -# ] -# for cs in coords: -# c1 = cs[0]; c2 = cs[1] -# print("Testing: "+str(c1)+" -> "+str(c2)) -# us = time.time() -# ld = c1+c2; g = geod.Inverse(c1[0], c1[1], c2[0], c2[1]) -# print("Lib computed in "+str(round((time.time()-us)*1e6, 3))+"us") -# us = time.time() -# eld = orthodromic_distance(c1,c2,ellipsoid=True) -# if eld: -# print("Own computed in "+str(round((time.time()-us)*1e6, 3))+"us") -# else: -# print("Own timed out in "+str(round((time.time()-us)*1e6, 3))+"us") -# ed_own = euclidian_distance(c1,c2,ellipsoid=True) -# sd_own = orthodromic_distance(c1,c2,ellipsoid=False) -# aa = azalt(c1,c2,ellipsoid=True) -# fac = 1 -# if eld: print("LibDiff = "+RNS.prettydistance(g['s12']-eld)+f" {fac*g['s12']-fac*eld}") -# print("Spherical = "+RNS.prettydistance(sd_own)+f" {fac*sd_own}") -# # print("EllipLib = "+RNS.prettydistance(g['s12'])+f" {fac*g['s12']}") -# if eld: print("Ellipsoid = "+RNS.prettydistance(eld)+f" {fac*eld}") -# print("Euclidian = "+RNS.prettydistance(ed_own)+f" {fac*ed_own}") -# print("AzAlt = "+f" {aa[0]} / {aa[1]}") -# print("") +def distance_to_horizon(c, ellipsoid=False): + if ellipsoid: + raise NotImplementedError("Distance to horizon on the ellipsoid is not yet implemented") + else: + # TODO: This is a only barely functional simplification. + # Need to calculate the geodesic distance to the horizon + # instead. + if len(c) >= 3: + r = mean_earth_radius + h = c[2] + return sqrt(pow((h+r),2) - r*r) + else: + return None + +def angle_to_horizon(c, ellipsoid=False): + if ellipsoid: + raise NotImplementedError("Angle to horizon on the ellipsoid is not yet implemented") + else: + r = mean_earth_radius + h = c[2] + return degrees(-acos(r/(r+h))) + +def radio_horizon(c1, c2, ellipsoid=False): + # dr = 4.12*(√h1 + √h2) + if ellipsoid: + raise NotImplementedError("Radio horizon on the ellipsoid is not yet implemented") + else: + h1 = c1[2] + h2 = c2[2] + ed = euclidian_distance(c1,c2) + rh1 = 1e3*4.12*(sqrt(h1)) + rh2 = 1e3*4.12*(sqrt(h2)) + rhc = 1e3*4.12*(sqrt(h1) + sqrt(h2)) + return (rh1, rh2, rhc, rhc > ed) + +def tests(): + import RNS + import numpy as np + from geographiclib.geodesic import Geodesic + geod = Geodesic.WGS84 + coords = [ + [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 0.0)], + [(51.2308, 4.38703, 0.0), (47.699437, 9.268651, 30.0*1e3)], + # [(51.230800, 4.38703, 0.0), (51.230801, 4.38703, 0.0)], + # [(35.3524, 135.0302, 100), (35.3532,135.0305, 500)], + # [(57.758793, 22.605194, 0.0), (43.048838, -9.241343, 0.0)], + # [(0.0, 0.0, 0.0), (0.0, 0.0, 0.0)], + # [(-90.0, 0.0, 0.0), (90.0, 0.0, 0.0)], + # [(-90.0, 0.0, 0.0), (78.0, 0.0, 0.0)], + # [(0.0, 0.0, 0.0), (0.5, 179.5, 0.0)], + # [(0.7, 0.0, 0.0), (0.0, -180.0, 0.0)], + ] + for cs in coords: + c1 = cs[0]; c2 = cs[1] + print("Testing: "+str(c1)+" -> "+str(c2)) + us = time.time() + ld = c1+c2; g = geod.Inverse(c1[0], c1[1], c2[0], c2[1]) + print("Lib computed in "+str(round((time.time()-us)*1e6, 3))+"us") + us = time.time() + eld = orthodromic_distance(c1,c2,ellipsoid=True) + if eld: + print("Own computed in "+str(round((time.time()-us)*1e6, 3))+"us") + else: + print("Own timed out in "+str(round((time.time()-us)*1e6, 3))+"us") + ed_own = euclidian_distance(c1,c2,ellipsoid=True) + sd_own = orthodromic_distance(c1,c2,ellipsoid=False) + aa = azalt(c1,c2,ellipsoid=True) + fac = 1 + if eld: print("LibDiff = "+RNS.prettydistance(g['s12']-eld)+f" {fac*g['s12']-fac*eld}") + print("Spherical = "+RNS.prettydistance(sd_own)+f" {fac*sd_own}") + # print("EllipLib = "+RNS.prettydistance(g['s12'])+f" {fac*g['s12']}") + if eld: print("Ellipsoid = "+RNS.prettydistance(eld)+f" {fac*eld}") + print("Euclidian = "+RNS.prettydistance(ed_own)+f" {fac*ed_own}") + print("AzAlt = "+f" {aa[0]} / {aa[1]}") + print("")