Added angle-to-horizon and radio horizon calculations
This commit is contained in:
Mark Qvist
parent
ecb5f0c38b
commit
9fe7632e8f
|
|
@ -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("")
|
||||
|
|
|
|||
Loading…
Reference in New Issue