From be491063ca25a88daea0ec2c3be57bce49a54a0c Mon Sep 17 00:00:00 2001 From: Pablo Jais Date: Wed, 31 Aug 2022 16:09:19 +0200 Subject: [PATCH 1/5] Use njit on more functions, to increase speed --- orbit_predictor/utils.py | 48 ++++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index 9f7dad2..c299d85 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -57,27 +57,36 @@ def _jit(f): DECEMBER_31TH_1999_MIDNIGHT_JD = 2451543.5 -def compose(*functions): - """Performs function composition with variadic arguments""" - return functools.reduce(lambda f, g: lambda *args: f(g(*args)), - functions, - lambda x: x) +@njit +def cos_d(angle_d): + return cos(radians(angle_d)) + + +@njit +def sin_d(angle_d): + return sin(radians(angle_d)) -cos_d = compose(cos, radians) -sin_d = compose(sin, radians) -atan2_d = compose(degrees, atan2) -asin_d = compose(degrees, asin) +@njit +def atan2_d(value1, value2): + return degrees(atan2(value1, value2)) + + +@njit +def asin_d(value): + return degrees(asin(value)) AzimuthElevation = namedtuple('AzimuthElevation', 'azimuth elevation') +@njit def euclidean_distance(*components): """Returns the norm of a vector""" - return sqrt(sum(c**2 for c in components)) + return sqrt(sum([c**2 for c in components])) +@njit def angle_between(a, b): """ Computes angle between two vectors in degrees. @@ -87,12 +96,13 @@ def angle_between(a, b): Naïve algorithm, see https://scicomp.stackexchange.com/q/27689/782. """ - return degrees(np.arccos(dot_product(a, b) / (vector_norm(a) * vector_norm(b)))) + return degrees(acos(dot_product(a, b) / (vector_norm(a) * vector_norm(b)))) +@njit def dot_product(a, b): """Computes dot product between two vectors writen as tuples or lists""" - return sum(ai * bj for ai, bj in zip(a, b)) + return sum([ai * bj for ai, bj in zip(a, b)]) def vector_diff(a, b): @@ -109,9 +119,10 @@ def cross_product(a, b): ) -def vector_norm(a): +@njit +def vector_norm(components): """Returns the norm of a vector""" - return euclidean_distance(*a) + return sqrt(sum([c**2 for c in components])) @njit @@ -266,6 +277,7 @@ def sun_azimuth_elevation(latitude_deg, longitude_deg, when=None): return AzimuthElevation(azimuth, elevation) +@njit def _sun_mean_ecliptic_elements(t_ut1): w = 282.9404 + 4.70935e-5 * t_ut1 # longitude of perihelion degrees eccentricity = 0.016709 - 1.151e-9 * t_ut1 # eccentricity @@ -276,6 +288,7 @@ def _sun_mean_ecliptic_elements(t_ut1): return w, M, L, eccentricity, oblecl +@njit def _sun_eci(w, M, L, eccentricity, oblecl): # auxiliary angle auxiliary_angle = M + degrees(eccentricity * sin_d(M) * (1 + eccentricity * cos_d(M))) @@ -335,6 +348,7 @@ def get_shadow(r, when_utc): return shadow(r_sun, ecef_to_eci(r, gmst)) +@njit def shadow(r_sun, r, r_p=R_E_MEAN_KM): """ Gives illumination of Earth satellite (2 for illuminated, 1 for penumbra, 0 for umbra). @@ -358,8 +372,9 @@ def shadow(r_sun, r, r_p=R_E_MEAN_KM): if dot_product(r_sun, r) < 0: angle = angle_between(-r_sun, r) - sat_horiz = vector_norm(r) * cos_d(angle) - sat_vert = vector_norm(r) * sin_d(angle) + norm_r = vector_norm(r) + sat_horiz = norm_r * cos_d(angle) + sat_vert = norm_r * sin_d(angle) x = r_p / sin(ALPHA_PEN) pen_vert = tan(ALPHA_PEN) * (x + sat_horiz) @@ -432,6 +447,7 @@ def eclipse_duration(beta, period, r_p=R_E_MEAN_KM): ) * period / pi +@njit def juliandate(utc_tuple): year, month, day, hour, minute, sec = utc_tuple[:6] if month <= 2: From f0a4be614b894fed84a2d26a02cc7008b037e673 Mon Sep 17 00:00:00 2001 From: Pablo Jais Date: Tue, 6 Sep 2022 18:37:56 +0200 Subject: [PATCH 2/5] Replace sum by np.sum to be compatible with older versions of numba --- orbit_predictor/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index c299d85..8af37c3 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -83,7 +83,7 @@ def asin_d(value): @njit def euclidean_distance(*components): """Returns the norm of a vector""" - return sqrt(sum([c**2 for c in components])) + return sqrt(np.sum([c**2 for c in components])) @njit @@ -102,7 +102,7 @@ def angle_between(a, b): @njit def dot_product(a, b): """Computes dot product between two vectors writen as tuples or lists""" - return sum([ai * bj for ai, bj in zip(a, b)]) + return np.sum([ai * bj for ai, bj in zip(a, b)]) def vector_diff(a, b): @@ -122,7 +122,7 @@ def cross_product(a, b): @njit def vector_norm(components): """Returns the norm of a vector""" - return sqrt(sum([c**2 for c in components])) + return sqrt(np.sum([c**2 for c in components])) @njit From 36180b6f76a7cd67e806dc489abe4bf48bebe0ad Mon Sep 17 00:00:00 2001 From: Pablo Jais Date: Tue, 6 Sep 2022 18:49:52 +0200 Subject: [PATCH 3/5] Add handling of perfect sate_pos and target alignment in off_nadir_deg --- orbit_predictor/predictors/pass_iterators.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/orbit_predictor/predictors/pass_iterators.py b/orbit_predictor/predictors/pass_iterators.py index 11e73e2..9a9e122 100644 --- a/orbit_predictor/predictors/pass_iterators.py +++ b/orbit_predictor/predictors/pass_iterators.py @@ -349,6 +349,8 @@ def off_nadir_deg(self): sign = dot / abs(dot) except ZeroDivisionError: sign = 1 + if np.isnan(sign): + sign = 1 return degrees(angle) * sign From ec645f16210e018f7ff2a80e9f6d757f39bfacd4 Mon Sep 17 00:00:00 2001 From: Pablo Jais Date: Tue, 6 Sep 2022 19:22:26 +0200 Subject: [PATCH 4/5] Revert to the version with builtin sum, but require a newer numba --- orbit_predictor/utils.py | 6 +++--- requirements.txt | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index 8af37c3..c299d85 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -83,7 +83,7 @@ def asin_d(value): @njit def euclidean_distance(*components): """Returns the norm of a vector""" - return sqrt(np.sum([c**2 for c in components])) + return sqrt(sum([c**2 for c in components])) @njit @@ -102,7 +102,7 @@ def angle_between(a, b): @njit def dot_product(a, b): """Computes dot product between two vectors writen as tuples or lists""" - return np.sum([ai * bj for ai, bj in zip(a, b)]) + return sum([ai * bj for ai, bj in zip(a, b)]) def vector_diff(a, b): @@ -122,7 +122,7 @@ def cross_product(a, b): @njit def vector_norm(components): """Returns the norm of a vector""" - return sqrt(np.sum([c**2 for c in components])) + return sqrt(sum([c**2 for c in components])) @njit diff --git a/requirements.txt b/requirements.txt index b47cc36..76ad9ad 100644 --- a/requirements.txt +++ b/requirements.txt @@ -9,3 +9,4 @@ requests>=2.9.1 sgp4>=1.4 numpy scipy +numba>=0.54.0 From 2d98e83edf9eb9bd92112a438c021b365adaffa9 Mon Sep 17 00:00:00 2001 From: Pablo Jais Date: Tue, 6 Sep 2022 19:26:32 +0200 Subject: [PATCH 5/5] Add numba requirement to setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 44b78ce..12abccb 100755 --- a/setup.py +++ b/setup.py @@ -40,7 +40,7 @@ ], extras_require={ "fast": [ - "numba>=0.38", + "numba>=0.54.0", "scipy>=0.16", ], "dev": [