Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use njit on more functions, to increase speed #125

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions orbit_predictor/predictors/pass_iterators.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
48 changes: 32 additions & 16 deletions orbit_predictor/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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):
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)))
Expand Down Expand Up @@ -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).
Expand All @@ -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)

Expand Down Expand Up @@ -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:
Expand Down
1 change: 1 addition & 0 deletions requirements.txt
Original file line number Diff line number Diff line change
Expand Up @@ -9,3 +9,4 @@ requests>=2.9.1
sgp4>=1.4
numpy
scipy
numba>=0.54.0
2 changes: 1 addition & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@
],
extras_require={
"fast": [
"numba>=0.38",
"numba>=0.54.0",
"scipy>=0.16",
],
"dev": [
Expand Down