From 9ef987a242661755f17d283cf7eec0bee0110610 Mon Sep 17 00:00:00 2001 From: gilgamezh Date: Fri, 10 Nov 2017 18:54:01 -0300 Subject: [PATCH] first public release --- .gitignore | 47 ++++ LICENSE | 21 ++ README.rst | 82 +++++++ examples/iss.tle | 3 + orbit_predictor/__init__.py | 0 orbit_predictor/accuratepredictor.py | 329 ++++++++++++++++++++++++++ orbit_predictor/coordinate_systems.py | 153 ++++++++++++ orbit_predictor/exceptions.py | 33 +++ orbit_predictor/locations.py | 261 ++++++++++++++++++++ orbit_predictor/predictors.py | 246 +++++++++++++++++++ orbit_predictor/sources.py | 165 +++++++++++++ orbit_predictor/utils.py | 208 ++++++++++++++++ orbit_predictor/version.py | 2 + requirements.txt | 9 + setup.cfg | 11 + setup.py | 52 ++++ test | 10 + tests/test_accurate_predictor.py | 281 ++++++++++++++++++++++ tests/test_azimuth_elevation.py | 104 ++++++++ tests/test_coordinate_systems.py | 109 +++++++++ tests/test_locations.py | 109 +++++++++ tests/test_position.py | 39 +++ tests/test_predicted_pass.py | 45 ++++ tests/test_sources.py | 184 ++++++++++++++ tests/test_sun_elevation.py | 215 +++++++++++++++++ tests/test_tle_predictor.py | 279 ++++++++++++++++++++++ 26 files changed, 2997 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 README.rst create mode 100644 examples/iss.tle create mode 100644 orbit_predictor/__init__.py create mode 100644 orbit_predictor/accuratepredictor.py create mode 100644 orbit_predictor/coordinate_systems.py create mode 100644 orbit_predictor/exceptions.py create mode 100644 orbit_predictor/locations.py create mode 100644 orbit_predictor/predictors.py create mode 100644 orbit_predictor/sources.py create mode 100644 orbit_predictor/utils.py create mode 100644 orbit_predictor/version.py create mode 100644 requirements.txt create mode 100644 setup.cfg create mode 100755 setup.py create mode 100755 test create mode 100644 tests/test_accurate_predictor.py create mode 100644 tests/test_azimuth_elevation.py create mode 100644 tests/test_coordinate_systems.py create mode 100644 tests/test_locations.py create mode 100644 tests/test_position.py create mode 100644 tests/test_predicted_pass.py create mode 100644 tests/test_sources.py create mode 100644 tests/test_sun_elevation.py create mode 100644 tests/test_tle_predictor.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0d67be0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,47 @@ +*.py[cod] +local_settings.py + +# C extensions +*.so + +# Packages +*.egg +*.egg-info +dist +build +eggs +parts +bin +var +sdist +develop-eggs +.installed.cfg +lib +lib64 +__pycache__ + +# Installer logs +pip-log.txt + +# Unit test / coverage reports +.coverage +.tox +nosetests.xml + +# Translations +*.mo + +# Mr Developer +.mr.developer.cfg +.project +.pydevproject +media +celerybeat-schedule +static + +#rope +.ropeproject + + +# hypothesis +.hypothesis diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..5d81999 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2017 Satellogic SA + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.rst b/README.rst new file mode 100644 index 0000000..1a44233 --- /dev/null +++ b/README.rst @@ -0,0 +1,82 @@ +Orbit Predictor +=============== + +Orbit Predictor is a Python library to propagate orbits of Earth-orbiting objects (satellites, ISS, +Santa Claus, etc) using `TLE (Two-Line Elements set) `_ + +Al the hard work is done by The Brandon Rhodes implementation of +`SGP4 `_. + +We can say *Orbit predictor* is kind of a "wrapper" for the python implementation of SGP4 + +To install it +------------- + +You can install orbit-predictor from pypi:: + + pip install orbit-predictor # WIP + +Use example +----------- + +When will be the ISS over Argentina? + +:: + + In [1]: from orbit_predictor.sources import EtcTLESource + + In [2]: from orbit_predictor.locations import ARG + + In [3]: source = EtcTLESource(filename="examples/iss.tle") + + In [4]: predictor = source.get_predictor("ISS") + + In [5]: predictor.get_next_pass(ARG) + Out[5]: + + In [6]: predicted_pass = _ + + In [7]: position = predictor.get_position(predicted_pass.aos) + + In [8]: ARG.is_visible(position) # Can I see the ISS from this location? + Out[8]: True + + In [9]: import datetime + + In [10]: position_delta = predictor.get_position(predicted_pass.los + datetime.timedelta(minutes=20)) + + In [11]: ARG.is_visible(position_delta) + Out[11]: False + + In [12]: tomorrow = datetime.datetime.utcnow() + datetime.timedelta(days=1) + + In [13]: predictor.get_next_pass(ARG, tomorrow, max_elevation_gt=20) + Out[13]: + + +`WSTLESource` needs the tle.satellogic.com service to be working. We are doing changes to have it public available. + + +Currently you have available these sources +------------------------------------------ + +- Memorytlesource: in memory storage. +- EtcTLESource: a uniq TLE is stored in `/etc/latest_tle` +- WSTLESource: It source is using the `TLE API. `_ + + +About HighAccuracyTLEPredictor +------------------------------ + +The default 'predictor' code is tunned to low CPU usage. (IE: a Satellite computer). The +error estimation is ~20 seconds. If you need more than that you can use the *HighAccuracyTLEPredictor* +passing `precise=True` to `get_predictor()`. + + +How to contribute +----------------- + +- Write pep8 complaint code. +- Wrap the code on 100 collumns. +- Always use a branch for each feature and Merge Proposals. +- Always run the tests before to push. (test implies pep8 validation) diff --git a/examples/iss.tle b/examples/iss.tle new file mode 100644 index 0000000..7e017f4 --- /dev/null +++ b/examples/iss.tle @@ -0,0 +1,3 @@ +ISS +1 25544U 98067A 17314.40254821 .00006490 00000-0 10525-3 0 9997 +2 25544 51.6429 29.4166 0004559 104.3372 354.3186 15.54111847 84492 diff --git a/orbit_predictor/__init__.py b/orbit_predictor/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/orbit_predictor/accuratepredictor.py b/orbit_predictor/accuratepredictor.py new file mode 100644 index 0000000..c046fd9 --- /dev/null +++ b/orbit_predictor/accuratepredictor.py @@ -0,0 +1,329 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +""" +Accurate Predictor +~~~~~~~~~~~~~~~~~~ + +Provides a faster an better predictor + +IMPORTANT!! +All calculations use radians, instead degrees +Be warned! + + + +Known Issues +~~~~~~~~~~~~ +In some cases not studied deeply, (because we don't have enough data) +ascending or descending point are not found and propagation fails. + +Code use some hacks to prevent multiple calculations, and function calls are +small as posible. + +Some stuff won't be trivial to understand, but comments and fixes are welcome + +""" +import datetime +import logging +import math +from datetime import timedelta +from math import degrees + +from sgp4 import ext, model +from sgp4.earth_gravity import wgs84 +from sgp4.io import twoline2rv +from sgp4.propagation import _gstime + +from orbit_predictor import coordinate_systems +from orbit_predictor.exceptions import NotReachable, PropagationError +from orbit_predictor.predictors import PredictedPass, TLEPredictor +from orbit_predictor.utils import lru_cache, reify + +logger = logging.getLogger(__name__) + +ONE_SECOND = datetime.timedelta(seconds=1) + +# Hack Zone be warned + + +@lru_cache(maxsize=365) +def jday_day(year, mon, day): + return (367.0 * year - + 7.0 * (year + ((mon + 9.0) // 12.0)) * 0.25 // 1.0 + + 275.0 * mon // 9.0 + + day + 1721013.5) + + +def jday(year, mon, day, hr, minute, sec): + base = jday_day(year, mon, day) + return base + ((sec / 60.0 + minute) / 60.0 + hr) / 24.0 + + +ext.jday = jday +model.jday = jday + +# finish hack zone + + +def round_datetime(dt): + return datetime.datetime(*dt.timetuple()[:6]) + + +class HighAccuracyTLEPredictor(TLEPredictor): + """A pass predictor with high accuracy on estimations""" + + def __init__(self, sate_id, source): + super(HighAccuracyTLEPredictor, self).__init__(sate_id, source) + + @reify + def tle(self): + return self.source.get_tle(self.sate_id, datetime.datetime.utcnow()) + + @reify + def propagator(self): + tle_line_1, tle_line_2 = self.tle.lines + return twoline2rv(tle_line_1, tle_line_2, wgs84) + + @reify + def mean_motion(self): + return self.propagator.no # this speed is in radians/minute + + @lru_cache(maxsize=3600 * 24 * 7) # Max cache, a week + def _propagate_only_position_ecef(self, timetuple): + """Return position in the given date using ECEF coordinate system.""" + position_eci, _ = self.propagator.propagate(*timetuple) + gmst = _gstime(jday(*timetuple)) + return coordinate_systems.eci_to_ecef(position_eci, gmst) + + def _propagate_ecef(self, when_utc): + """Return position and velocity in the given date using ECEF coordinate system.""" + timetuple = (when_utc.year, when_utc.month, when_utc.day, + when_utc.hour, when_utc.minute, when_utc.second) + + position_eci, velocity_eci = self.propagator.propagate(*timetuple) + gmst = _gstime(jday(*timetuple)) + position_ecef = coordinate_systems.eci_to_ecef(position_eci, gmst) + velocity_ecef = coordinate_systems.eci_to_ecef(velocity_eci, gmst) + return (position_ecef, velocity_ecef) + + def get_only_position(self, when_utc): + """Return a tuple in ECEF coordinate system + + Code is optimized, dont complain too much! + """ + timetuple = (when_utc.year, when_utc.month, when_utc.day, + when_utc.hour, when_utc.minute, when_utc.second) + return self._propagate_only_position_ecef(timetuple) + + def passes_over(self, location, when_utc, limit_date=None, max_elevation_gt=0, aos_at_dg=0): + return LocationPredictor(location, self, when_utc, limit_date, + max_elevation_gt, aos_at_dg) + + def get_next_pass(self, location, when_utc=None, max_elevation_gt=5, + aos_at_dg=0, limit_date=None): + """Implements same api as standard predictor""" + if when_utc is None: + when_utc = datetime.datetime.utcnow() + + for pass_ in self.passes_over(location, when_utc, limit_date, + max_elevation_gt=max_elevation_gt, + aos_at_dg=aos_at_dg): + return pass_ + else: + raise NotReachable('Propagation limit date exceded') + + +class AccuratePredictedPass: + + def __init__(self, aos, tca, los, max_elevation): + self.aos = round_datetime(aos) if aos is not None else None + self.tca = round_datetime(tca) + self.los = round_datetime(los) if los is not None else None + self.max_elevation = max_elevation + + @property + def valid(self): + return self.max_elevation > 0 and self.aos is not None and self.los is not None + + @reify + def max_elevation_deg(self): + return degrees(self.max_elevation) + + @reify + def duration(self): + return self.los - self.aos + + +class LocationPredictor(object): + """Predicts passes over a given location + Exposes an iterable interface + """ + + def __init__(self, location, propagator, start_date, limit_date=None, + max_elevation_gt=0, aos_at_dg=0): + self.location = location + self.propagator = propagator + self.start_date = start_date + self.limit_date = limit_date + + self.max_elevation_gt = math.radians(max([max_elevation_gt, aos_at_dg])) + self.aos_at = math.radians(aos_at_dg) + + def __iter__(self): + """Returns one pass each time""" + current_date = self.start_date + while True: + if self.is_ascending(current_date): + # we need a descending point + ascending_date = current_date + descending_date = self._find_nearest_descending(ascending_date) + pass_ = self._refine_pass(ascending_date, descending_date) + if pass_.valid: + yield self._build_predicted_pass(pass_) + + if self.limit_date is not None and current_date > self.limit_date: + break + + current_date = pass_.tca + self._orbit_step(0.6) + + else: + current_date = self._find_nearest_ascending(current_date) + + def _build_predicted_pass(self, accuratepass): + """Returns a classic predicted pass""" + tca_position = self.propagator.get_position(accuratepass.tca) + + return PredictedPass(self.location, self.propagator.sate_id, + max_elevation_deg=accuratepass.max_elevation_deg, + aos=accuratepass.aos, + los=accuratepass.los, + duration_s=accuratepass.duration.total_seconds(), + max_elevation_position=tca_position, + max_elevation_date=accuratepass.tca, + ) + + def _find_nearest_descending(self, ascending_date): + for candidate in self._sample_points(ascending_date): + if not self.is_ascending(candidate): + return candidate + else: + logger.error('Could not find a descending pass over %s start date: %s - TLE: %s', + self.location, ascending_date, self.propagator.tle) + raise PropagationError("Can not find an descending phase") + + def _find_nearest_ascending(self, descending_date): + for candidate in self._sample_points(descending_date): + if self.is_ascending(candidate): + return candidate + else: + logger.error('Could not find an ascending pass over %s start date: %s - TLE: %s', + self.location, descending_date, self.propagator.tle) + raise PropagationError('Can not find an ascending phase') + + def _sample_points(self, date): + """Helper method to found ascending or descending phases of elevation""" + start = date + end = date + self._orbit_step(0.99) + mid = self.midpoint(start, end) + mid_right = self.midpoint(mid, end) + mid_left = self.midpoint(start, mid) + + return [end, mid, mid_right, mid_left] + + def _refine_pass(self, ascending_date, descending_date): + tca = self._find_tca(ascending_date, descending_date) + elevation = self._elevation_at(tca) + + if elevation > self.max_elevation_gt: + aos = self._find_aos(tca) + los = self._find_los(tca) + else: + aos = los = None + + return AccuratePredictedPass(aos, tca, los, elevation) + + def _find_tca(self, ascending_date, descending_date): + while not self._precision_reached(ascending_date, descending_date): + midpoint = self.midpoint(ascending_date, descending_date) + if self.is_ascending(midpoint): + ascending_date = midpoint + else: + descending_date = midpoint + + return ascending_date + + def _precision_reached(self, start, end): + return end - start <= ONE_SECOND + + @staticmethod + def midpoint(start, end): + """Returns the midpoint between two dates""" + return start + (end - start) / 2 + + def _elevation_at(self, when_utc): + position = self.propagator.get_only_position(when_utc) + return self.location.elevation_for(position) + + def is_passing(self, when_utc): + """Returns a boolean indicating if satellite is actually visible""" + return bool(self._elevation_at(when_utc)) + + def is_ascending(self, when_utc): + """Check is elevation is ascending or descending on a given point""" + elevation = self._elevation_at(when_utc) + next_elevation = self._elevation_at(when_utc + ONE_SECOND) + return elevation <= next_elevation + + def _orbit_step(self, size): + """Returns a time step, that will make the satellite advance a given number of orbits""" + step_in_radians = size * 2 * math.pi + seconds = (step_in_radians / self.propagator.mean_motion) * 60 + return timedelta(seconds=seconds) + + def _find_aos(self, tca): + end = tca + start = tca - self._orbit_step(0.34) # On third of the orbit + elevation = self._elevation_at(start) + assert elevation < 0 + while not self._precision_reached(start, end): + midpoint = self.midpoint(start, end) + elevation = self._elevation_at(midpoint) + if elevation < self.aos_at: + start = midpoint + else: + end = midpoint + return end + + def _find_los(self, tca): + start = tca + end = tca + self._orbit_step(0.34) + while not self._precision_reached(start, end): + midpoint = self.midpoint(start, end) + elevation = self._elevation_at(midpoint) + + if elevation < self.aos_at: + end = midpoint + else: + start = midpoint + + return start diff --git a/orbit_predictor/coordinate_systems.py b/orbit_predictor/coordinate_systems.py new file mode 100644 index 0000000..a435641 --- /dev/null +++ b/orbit_predictor/coordinate_systems.py @@ -0,0 +1,153 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +from math import asin, atan, atan2, cos, degrees, pi, radians, sin, sqrt + + +def llh_to_ecef(lat_deg, lon_deg, h_km): + """ + Latitude is geodetic, height is above ellipsoid. Output is in km. + Formula from http://mathforum.org/library/drmath/view/51832.html + """ + f = 1 / 298.257224 + a = 6378.137 + lat_rad = radians(lat_deg) + lon_rad = radians(lon_deg) + cos_lat = cos(lat_rad) + sin_lat = sin(lat_rad) + C = 1 / sqrt(cos_lat ** 2 + (1 - f)**2 * sin_lat ** 2) + S = (1 - f) ** 2 * C + k1 = a * C + h_km + return (k1 * cos_lat * cos(lon_rad), k1 * cos_lat * sin(lon_rad), (a * S + h_km) * sin_lat) + + +# TODO: Same transformation as llh_to_ecef +def geodetic_to_ecef(lat, lon, height_km): + a = 6378.137 + b = 6356.7523142 + f = (a - b) / a + e2 = ((2 * f) - (f * f)) + normal = a / sqrt(1. - (e2 * (sin(lat) * sin(lat)))) + x = (normal + height_km) * cos(lat) * cos(lon) + y = (normal + height_km) * cos(lat) * sin(lon) + z = ((normal * (1. - e2)) + height_km) * sin(lat) + return x, y, z + + +def ecef_to_llh(ecef_km): + # WGS-84 ellipsoid parameters */ + a = 6378.1370 + b = 6356.752314 + + p = sqrt(ecef_km[0] ** 2 + ecef_km[1] ** 2) + thet = atan(ecef_km[2] * a / (p * b)) + esq = 1.0 - (b / a) ** 2 + epsq = (a / b) ** 2 - 1.0 + + lat = atan((ecef_km[2] + epsq * b * sin(thet) ** 3) / (p - esq * a * cos(thet) ** 3)) + lon = atan2(ecef_km[1], ecef_km[0]) + n = a * a / sqrt(a * a * cos(lat) ** 2 + b ** 2 * sin(lat) ** 2) + h = p / cos(lat) - n + + lat = degrees(lat) + lon = degrees(lon) + return lat, lon, h + + +def eci_to_ecef(eci_coords, gmst): + # ccar.colorado.edu/ASEN5070/handouts/coordsys.doc + # + # [X] [C -S 0][X] + # [Y] = [S C 0][Y] + # [Z]eci [0 0 1][Z]ecef + # + # + # Inverse: + # [X] [C S 0][X] + # [Y] = [-S C 0][Y] + # [Z]ecef [0 0 1][Z]e + sin_gmst = sin(gmst) + cos_gmst = cos(gmst) + eci_x, eci_y, eci_z = eci_coords + x = (eci_x * cos_gmst) + (eci_y * sin_gmst) + y = (eci_x * (-sin_gmst)) + (eci_y * cos_gmst) + z = eci_z + return x, y, z + + +def ecef_to_eci(eci_coords, gmst): + # ccar.colorado.edu/ASEN5070/handouts/coordsys.doc + # + # [X] [C -S 0][X] + # [Y] = [S C 0][Y] + # [Z]eci [0 0 1][Z]ecef + # + # + # Inverse: + # [X] [C S 0][X] + # [Y] = [-S C 0][Y] + # [Z]ecef [0 0 1][Z]e + x = (eci_coords[0] * cos(gmst)) - (eci_coords[1] * sin(gmst)) + y = (eci_coords[0] * (sin(gmst))) + (eci_coords[1] * cos(gmst)) + z = eci_coords[2] + return x, y, z + + +def horizon_to_az_elev(top_s, top_e, top_z): + range_sat = sqrt((top_s * top_s) + (top_e * top_e) + (top_z * top_z)) + elevation = asin(top_z / range_sat) + azimuth = atan2(-top_e, top_s) + pi + return azimuth, elevation + + +def to_horizon(observer_pos_lat_rad, observer_pos_long_rad, observer_pos_ecef, object_coords_ecef): + # http://www.celestrak.com/columns/v02n02/ + # TS Kelso's method, except I'm using ECF frame + # and he uses ECI. + + rx = object_coords_ecef[0] - observer_pos_ecef[0] + ry = object_coords_ecef[1] - observer_pos_ecef[1] + rz = object_coords_ecef[2] - observer_pos_ecef[2] + + sin_observer_lat = sin(observer_pos_lat_rad) + sin_observer_long = sin(observer_pos_long_rad) + cos_observer_lat = cos(observer_pos_lat_rad) + cos_observer_long = cos(observer_pos_long_rad) + + top_s = ((sin_observer_lat * cos_observer_long * rx) + + (sin_observer_lat * sin_observer_long * ry) - + (cos_observer_lat * rz)) + top_e = -sin_observer_long * rx + cos_observer_long * ry + top_z = ((cos_observer_lat * cos_observer_long * rx) + + (cos_observer_lat * sin_observer_long * ry) + + (sin_observer_lat * rz)) + + return top_s, top_e, top_z + + +def deg_to_dms(deg): + d = int(deg) + md = abs(deg - d) * 60 + m = int(md) + sd = (md - m) * 60 + return [d, m, sd] diff --git a/orbit_predictor/exceptions.py b/orbit_predictor/exceptions.py new file mode 100644 index 0000000..43218c8 --- /dev/null +++ b/orbit_predictor/exceptions.py @@ -0,0 +1,33 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +""" +Exceptions for orbit predictor +""" + + +class NotReachable(Exception): + """Raised when a pass is not reachable on a time window""" + + +class PropagationError(Exception): + """Raised when a calculation issue is found""" diff --git a/orbit_predictor/locations.py b/orbit_predictor/locations.py new file mode 100644 index 0000000..c128e08 --- /dev/null +++ b/orbit_predictor/locations.py @@ -0,0 +1,261 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import datetime +from math import asin, cos, degrees, radians, sin, sqrt + +from orbit_predictor import coordinate_systems +from orbit_predictor.utils import reify, sun_azimuth_elevation + +LIGHT_SPEED_KMS = 299792.458 + + +class Location(object): + def __init__(self, name, latitude_deg, longitude_deg, elevation_m): + self.name = name + self.latitude_deg = latitude_deg + self.longitude_deg = longitude_deg + self.elevation_m = elevation_m + self.position_ecef = coordinate_systems.geodetic_to_ecef( + radians(latitude_deg), + radians(longitude_deg), + elevation_m / 1000.) + self.position_llh = latitude_deg, longitude_deg, elevation_m + + def __eq__(self, other): + return all([issubclass(other.__class__, Location), + self.name == other.name, + self.latitude_deg == other.latitude_deg, + self.longitude_deg == other.longitude_deg, + self.elevation_m == other.elevation_m]) + + def __repr__(self): + return "".format(self.name) + + def __str__(self): + return self.name + + @reify + def latitude_rad(self): + return radians(self.latitude_deg) + + @reify + def longitude_rad(self): + return radians(self.longitude_deg) + + @reify + def _cached_elevation_calculation_data(self): + sin_lat, sin_long = sin(self.latitude_rad), sin(self.longitude_rad) + cos_lat, cos_long = cos(self.latitude_rad), cos(self.longitude_rad) + return (cos_lat * cos_long, + cos_lat * sin_long, + sin_lat) + + def sun_elevation_on_earth(self, when_utc=None): + """Return Sun elevation on Earth of location at when_utc.""" + if when_utc is None: + when_utc = datetime.datetime.utcnow() + _, elevation = sun_azimuth_elevation(self.latitude_deg, self.longitude_deg, when_utc) + return elevation + + def elevation_for(self, position): + """Returns elevation to given position in radians + + calculation is made inline to have better performance + """ + observer_pos_ecef = self.position_ecef + object_coords_ecef = position + + rx = object_coords_ecef[0] - observer_pos_ecef[0] + ry = object_coords_ecef[1] - observer_pos_ecef[1] + rz = object_coords_ecef[2] - observer_pos_ecef[2] + + a, b, c = self._cached_elevation_calculation_data + + top_z = (a * rx) + (b * ry) + (c * rz) + + range_sat = sqrt((rx * rx) + (ry * ry) + (rz * rz)) + + return asin(top_z / range_sat) + + def get_azimuth_elev(self, position): + """Return azimuth and elevation of position_ecef from the current Location instance.""" + + top = coordinate_systems.to_horizon(self.latitude_rad, self.longitude_rad, + self.position_ecef, position.position_ecef) + + return coordinate_systems.horizon_to_az_elev(*top) + + def get_azimuth_elev_deg(self, position): + """Idem that get_azimuth_elev() but using degrees.""" + az, el = self.get_azimuth_elev(position) + return degrees(az), degrees(el) + + def is_visible(self, position, elevation=0): + """Return True if the Satellite if visible from the current instance.""" + _, elev_deg = self.get_azimuth_elev_deg(position) + return elev_deg >= elevation + + def slant_range_km(self, position_ecef): + """Distance to the satellite in straight line""" + pos = position_ecef + loc = self.position_ecef + return sqrt((pos[0]-loc[0])**2 + (pos[1]-loc[1])**2 + (pos[2]-loc[2])**2) + + def slant_range_velocity_kms(self, position): + """Velocity the satellite from location's point of view""" + pos = position.position_ecef + vel = position.velocity_ecef + + current_range = self.slant_range_km(pos) + next_pos = (pos[0]+vel[0], pos[1]+vel[1], pos[2]+vel[2]) + next_range = self.slant_range_km(next_pos) + + return next_range - current_range + + def doppler_factor(self, position): + """Doppler effect factor relative to 1""" + range_rate = self.slant_range_velocity_kms(position) + return 1. + (range_rate / LIGHT_SPEED_KMS) + + +# A hardcoded list of locations. Some of them are satellite groundstations or HAM +AFRICA1 = Location( + "AFRICA1", latitude_deg=-4.2937711, longitude_deg=15.4493049, elevation_m=266.00) +AFRICA2 = Location( + "AFRICA2", latitude_deg=-19.9243839, longitude_deg=23.439418, elevation_m=939.12) +AFRICA3 = Location( + "AFRICA3", latitude_deg=-26.0317764, longitude_deg=28.254681, elevation_m=1617.62) +AFRICA4 = Location( + "AFRICA4", latitude_deg=0.3979327, longitude_deg=32.5021788, elevation_m=1165.38) +AFRICA5 = Location( + "AFRICA5", latitude_deg=-1.2960418, longitude_deg=36.9340893, elevation_m=1599.74) +AMERICA1 = Location( + "AMERICA1", latitude_deg=40.6599903, longitude_deg=-74.1713736, elevation_m=10.46) +AMERICA10 = Location( + "AMERICA10", latitude_deg=34.0863943, longitude_deg=-118.0329261, elevation_m=88.67) +AMERICA11 = Location( + "AMERICA11", latitude_deg=37.9916467, longitude_deg=-122.0559013, elevation_m=6.53) +AMERICA2 = Location( + "AMERICA2", latitude_deg=38.9054965, longitude_deg=-77.0230685, elevation_m=25.25) +AMERICA3 = Location( + "AMERICA3", latitude_deg=33.7800684, longitude_deg=-84.5208486, elevation_m=245.74) +AMERICA4 = Location( + "AMERICA4", latitude_deg=29.9414947, longitude_deg=-90.0633866, elevation_m=3.64) +AMERICA5 = Location( + "AMERICA5", latitude_deg=29.9865571, longitude_deg=-95.3423456, elevation_m=29.35) +AMERICA6 = Location( + "AMERICA6", latitude_deg=19.4361691, longitude_deg=-99.0719249, elevation_m=2224.95) +AMERICA7 = Location( + "AMERICA7", latitude_deg=20.5216683, longitude_deg=-103.310728, elevation_m=1530.03) +AMERICA8 = Location( + "AMERICA8", latitude_deg=35.0401972, longitude_deg=-106.6090026, elevation_m=1619.52) +AMERICA9 = Location( + "AMERICA9", latitude_deg=33.6928889, longitude_deg=-112.078808, elevation_m=450.69) +ARG = Location("ARG", latitude_deg=-31.2884, longitude_deg=-64.2032868, elevation_m=492.96) +ASIA1 = Location("ASIA1", latitude_deg=32.0092853, longitude_deg=34.8945777, elevation_m=38.03) +ASIA10 = Location("ASIA10", latitude_deg=23.8450823, longitude_deg=90.4016501, elevation_m=11.71) +ASIA11 = Location("ASIA11", latitude_deg=16.9069935, longitude_deg=96.1342117, elevation_m=24.81) +ASIA12 = Location("ASIA12", latitude_deg=13.9125333, longitude_deg=100.6068365, elevation_m=6.22) +ASIA13 = Location("ASIA13", latitude_deg=21.2137962, longitude_deg=105.805638, elevation_m=12.45) +ASIA14 = Location("ASIA14", latitude_deg=23.3924882, longitude_deg=113.2990193, elevation_m=5.97) +ASIA15 = Location("ASIA15", latitude_deg=31.7392443, longitude_deg=118.8733768, elevation_m=13.34) +ASIA16 = Location("ASIA16", latitude_deg=37.5602464, longitude_deg=126.7909134, elevation_m=20.35) +ASIA17 = Location("ASIA17", latitude_deg=33.5846715, longitude_deg=130.4510901, elevation_m=5.67) +ASIA18 = Location("ASIA18", latitude_deg=34.7854932, longitude_deg=135.4384313, elevation_m=10.01) +ASIA19 = Location("ASIA19", latitude_deg=35.7120096, longitude_deg=139.4033569, elevation_m=91.00) +ASIA2 = Location("ASIA2", latitude_deg=24.5536664, longitude_deg=39.7053953, elevation_m=638.85) +ASIA3 = Location("ASIA3", latitude_deg=33.2621459, longitude_deg=44.234124, elevation_m=36.05) +ASIA4 = Location("ASIA4", latitude_deg=35.6883245, longitude_deg=51.3143664, elevation_m=1185.27) +ASIA5 = Location("ASIA5", latitude_deg=36.2320987, longitude_deg=59.6430435, elevation_m=996.20) +ASIA6 = Location("ASIA6", latitude_deg=31.628871, longitude_deg=65.7371749, elevation_m=1014.35) +ASIA7 = Location("ASIA7", latitude_deg=33.6187486, longitude_deg=73.0960301, elevation_m=502.00) +ASIA8 = Location("ASIA8", latitude_deg=28.5543983, longitude_deg=77.086455, elevation_m=226.08) +ASIA9 = Location("ASIA9", latitude_deg=12.950055, longitude_deg=77.66856, elevation_m=889.07) +AUSTRALIA1 = Location( + "AUSTRALIA1", latitude_deg=-31.9170947, longitude_deg=115.970206, elevation_m=12.73) +AUSTRALIA2 = Location( + "AUSTRALIA2", latitude_deg=-17.8184141, longitude_deg=122.2364966, elevation_m=28.67) +AUSTRALIA5 = Location( + "AUSTRALIA5", latitude_deg=-34.7794086, longitude_deg=138.6370729, elevation_m=22.43) +AUSTRALIA6 = Location( + "AUSTRALIA6", latitude_deg=-36.7103328, longitude_deg=144.3303179, elevation_m=197.96) +AUSTRALIA7 = Location( + "AUSTRALIA7", latitude_deg=-34.0096484, longitude_deg=150.6926073, elevation_m=74.77) +BA1 = Location("BA1", latitude_deg=-34.5561944, longitude_deg=-58.41368, elevation_m=7.02) +CHILE = Location("CHILE", latitude_deg=-33.3631552, longitude_deg=-70.7904123, elevation_m=477.00) +EASTER_ISLAND = Location( + "EASTER_ISLAND", latitude_deg=-27.0578009, longitude_deg=-109.3817317, elevation_m=61.69) +EUROPA1 = Location("EUROPA1", latitude_deg=41.2486859, longitude_deg=-8.6813677, elevation_m=56.44) +EUROPA10 = Location("EUROPA10", latitude_deg=45.7274069, longitude_deg=65.37, elevation_m=122.98) +EUROPA11 = Location( + "EUROPA11", latitude_deg=45.4452575, longitude_deg=9.2767394, elevation_m=106.02) +EUROPA12 = Location( + "EUROPA12", latitude_deg=48.3534778, longitude_deg=11.7864782, elevation_m=447.39) +EUROPA13 = Location( + "EUROPA13", latitude_deg=42.4310529, longitude_deg=14.1828016, elevation_m=9.81) +EUROPA14 = Location( + "EUROPA14", latitude_deg=41.1150865, longitude_deg=16.8624173, elevation_m=13.89) +EUROPA15 = Location("EUROPA15", latitude_deg=37.9364517, longitude_deg=23.94452, elevation_m=80.00) +EUROPA16 = Location( + "EUROPA16", latitude_deg=38.2925088, longitude_deg=27.1556125, elevation_m=119.68) +EUROPA17 = Location( + "EUROPA17", latitude_deg=35.1544144, longitude_deg=33.3585865, elevation_m=172.25) +EUROPA3 = Location("EUROPA3", latitude_deg=37.4189722, longitude_deg=-5.8929429, elevation_m=27.52) +EUROPA5 = Location( + "EUROPA5", latitude_deg=40.4915238, longitude_deg=-3.5677712, elevation_m=597.39) +EUROPA7 = Location("EUROPA7", latitude_deg=39.4892396, longitude_deg=-0.4819177, elevation_m=60.64) +EUROPA9 = Location("EUROPA9", latitude_deg=49.0067717, longitude_deg=2.5529958, elevation_m=102.37) +MADAGASCAR1 = Location( + "MADAGASCAR1", latitude_deg=15.4967687, longitude_deg=44.2171958, elevation_m=2186.00) +MADAGASCAR2 = Location( + "MADAGASCAR2", latitude_deg=-18.7825536, longitude_deg=47.4800904, elevation_m=1260.62) +NZ1 = Location("NZ1", latitude_deg=-44.7149065, longitude_deg=169.2468643, elevation_m=339.58) +NZ2 = Location("NZ2", latitude_deg=-36.5886632, longitude_deg=174.8717244, elevation_m=0.00) +RIO = Location("RIO", latitude_deg=-22.910590, longitude_deg=-43.188958, elevation_m=16.92) +USA = Location("USA", latitude_deg=40.24, longitude_deg=-101.9, elevation_m=1100) +australia = Location('australia', latitude_deg=-25.1, longitude_deg=134.5, elevation_m=290) +brazil = Location("brazil", latitude_deg=-11.2, longitude_deg=-54.66, elevation_m=310) +blq_leafline = Location('blq_leafline', latitude_deg=45.59, longitude_deg=9.361, elevation_m=194) +brc = Location("BRC", latitude_deg=-41.1424, longitude_deg=-71.1530, elevation_m=800) +cba_conae = Location( + 'cba_conae', latitude_deg=-31.52407, longitude_deg=-64.46352, elevation_m=730) +central_america = Location( + "central_america", latitude_deg=11.17, longitude_deg=-87.23, elevation_m=310) +central_argentina = Location( + 'central_argentina', latitude_deg=-35.75, longitude_deg=-63.9, elevation_m=133) +china = Location('china', latitude_deg=35.4, longitude_deg=110, elevation_m=1000) +eastern_russia = Location('eastern_russia', latitude_deg=66, longitude_deg=145, elevation_m=650) +france = Location('france', latitude_deg=46.4, longitude_deg=2.75, elevation_m=300) +germany = Location("ALEMANIA", latitude_deg=52.515083, longitude_deg=13.323723, elevation_m=30) +india = Location('india', latitude_deg=23.5, longitude_deg=78.5, elevation_m=550) +moscu = Location('moscu', latitude_deg=55.7, longitude_deg=37.5, elevation_m=137) +niger = Location('niger', latitude_deg=20, longitude_deg=12.5, elevation_m=430) +riogrande = Location("RIOGRANDE", latitude_deg=-53.8, longitude_deg=-67.75, elevation_m=30) +svalbard = Location("SVALBARD", latitude_deg=78.229772, longitude_deg=15.407786, elevation_m=501) +# HAM +DK3WN = Location("DK3WN", latitude_deg=49.7317, longitude_deg=8.9548, elevation_m=138) +JA6PL = Location("JA6PL", latitude_deg=33.8685, longitude_deg=130.719, elevation_m=10) +LU4EOU = Location("LU4EOU", latitude_deg=-38.7243008, longitude_deg=-62.2449124, elevation_m=10) +# tortu1 == tortu2 == LU1CGB +tortu1 = Location("TORTU1", latitude_deg=-34.4791, longitude_deg=-58.7916, elevation_m=10) +tortu2 = Location("TORTU2", latitude_deg=-34.4791, longitude_deg=-58.7916, elevation_m=10) diff --git a/orbit_predictor/predictors.py b/orbit_predictor/predictors.py new file mode 100644 index 0000000..7e872d9 --- /dev/null +++ b/orbit_predictor/predictors.py @@ -0,0 +1,246 @@ +# -*- coding: utf-8 -*- +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +import datetime +import logging +import warnings +from collections import namedtuple +from math import acos, degrees + +from sgp4.earth_gravity import wgs84 +from sgp4.ext import jday +from sgp4.io import twoline2rv +from sgp4.propagation import _gstime + +from orbit_predictor import coordinate_systems +from orbit_predictor.exceptions import NotReachable +from orbit_predictor.locations import Location +from orbit_predictor.utils import ( + cross_product, + dot_product, + reify, + vector_diff, + vector_norm +) + +logger = logging.getLogger(__name__) + + +class Position(namedtuple( + "Position", ['when_utc', 'position_ecef', 'velocity_ecef', 'error_estimate'])): + + @reify + def position_llh(self): + return coordinate_systems.ecef_to_llh(self.position_ecef) + + +class PredictedPass(object): + def __init__(self, location, sate_id, + max_elevation_deg, + aos, los, duration_s, + max_elevation_position=None, + max_elevation_date=None): + self.location = location + self.sate_id = sate_id + self.max_elevation_position = max_elevation_position + self.max_elevation_date = max_elevation_date + self.max_elevation_deg = max_elevation_deg + self.aos = aos + self.los = los + self.duration_s = duration_s + + @property + def midpoint(self): + """Returns a datetime of the midpoint of the pass""" + return self.aos + (self.los - self.aos) / 2 + + def __repr__(self): + return "".format(self.sate_id, self.location, self.aos) + + def __eq__(self, other): + return all([issubclass(other.__class__, PredictedPass), + self.location == other.location, + self.sate_id == other.sate_id, + self.max_elevation_position == other.max_elevation_position, + self.max_elevation_date == other.max_elevation_date, + self.max_elevation_deg == other.max_elevation_deg, + self.aos == other.aos, + self.los == other.los, + self.duration_s == other.duration_s]) + + def get_off_nadir_angle(self): + warnings.warn("This method is deprecated!", DeprecationWarning) + return self.off_nadir_deg + + @reify + def off_nadir_deg(self): + """Computes off-nadir angle calculation + + Given satellite position ``sate_pos``, velocity ``sate_vel``, and + location ``target`` in a common frame, off-nadir angle ``off_nadir_angle`` + is given by: + t2b = sate_pos - target + cos(off_nadir_angle) = (target · t2b) # Vectorial dot product + _____________________ + || target || || t2b|| + + Sign for the rotation is calculated this way + + cross = target ⨯ sate_pos + sign = cross · sate_vel + ____________________ + | cross · sate_vel | + """ + sate_pos = self.max_elevation_position.position_ecef + sate_vel = self.max_elevation_position.velocity_ecef + target = self.location.position_ecef + t2b = vector_diff(sate_pos, target) + angle = acos( + dot_product(target, t2b) / (vector_norm(target) * vector_norm(t2b)) + ) + + cross = cross_product(target, sate_pos) + dot = dot_product(cross, sate_vel) + try: + sign = dot / abs(dot) + except ZeroDivisionError: + sign = 1 + + return degrees(angle) * sign + + +class Predictor(object): + + def __init__(self, source, sate_id): + self.source = source + self.sate_id = sate_id + + def get_position(self, when_utc=None): + raise NotImplementedError("You have to implement it!") + + +class TLEPredictor(Predictor): + + def __init__(self, sate_id, source): + super(TLEPredictor, self).__init__(source, sate_id) + self._iterations = 0 + + def _propagate_eci(self, when_utc=None): + """Return position and velocity in the given date using ECI coordinate system.""" + tle = self.source.get_tle(self.sate_id, when_utc) + logger.debug("Propagating using ECI. sate_id: %s, when_utc: %s, tle: %s", + self.sate_id, when_utc, tle) + tle_line_1, tle_line_2 = tle.lines + sgp4_sate = twoline2rv(tle_line_1, tle_line_2, wgs84) + timetuple = when_utc.timetuple()[:6] + position_eci, velocity_eci = sgp4_sate.propagate(*timetuple) + return (position_eci, velocity_eci) + + def _gstime_from_datetime(self, when_utc): + timetuple = when_utc.timetuple()[:6] + return _gstime(jday(*timetuple)) + + def _propagate_ecef(self, when_utc=None): + """Return position and velocity in the given date using ECEF coordinate system.""" + position_eci, velocity_eci = self._propagate_eci(when_utc) + gmst = self._gstime_from_datetime(when_utc) + position_ecef = coordinate_systems.eci_to_ecef(position_eci, gmst) + velocity_ecef = coordinate_systems.eci_to_ecef(velocity_eci, gmst) + return (position_ecef, velocity_ecef) + + def get_position(self, when_utc=None): + """Return a Position namedtuple in ECEF coordinate system""" + if when_utc is None: + when_utc = datetime.datetime.utcnow() + position_ecef, velocity_ecef = self._propagate_ecef(when_utc) + + return Position(when_utc=when_utc, position_ecef=position_ecef, + velocity_ecef=velocity_ecef, error_estimate=None) + + def get_next_pass(self, location, when_utc=None, max_elevation_gt=5, + aos_at_dg=0, limit_date=None): + """Return a PredictedPass instance with the data of the next pass over the given location + + locattion_llh: point on Earth we want to see from the satellite. + when_utc: datetime UTC. + max_elevation_gt: filter passings with max_elevation under it. + aos_at_dg: This is if we want to start the pass at a specific elevation. + """ + if when_utc is None: + when_utc = datetime.datetime.utcnow() + if max_elevation_gt < aos_at_dg: + max_elevation_gt = aos_at_dg + pass_ = self._get_next_pass(location, when_utc, aos_at_dg, limit_date) + while pass_.max_elevation_deg < max_elevation_gt: + pass_ = self._get_next_pass( + location, pass_.los, aos_at_dg, limit_date) # when_utc is changed! + return pass_ + + def _get_next_pass(self, location, when_utc, aos_at_dg=0, limit_date=None): + if not isinstance(location, Location): + raise TypeError("location must be a Location instance") + + pass_ = PredictedPass(location=location, sate_id=self.sate_id, aos=None, los=None, + max_elevation_date=None, max_elevation_position=None, + max_elevation_deg=0, duration_s=0) + + seconds = 0 + self._iterations = 0 + while True: + # to optimize the increment in seconds must be inverse proportional to + # the distance of 0 elevation + date = when_utc + datetime.timedelta(seconds=seconds) + + if limit_date is not None and date > limit_date: + raise NotReachable('Propagation limit date exceded') + + elev_pos = self.get_position(date) + _, elev = location.get_azimuth_elev(elev_pos) + elev_deg = degrees(elev) + + if elev_deg > pass_.max_elevation_deg: + pass_.max_elevation_position = elev_pos + pass_.max_elevation_date = date + pass_.max_elevation_deg = elev_deg + + if elev_deg > aos_at_dg and pass_.aos is None: + pass_.aos = date + if pass_.aos and elev_deg < aos_at_dg: + pass_.los = date + pass_.duration_s = (pass_.los - pass_.aos).total_seconds() + break + + if elev_deg < -2: + delta_s = abs(elev_deg) * 15 + 10 + else: + delta_s = 20 + + seconds += delta_s + self._iterations += 1 + + return pass_ + + +class GPSPredictor(Predictor): + pass diff --git a/orbit_predictor/sources.py b/orbit_predictor/sources.py new file mode 100644 index 0000000..2e047c1 --- /dev/null +++ b/orbit_predictor/sources.py @@ -0,0 +1,165 @@ +# -*- coding: utf-8 -*- +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import logging +from collections import defaultdict, namedtuple + +import requests + +from orbit_predictor.accuratepredictor import HighAccuracyTLEPredictor +from orbit_predictor.predictors import TLEPredictor + +try: + from urllib import parse as urlparse + from urllib.parse import urlencode +except ImportError: + # pytyhon2 support. + import urlparse # NOQA + from urllib import urlencode # NOQA + +logger = logging.getLogger(__name__) + +TLE = namedtuple('TLE', + ['sate_id', 'lines', 'date']) + + +class GPSSource(object): + def get_position_ecef(self, sate_id, when_utc): + raise NotImplementedError("You have to implement it.") + + +class TLESource(object): + + def add_tle(self, sate_id, tle, epoch): + raise NotImplementedError("You have to implement it.") + + def _get_tle(self, sate_id, date): + raise NotImplementedError("You have to implement it.") + + def get_tle(self, sate_id, date): + logger.debug("searching a TLE for %s, date: %s", sate_id, date) + lines = self._get_tle(sate_id, date) + return TLE(sate_id=sate_id, date=date, lines=lines) + + def get_predictor(self, sate_id, precise=False): + """Return a Predictor instance using the current storage.""" + if precise: + return HighAccuracyTLEPredictor(sate_id, self) + + return TLEPredictor(sate_id, self) + + +class MemoryTLESource(TLESource): + def __init__(self): + self.tles = defaultdict(lambda: set()) + + def add_tle(self, sate_id, tle, epoch): + self.tles[sate_id].add((epoch, tle)) + + def _get_tle(self, sate_id, date): + candidates = self.tles[sate_id] + + winner = None + winner_dt = float("inf") + + for epoch, candidate in candidates: + c_dt = abs((epoch - date).total_seconds()) + if c_dt < winner_dt: + winner = candidate + winner_dt = c_dt + + if winner is None: + raise LookupError("no tles in storage") + + return winner + + +class EtcTLESource(TLESource): + def __init__(self, filename="/etc/latest_tle"): + self.filename = filename + + def add_tle(self, sate_id, tle, epoch): + with open(self.filename, "w") as fd: + fd.write(sate_id + "\n") + for l in tle: + fd.write(l + "\n") + + def _get_tle(self, sate_id, date): + with open(self.filename) as fd: + data = fd.read() + lines = data.split("\n") + if not lines[0] == sate_id: + raise LookupError("Stored satellite id not") + return tuple(lines[1:3]) + + +class WSTLESource(TLESource): + + def __init__(self, url): + self.url = url + self.cache = MemoryTLESource() + + def add_tle(self, *args): + raise ValueError("You can't add TLEs. The service has his own update task.") + + def _get_tle(self, sate_id, date): + # first lookup on cache + try: + lines_from_cache = self.cache._get_tle(sate_id, date) + except LookupError: + pass + else: + return lines_from_cache + + lines = self.get_tle_for_date(sate_id, date) + # save on cache + self.cache.add_tle(sate_id, lines, date) + return lines + + def get_last_update(self, sate_id): + return self._fetch_tle("api/tle/last/", sate_id) + + def get_tle_for_date(self, sate_id, date): + return self._fetch_tle("api/tle/closest/", sate_id, date) + + def _fetch_tle(self, path, sate_id, date=None): + url = urlparse.urljoin(self.url, path) + url = urlparse.urlparse(url) + qargs = {'satellite_number': sate_id} + if date is not None: + date_str = date.strftime("%Y-%m-%d") + qargs['date'] = date_str + + query_string = urlencode(qargs) + url = urlparse.urlunsplit((url.scheme, url.netloc, url.path, query_string, url.fragment)) + headers = {'user-agent': 'orbit-predictor', 'Accept': 'application/json'} + try: + response = requests.get(url, headers=headers) + except requests.exceptions.RequestException as error: + logger.error("Exception requesting TLE: %s", error) + raise + if response.ok and 'lines' in response.json(): + lines = tuple(response.json()['lines']) + return lines + else: + raise ValueError("Error requesting TLE: %s", response.text) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py new file mode 100644 index 0000000..142edd2 --- /dev/null +++ b/orbit_predictor/utils.py @@ -0,0 +1,208 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import functools +from collections import namedtuple +from datetime import datetime +from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt + +try: + from functools import lru_cache +except ImportError: + class lru_cache(object): + """dummy function for python 2""" + def __init__(self, *args, **kwargs): + pass + + def __call__(self, f): + return f + +# This function was ported from its Matlab equivalent here: +# http://www.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation + +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) + + +cos_d = compose(cos, radians) +sin_d = compose(sin, radians) +atan2_d = compose(degrees, atan2) +asin_d = compose(degrees, asin) + + +AzimuthElevation = namedtuple('AzimuthElevation', 'azimuth elevation') + + +def euclidean_distance(*components): + """Returns the norm of a vector""" + return sqrt(sum(c**2 for c in components)) + + +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)) + + +def vector_diff(a, b): + """Computes difference between two vectors""" + return tuple((ai - bi) for ai, bi in zip(a, b)) + + +def cross_product(a, b): + """Computes cross product between two vectors""" + return ( + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0] + ) + + +def vector_norm(a): + """Returns the norm of a vector""" + return euclidean_distance(*a) + + +def sun_azimuth_elevation(latitude_deg, longitude_deg, when=None): + """ + Return (azimuth, elevation) of the Sun at ground point + + :param latitude_deg: a float number representing latitude on degrees + :param longitude_deg: a float number representing longitude on degrees + :param when: a ``datetime.datetime`` object in utc, if not provided, + utcnow() is used + :returns: an ``AzimuthElevation`` namedtuple + """ + if when is None: + when = datetime.utcnow() + + utc_time_tuple = when.timetuple() + jd = juliandate(utc_time_tuple) + date = jd - DECEMBER_31TH_1999_MIDNIGHT_JD + + w = 282.9404 + 4.70935e-5 * date # longitude of perihelion degrees + eccentricity = 0.016709 - 1.151e-9 * date # eccentricity + M = (356.0470 + 0.9856002585 * date) % 360 # mean anomaly degrees + L = w + M # Sun's mean longitude degrees + oblecl = 23.4393 - 3.563e-7 * date # Sun's obliquity of the ecliptic + + # auxiliary angle + auxiliary_angle = M + degrees(eccentricity * sin_d(M) * (1 + eccentricity * cos_d(M))) + + # rectangular coordinates in the plane of the ecliptic (x axis toward perhilion) + x = cos_d(auxiliary_angle) - eccentricity + y = sin_d(auxiliary_angle) * sqrt(1 - eccentricity**2) + + # find the distance and true anomaly + r = euclidean_distance(x, y) + v = atan2_d(y, x) + + # find the longitude of the sun + sun_lon = v + w + + # compute the ecliptic rectangular coordinates + xeclip = r * cos_d(sun_lon) + yeclip = r * sin_d(sun_lon) + zeclip = 0.0 + + # rotate these coordinates to equitorial rectangular coordinates + xequat = xeclip + yequat = yeclip * cos_d(oblecl) + zeclip * sin_d(oblecl) + zequat = yeclip * sin_d(23.4406) + zeclip * cos_d(oblecl) + + # convert equatorial rectangular coordinates to RA and Decl: + r = euclidean_distance(xequat, yequat, zequat) + RA = atan2_d(yequat, xequat) + delta = asin_d(zequat/r) + + # Following the RA DEC to Az Alt conversion sequence explained here: + # http://www.stargazing.net/kepler/altaz.html + + sidereal = sidereal_time(utc_time_tuple, longitude_deg, L) + + # Replace RA with hour angle HA + HA = sidereal * 15 - RA + + # convert to rectangular coordinate system + x = cos_d(HA) * cos_d(delta) + y = sin_d(HA) * cos_d(delta) + z = sin_d(delta) + + # rotate this along an axis going east-west. + xhor = x * cos_d(90 - latitude_deg) - z * sin_d(90 - latitude_deg) + yhor = y + zhor = x * sin_d(90 - latitude_deg) + z * cos_d(90 - latitude_deg) + + # Find the h and AZ + azimuth = atan2_d(yhor, xhor) + 180 + elevation = asin_d(zhor) + + return AzimuthElevation(azimuth, elevation) + + +def juliandate(utc_tuple): + year, month, day, hour, minute, sec = utc_tuple[:6] + if month <= 2: + year -= 1 + month += 12 + + return (floor(365.25*(year + 4716.0)) + floor(30.6001*(month+1.0)) + 2.0 - + floor(year / 100.0) + floor(floor(year / 100.0) / 4.0) + day - 1524.5 + + (hour + minute / 60.0 + sec / 3600.0) / 24.0) + + +def sidereal_time(utc_tuple, local_lon, sun_lon): + # Find the J2000 value + # J2000 = jd - 2451545.0; + UTH = utc_tuple.tm_hour + utc_tuple.tm_min / 60.0 + utc_tuple.tm_sec / 3600.0 + + # Calculate local siderial time + GMST0 = ((sun_lon + 180) % 360) / 15 + return GMST0 + UTH + local_lon / 15 + + +class reify(object): + """ + Use as a class method decorator. It operates almost exactly like the + Python ``@property`` decorator, but it puts the result of the method it + decorates into the instance dict after the first call, effectively + replacing the function it decorates with an instance variable. It is, in + Python parlance, a non-data descriptor. + + Taken from: http://docs.pylonsproject.org/projects/pyramid/en/latest/api/decorator.html + """ + + def __init__(self, wrapped): + self.wrapped = wrapped + functools.update_wrapper(self, wrapped) + + def __get__(self, inst, objtype=None): + if inst is None: + return self + val = self.wrapped(inst) + setattr(inst, self.wrapped.__name__, val) + return val diff --git a/orbit_predictor/version.py b/orbit_predictor/version.py new file mode 100644 index 0000000..21d456e --- /dev/null +++ b/orbit_predictor/version.py @@ -0,0 +1,2 @@ +# https://www.python.org/dev/peps/pep-0440/ +__version__ = '1.8.4' diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..766ac8f --- /dev/null +++ b/requirements.txt @@ -0,0 +1,9 @@ +flake8 +hypothesis +hypothesis[datetime] +logassert +mock +pytest +pytest-cov +requests>=2.9.1 +sgp4>=1.4 diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..ef28ae2 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,11 @@ +[flake8] +max-line-length = 99 +exclude = + .git, + __pycache__, + .ropeproject, + .fades + +[isort] +line_length = 79 +multi_line_output = 3 diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..219e273 --- /dev/null +++ b/setup.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python3 +import os.path +from setuptools import setup, find_packages + +# Copyright 2017 Satellogic SA. + + +# https://packaging.python.org/guides/single-sourcing-package-version/ +version = {} +with open(os.path.join("orbit_predictor", "version.py")) as fp: + exec(fp.read(), version) + + +setup( + name='orbit-predictor', + version=version["__version__"], + author='Satellogic SA', + author_email='oss@satellogic.com', + description='Python library to propagate satellite orbits.', + long_description=open('README.rst').read(), + packages=find_packages(exclude=["tests"]), + license="MIT", + url='https://github.com/satellogic/orbit-predictor', + # Keywords to get found easily on PyPI results,etc. + keywords="orbit, sgp4, TLE, space, satellites", + classifiers=[ + 'Development Status :: 5 - Production/Stable', + 'Environment :: Console', + 'Intended Audience :: Science/Research', + 'Operating System :: OS Independent', + 'Programming Language :: Python', + 'Topic :: Software Development :: Libraries :: Python Modules', + 'Topic :: Utilities', + 'License :: OSI Approved :: MIT License' + 'Programming Language :: Python :: 2.7', + 'Programming Language :: Python :: 3', + ], + install_requires=[ + 'sgp4', + 'requests', + ], + tests_require=['logassert', 'flake8', 'hypothesis', 'mock', 'hypothesis[datetime]'], + extras_require={ + "dev": [ + "hypothesis", + "flake8", + "hypothesis[datetime]", + "mock", + "logassert", + ] + } +) diff --git a/test b/test new file mode 100755 index 0000000..c1524b9 --- /dev/null +++ b/test @@ -0,0 +1,10 @@ +#!/bin/bash +# +# Copyright 2017 Satellogic +set -e +FADES='fades -r requirements.txt' +FADES2='fades -r requirements.txt -p python2 -d mock' + +$FADES -x flake8 +$FADES -x python3 -m unittest discover -v -b tests +$FADES2 -x python -m unittest discover -v -b tests diff --git a/tests/test_accurate_predictor.py b/tests/test_accurate_predictor.py new file mode 100644 index 0000000..8816d8e --- /dev/null +++ b/tests/test_accurate_predictor.py @@ -0,0 +1,281 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import time +from datetime import datetime, timedelta +from unittest import TestCase + +import logassert +from hypothesis import example, given, settings +from hypothesis.extra.datetime import datetimes +from hypothesis.strategies import floats, tuples + +from orbit_predictor.accuratepredictor import ( + ONE_SECOND, + HighAccuracyTLEPredictor +) +from orbit_predictor.exceptions import PropagationError +from orbit_predictor.locations import Location, svalbard, tortu1 +from orbit_predictor.predictors import TLEPredictor +from orbit_predictor.sources import MemoryTLESource + +try: + from unittest import mock # py3 +except ImportError: + import mock # py2 + + +SATE_ID = '41558U' # newsat 1 +LINES = ( + '1 41558U 16033C 17065.21129769 .00002236 00000-0 88307-4 0 9995', + '2 41558 97.4729 144.7611 0014207 16.2820 343.8872 15.26500433 42718', +) + +BUGSAT_SATE_ID = 'BUGSAT-1' +BUGSAT1_TLE_LINES = ( + "1 40014U 14033E 14294.41438078 .00003468 00000-0 34565-3 0 3930", + "2 40014 97.9781 190.6418 0032692 299.0467 60.7524 14.91878099 18425") + + +class AccuratePredictorTests(TestCase): + + def setUp(self): + # Source + self.db = MemoryTLESource() + self.start = datetime(2017, 3, 6, 7, 51) + self.db.add_tle(SATE_ID, LINES, self.start) + # Predictor + self.predictor = HighAccuracyTLEPredictor(SATE_ID, self.db) + self.old_predictor = TLEPredictor(SATE_ID, self.db) + self.end = self.start + timedelta(days=5) + + def all_passes_old_predictor(self, location, start, end): + while True: + try: + pass_ = self.old_predictor.get_next_pass(location, + when_utc=start, + limit_date=end) + start = pass_.los + yield pass_ + except Exception: + break + + def assertEqualOrGreaterPassesAmount(self, location): + """Compare propagators and check no passes lost and no performance degradation""" + t0 = time.time() + old_passes = list( + self.all_passes_old_predictor(location, self.start, self.end) + ) + t1 = time.time() + predicted_passes = list( + self.predictor.passes_over(location, self.start, self.end) + ) + t2 = time.time() + self.assertGreaterEqual(len(predicted_passes), len(old_passes), + 'We are loosing passes') + self.assertLessEqual(t2 - t1, t1 - t0, 'Performance is degraded') + + def test_accurate_predictor_find_more_or_equal_passes_amount(self): + self.assertEqualOrGreaterPassesAmount( + Location('bad-case-1', 11.937501570612568, + -55.35189435098657, 1780.674044538666) + ) + self.assertEqualOrGreaterPassesAmount(tortu1) + self.assertEqualOrGreaterPassesAmount(svalbard) + self.assertEqualOrGreaterPassesAmount( + Location('bad-case-2', -11.011509137116818, + 123.29554733688798, 1451.5695915302097) + ) + self.assertEqualOrGreaterPassesAmount( + Location('bad-case-3', 10.20803236163988, + 138.01236517021056, 4967.661890730469) + ) + self.assertEqualOrGreaterPassesAmount( + Location('less passes', -82.41515032683046, + -33.712555446065664, 4417.427841452149) + ) + + def test_predicted_passes_are_equal_between_executions(self): + location = Location('bad-case-1', 11.937501570612568, + -55.35189435098657, 1780.674044538666) + first_set = list( + self.predictor.passes_over(location, self.start, self.end)) + second_set = list( + self.predictor.passes_over(location, self.start + timedelta(seconds=3), self.end) + ) + + self.assertEqual(first_set, second_set) + + def test_predicted_passes_have_elevation_positive_and_visible_on_date(self): + end = self.start + timedelta(days=60) + for pass_ in self.predictor.passes_over(svalbard, self.start, end): + self.assertGreater(pass_.max_elevation_deg, 0) + position = self.old_predictor.get_position(pass_.max_elevation_date) + svalbard.is_visible(position) + self.assertGreaterEqual(pass_.off_nadir_deg, -90) + self.assertLessEqual(pass_.off_nadir_deg, 90) + + def test_predicted_passes_off_nadir_angle_works(self): + start = datetime(2017, 3, 6, 13, 30) + end = start + timedelta(hours=1) + location = Location('bad-case-1', 11.937501570612568, + -55.35189435098657, 1780.674044538666) + + pass_ = self.predictor.get_next_pass(location, when_utc=start, limit_date=end) + self.assertGreaterEqual(0, pass_.off_nadir_deg) + + @given(start=datetimes(min_year=2017, max_year=2020, timezones=[]), + location=tuples( + floats(min_value=-90, max_value=90), + floats(min_value=0, max_value=180), + floats(min_value=-200, max_value=9000) + )) + @settings(max_examples=1000000, timeout=300) + @example(start=datetime(2017, 1, 26, 11, 51, 51), + location=(-37.69358328273305, 153.96875, 0.0)) + def test_pass_is_always_returned(self, start, location): + location = Location('bad-case-1', *location) + pass_ = self.predictor.get_next_pass(location, start) + self.assertGreater(pass_.max_elevation_deg, 0) + + def test_aos_deg_can_be_used_in_get_next_pass(self): + start = datetime(2017, 3, 6, 13, 30) + end = start + timedelta(hours=1) + location = Location('bad-case-1', 11.937501570612568, + -55.35189435098657, 1780.674044538666) + complete_pass = self.predictor.get_next_pass(location, when_utc=start, + limit_date=end) + + pass_with_aos = self.predictor.get_next_pass(location, when_utc=start, + limit_date=end, + aos_at_dg=5) + + self.assertGreater(pass_with_aos.aos, complete_pass.aos) + self.assertLess(pass_with_aos.aos, complete_pass.max_elevation_date) + self.assertAlmostEqual(pass_with_aos.max_elevation_date, + complete_pass.max_elevation_date, + delta=timedelta(seconds=1)) + + self.assertGreater(pass_with_aos.los, complete_pass.max_elevation_date) + self.assertLess(pass_with_aos.los, complete_pass.los) + + position = self.old_predictor.get_position(pass_with_aos.aos) + _, elev = location.get_azimuth_elev_deg(position) + + self.assertAlmostEqual(elev, 5, delta=0.1) + + position = self.old_predictor.get_position(pass_with_aos.los) + _, elev = location.get_azimuth_elev_deg(position) + + self.assertAlmostEqual(elev, 5, delta=0.1) + + def test_predicted_passes_whit_aos(self): + end = self.start + timedelta(days=60) + for pass_ in self.predictor.passes_over(svalbard, self.start, end, aos_at_dg=5): + self.assertGreater(pass_.max_elevation_deg, 5) + position = self.old_predictor.get_position(pass_.aos) + _, elev = svalbard.get_azimuth_elev_deg(position) + self.assertAlmostEqual(elev, 5, delta=0.1) + + +class AccurateVsGpredictTests(TestCase): + + def setUp(self): + # Source + self.db = MemoryTLESource() + self.db.add_tle(BUGSAT_SATE_ID, BUGSAT1_TLE_LINES, datetime.now()) + # Predictor + self.predictor = HighAccuracyTLEPredictor(BUGSAT_SATE_ID, self.db) + + def test_get_next_pass_with_gpredict_data(self): + GPREDICT_DATA = """ + ------------------------------------------------------------------------------------------------- + AOS TCA LOS Duration Max El AOS Az LOS Az + ------------------------------------------------------------------------------------------------- + 2014/10/23 01:27:09 2014/10/23 01:33:03 2014/10/23 01:38:57 00:11:47 25.85 40.28 177.59 + 2014/10/23 03:02:44 2014/10/23 03:08:31 2014/10/23 03:14:17 00:11:32 20.55 341.35 209.65 + 2014/10/23 14:48:23 2014/10/23 14:54:39 2014/10/23 15:00:55 00:12:31 75.31 166.30 350.27 + 2014/10/23 16:25:19 2014/10/23 16:29:32 2014/10/23 16:33:46 00:08:27 7.14 200.60 287.00 + 2014/10/24 01:35:34 2014/10/24 01:41:37 2014/10/24 01:47:39 00:12:05 32.20 34.97 180.38 + 2014/10/24 03:11:40 2014/10/24 03:17:11 2014/10/24 03:22:42 00:11:02 16.30 335.44 213.21 + 2014/10/24 14:57:00 2014/10/24 15:03:16 2014/10/24 15:09:32 00:12:32 84.30 169.06 345.11 + 2014/10/24 16:34:18 2014/10/24 16:38:02 2014/10/24 16:41:45 00:07:27 5.09 205.18 279.57 + 2014/10/25 01:44:01 2014/10/25 01:50:11 2014/10/25 01:56:20 00:12:19 40.61 29.75 183.12 + 2014/10/25 03:20:39 2014/10/25 03:25:51 2014/10/25 03:31:04 00:10:25 12.78 329.21 217.10""" # NOQA + + for line in GPREDICT_DATA.splitlines()[4:]: + line_parts = line.split() + aos = datetime.strptime(" ".join(line_parts[:2]), '%Y/%m/%d %H:%M:%S') + max_elevation_date = datetime.strptime(" ".join(line_parts[2:4]), + '%Y/%m/%d %H:%M:%S') + los = datetime.strptime(" ".join(line_parts[4:6]), '%Y/%m/%d %H:%M:%S') + duration = datetime.strptime(line_parts[6], '%H:%M:%S') + duration_s = timedelta( + minutes=duration.minute, seconds=duration.second).total_seconds() + max_elev_deg = float(line_parts[7]) + + try: + date = pass_.los # NOQA + except UnboundLocalError: + date = datetime.strptime( + "2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + + pass_ = self.predictor.get_next_pass(tortu1, date) + self.assertAlmostEqual(pass_.aos, aos, delta=ONE_SECOND) + self.assertAlmostEqual(pass_.los, los, delta=ONE_SECOND) + self.assertAlmostEqual(pass_.max_elevation_date, max_elevation_date, delta=ONE_SECOND) + self.assertAlmostEqual(pass_.duration_s, duration_s, delta=1) + self.assertAlmostEqual(pass_.max_elevation_deg, max_elev_deg, delta=0.05) + + +class AccuratePredictorCalculationErrorTests(TestCase): + """Check that we can learn from calculation errors and provide patches for corner cases""" + + def setUp(self): + # Source + self.db = MemoryTLESource() + self.db.add_tle(BUGSAT_SATE_ID, BUGSAT1_TLE_LINES, datetime.now()) + # Predictor + self.predictor = HighAccuracyTLEPredictor(BUGSAT_SATE_ID, self.db) + self.is_ascending_mock = self._patch( + 'orbit_predictor.accuratepredictor.LocationPredictor.is_ascending') + self.start = datetime(2017, 3, 6, 7, 51) + logassert.setup(self, 'orbit_predictor.accuratepredictor') + + def _patch(self, *args, **kwargs): + patcher = mock.patch(*args, **kwargs) + self.addCleanup(patcher.stop) + return patcher.start() + + def test_ascending_failure(self): + self.is_ascending_mock.return_value = False + with self.assertRaises(PropagationError): + self.predictor.get_next_pass(svalbard, self.start) + + self.assertLoggedError(str(svalbard), str(self.start), *BUGSAT1_TLE_LINES) + + def test_descending_failure(self): + self.is_ascending_mock.return_value = True + with self.assertRaises(PropagationError): + self.predictor.get_next_pass(svalbard, self.start) + + self.assertLoggedError(str(svalbard), str(self.start), *BUGSAT1_TLE_LINES) diff --git a/tests/test_azimuth_elevation.py b/tests/test_azimuth_elevation.py new file mode 100644 index 0000000..5d09bd5 --- /dev/null +++ b/tests/test_azimuth_elevation.py @@ -0,0 +1,104 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import unittest +from datetime import datetime, timedelta + +from orbit_predictor.utils import sun_azimuth_elevation + + +class AzimuthElevationTests(unittest.TestCase): + """Ensures calculations are ok, comparing against: + http://www.esrl.noaa.gov/gmd/grad/solcalc/ + all dates are UTC + """ + def setUp(self): + self.coords = (0, 0) + self.base_date = datetime(2016, 9, 8) + + def test_return_value(self): + obj = sun_azimuth_elevation(*self.coords, when=self.base_date) + self.assertTrue(hasattr(obj, 'azimuth')) + self.assertTrue(hasattr(obj, 'elevation')) + + def test_8(self): + date = self.base_date + timedelta(hours=8) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 83.64, delta=0.5) + self.assertAlmostEqual(elevation, 30.67, delta=0.5) + + def test_9(self): + date = self.base_date + timedelta(hours=9) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 82.21, delta=0.5) + self.assertAlmostEqual(elevation, 45.36, delta=0.5) + + def test_10(self): + date = self.base_date + timedelta(hours=10) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 79, delta=0.5) + self.assertAlmostEqual(elevation, 60.16, delta=0.5) + + def test_11(self): + date = self.base_date + timedelta(hours=11) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 69.05, delta=0.5) + self.assertAlmostEqual(elevation, 74.64, delta=0.5) + + def test_12(self): + date = self.base_date + timedelta(hours=12) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 353.54, delta=0.5) + self.assertAlmostEqual(elevation, 84.55, delta=0.5) + + def test_13(self): + date = self.base_date + timedelta(hours=13) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 289.36, delta=0.5) + self.assertAlmostEqual(elevation, 73.5, delta=0.5) + + def test_14(self): + date = self.base_date + timedelta(hours=14) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 280.49, delta=0.5) + self.assertAlmostEqual(elevation, 58.96, delta=0.5) + + def test_15(self): + date = self.base_date + timedelta(hours=15) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 277.49, delta=0.5) + self.assertAlmostEqual(elevation, 44.14, delta=0.5) + + def test_16(self): + date = self.base_date + timedelta(hours=16) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 276.14, delta=0.5) + self.assertAlmostEqual(elevation, 29.26, delta=0.5) + + def test_19(self): + date = self.base_date + timedelta(hours=19) + azimuth, elevation = sun_azimuth_elevation(*self.coords, when=date) + self.assertAlmostEqual(azimuth, 275.51, delta=0.5) + self.assertAlmostEqual(elevation, -15.55, delta=0.5) + + def test_default_when(self): + sun_azimuth_elevation(0, 0) diff --git a/tests/test_coordinate_systems.py b/tests/test_coordinate_systems.py new file mode 100644 index 0000000..448d850 --- /dev/null +++ b/tests/test_coordinate_systems.py @@ -0,0 +1,109 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import unittest + +from orbit_predictor import coordinate_systems as coords + +buenos_aires_llh = (-34.628284, -58.436045, 0.0) +lujan_llh = (-34.561207, -59.114542, 0.0) +buenos_aires_lujan_distance = 64.0 +north_pole_llh = (84.0, 0.0, 0.0) +sixty_four_km_away_from_north_pole_llh = (83.53837, -3.230163, 0.0) + + +class TestLLHToECEF(unittest.TestCase): + + def test_buenos_aires(self): + ecef = coords.llh_to_ecef(buenos_aires_llh[0], buenos_aires_llh[1], .025) + self.assertAlmostEqual(ecef[0], 2750.19, delta=.2) + self.assertAlmostEqual(ecef[1], -4476.679, delta=.2) + self.assertAlmostEqual(ecef[2], -3604.011, delta=.2) + + def test_buenos_aires_high(self): + ecef = coords.llh_to_ecef(buenos_aires_llh[0], buenos_aires_llh[1], 400.0) + self.assertAlmostEqual(ecef[0], 2922.48, delta=.2) + self.assertAlmostEqual(ecef[1], -4757.126, delta=.2) + self.assertAlmostEqual(ecef[2], -3831.311, delta=.2) + + def test_lisbon(self): + ecef = coords.llh_to_ecef(38.716666, 9.16666, .002) + self.assertAlmostEqual(ecef[0], 4919.4248, delta=.2) + self.assertAlmostEqual(ecef[1], 793.8355, delta=.2) + self.assertAlmostEqual(ecef[2], 3967.8253, delta=.2) + + def test_lisbon_high(self): + ecef = coords.llh_to_ecef(38.716666, 9.16666, 450) + self.assertAlmostEqual(ecef[0], 5266.051, delta=.2) + self.assertAlmostEqual(ecef[1], 849.7698, delta=.2) + self.assertAlmostEqual(ecef[2], 4249.2854, delta=.2) + + def test_north_pole(self): + ecef = coords.llh_to_ecef(89.99, -70, .025) + self.assertAlmostEqual(ecef[0], .382, delta=.2) + self.assertAlmostEqual(ecef[1], -1.0495, delta=.2) + self.assertAlmostEqual(ecef[2], 6356.7772, delta=.2) + + def test_south_pole(self): + ecef = coords.llh_to_ecef(-89.99, -70, .025) + self.assertAlmostEqual(ecef[0], .382, delta=.2) + self.assertAlmostEqual(ecef[1], -1.0495, delta=.2) + self.assertAlmostEqual(ecef[2], -6356.7772, delta=.2) + + +class TestECEFToLLH(unittest.TestCase): + + def test_buenos_aires(self): + llh = coords.ecef_to_llh((2750.19, -4476.679, -3604.011)) + self.assertAlmostEqual(llh[0], buenos_aires_llh[0], delta=.2) + self.assertAlmostEqual(llh[1], buenos_aires_llh[1], delta=.2) + self.assertAlmostEqual(llh[2], 0, delta=.2) + + def test_buenos_aires_high(self): + llh = coords.ecef_to_llh((2922.48, -4757.126, -3831.311)) + self.assertAlmostEqual(llh[0], buenos_aires_llh[0], delta=.2) + self.assertAlmostEqual(llh[1], buenos_aires_llh[1], delta=.2) + self.assertAlmostEqual(llh[2], 400, delta=.2) + + def test_lisbon(self): + llh = coords.ecef_to_llh((4919.4248, 793.8355, 3967.8253)) + self.assertAlmostEqual(llh[0], 38.716666, delta=.2) + self.assertAlmostEqual(llh[1], 9.16666, delta=.2) + self.assertAlmostEqual(llh[2], .002, delta=.2) + + def test_lisbon_high(self): + llh = coords.ecef_to_llh((5266.051, 849.7698, 4249.2854)) + self.assertAlmostEqual(llh[0], 38.716666, delta=.2) + self.assertAlmostEqual(llh[1], 9.16666, delta=.2) + self.assertAlmostEqual(llh[2], 450, delta=.2) + + def test_north_pole(self): + llh = coords.ecef_to_llh((.382, -1.0495, 6356.7772)) + self.assertAlmostEqual(llh[0], 89.99, delta=.2) + self.assertAlmostEqual(llh[1], -70, delta=.2) + self.assertAlmostEqual(llh[2], .025, delta=.2) + + def test_south_pole(self): + llh = coords.ecef_to_llh((.382, -1.0495, -6356.7772)) + self.assertAlmostEqual(llh[0], -89.99, delta=.2) + self.assertAlmostEqual(llh[1], -70, delta=.2) + self.assertAlmostEqual(llh[2], .025, delta=.2) diff --git a/tests/test_locations.py b/tests/test_locations.py new file mode 100644 index 0000000..652c966 --- /dev/null +++ b/tests/test_locations.py @@ -0,0 +1,109 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import datetime +import unittest +from math import degrees + +from orbit_predictor.locations import Location, tortu1 +from orbit_predictor.predictors import TLEPredictor +from orbit_predictor.sources import MemoryTLESource + +SATE_ID = 'BUGSAT-1' +BUGSAT1_TLE_LINES = ("1 40014U 14033E 14294.41438078 .00003468 00000-0 34565-3 0 3930", + "2 40014 97.9781 190.6418 0032692 299.0467 60.7524 14.91878099 18425") + + +class LocationTestCase(unittest.TestCase): + + def setUp(self): + # Source + self.db = MemoryTLESource() + self.db.add_tle(SATE_ID, BUGSAT1_TLE_LINES, datetime.datetime.now()) + # Predictor + self.predictor = TLEPredictor(SATE_ID, self.db) + date = datetime.datetime.strptime("2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + self.next_pass = self.predictor.get_next_pass(tortu1, when_utc=date) + + def test_compare_eq(self): + l1 = Location(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location1") + l2 = Location(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location1") + + self.assertEqual(l1, l2) + self.assertEqual(l2, l1) + + def test_compare_no_eq(self): + l1 = Location(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location_other") + l2 = Location(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location1") + + self.assertNotEqual(l1, l2) + self.assertNotEqual(l2, l1) + + def test_compare_eq_subclass(self): + + class SubLocation(Location): + pass + + l1 = Location(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location1") + l2 = SubLocation(latitude_deg=1, longitude_deg=2, elevation_m=3, name="location1") + + self.assertEqual(l1, l2) + self.assertEqual(l2, l1) + + def test_get_azimuth_elev(self): + date = datetime.datetime.strptime("2014-10-21 22:47:29.147740", '%Y-%m-%d %H:%M:%S.%f') + azimuth, elevation = tortu1.get_azimuth_elev(self.predictor.get_position(date)) + + self.assertAlmostEqual(degrees(azimuth), 245.1, delta=0.1) + self.assertAlmostEqual(degrees(elevation), -53.6, delta=0.1) + + def test_get_azimuth_elev_deg(self): + date = datetime.datetime.strptime("2014-10-21 22:47:29.147740", '%Y-%m-%d %H:%M:%S.%f') + azimuth, elevation = tortu1.get_azimuth_elev_deg(self.predictor.get_position(date)) + + self.assertAlmostEqual(azimuth, 245.1, delta=0.1) + self.assertAlmostEqual(elevation, -53.6, delta=0.1) + + def test_is_visible(self): + position = self.predictor.get_position(self.next_pass.aos) + self.assertTrue(tortu1.is_visible(position)) + + def test_no_visible(self): + position = self.predictor.get_position(self.next_pass.los + datetime.timedelta(minutes=10)) + self.assertFalse(tortu1.is_visible(position)) + + def test_is_visible_with_deg(self): + position = self.predictor.get_position(self.next_pass.aos + datetime.timedelta(minutes=4)) + # 21 deg + self.assertTrue(tortu1.is_visible(position, elevation=4)) + + def test_no_visible_with_deg(self): + position = self.predictor.get_position(self.next_pass.aos + datetime.timedelta(minutes=4)) + # 21 deg + self.assertFalse(tortu1.is_visible(position, elevation=30)) + + def test_doppler_factor(self): + date = datetime.datetime.strptime("2014-10-21 23:06:11.132438", '%Y-%m-%d %H:%M:%S.%f') + position = self.predictor.get_position(date) + doppler_factor = tortu1.doppler_factor(position) + + self.assertAlmostEqual((2 - doppler_factor)*437.445e6, 437.445796e6, delta=100) diff --git a/tests/test_position.py b/tests/test_position.py new file mode 100644 index 0000000..9dd35fd --- /dev/null +++ b/tests/test_position.py @@ -0,0 +1,39 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +import unittest + +from orbit_predictor import coordinate_systems +from orbit_predictor.predictors import Position + + +class PositionTestCase(unittest.TestCase): + def test_position_convertion(self): + position_ecef = (2750.19, -4476.679, -3604.011) + + p = Position(None, position_ecef, None, None) + ecef_est = coordinate_systems.llh_to_ecef(*p.position_llh) + + self.assertAlmostEqual(position_ecef[0], ecef_est[0], 3) + self.assertAlmostEqual(position_ecef[1], ecef_est[1], 3) + self.assertAlmostEqual(position_ecef[2], ecef_est[2], 3) diff --git a/tests/test_predicted_pass.py b/tests/test_predicted_pass.py new file mode 100644 index 0000000..2432193 --- /dev/null +++ b/tests/test_predicted_pass.py @@ -0,0 +1,45 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import unittest +from datetime import datetime, timedelta + +from orbit_predictor.locations import tortu1 +from orbit_predictor.predictors import PredictedPass + + +class PredictedPassTests(unittest.TestCase): + + def test_midpoint(self): + aos = datetime.utcnow() + max_elevation_date = aos + timedelta(minutes=5) + los = aos + timedelta(minutes=10) + pass_ = PredictedPass(sate_id=1, + location=tortu1, + aos=aos, + los=los, + duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=None, + max_elevation_deg=10) + + self.assertEqual(pass_.midpoint, aos + timedelta(minutes=5)) diff --git a/tests/test_sources.py b/tests/test_sources.py new file mode 100644 index 0000000..4e3fa2d --- /dev/null +++ b/tests/test_sources.py @@ -0,0 +1,184 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import datetime +import os +import shutil +import tempfile +import unittest + +from orbit_predictor import sources +from orbit_predictor.accuratepredictor import HighAccuracyTLEPredictor +from orbit_predictor.predictors import TLEPredictor + +try: + from unittest.mock import Mock, patch +except ImportError: + from mock import Mock, patch # Python2 + +try: + from urllib import parse as urlparse +except ImportError: + import urlparse # Python2 + + +SATE_ID = "AAUSAT-II" +SAMPLE_TLE = ("1 32788U 08021F 15227.82608814 .00001480 00000-0 15110-3 0 9997", + "2 32788 97.6474 275.2739 0011863 204.9398 155.1249 14.92031413395491") + +SAMPLE_TLE2 = ("1 32791U 08021J 15228.17298173 .00001340 00000-0 14806-3 0 9999", + "2 32791 97.6462 271.6584 0012961 215.4867 144.5490 14.88966377395242") + + +class TestMemoryTLESource(unittest.TestCase): + def setUp(self): + self.db = sources.MemoryTLESource() + + def test_add_tle(self): + self.db.add_tle(SATE_ID, SAMPLE_TLE, datetime.datetime.now()) + tle = self.db._get_tle(SATE_ID, datetime.datetime.now()) + self.assertEqual(tle, SAMPLE_TLE) + + def test_add_tle_twice(self): + self.db.add_tle(SATE_ID, SAMPLE_TLE, datetime.datetime.now()) + self.db.add_tle(SATE_ID, SAMPLE_TLE2, datetime.datetime.now()) + tle = self.db._get_tle(SATE_ID, datetime.datetime.now()) + self.assertEqual(tle, SAMPLE_TLE2) + + def test_add_tle_two_id(self): + self.db.add_tle(SATE_ID, SAMPLE_TLE, datetime.datetime.now()) + self.db.add_tle("fake_id", SAMPLE_TLE2, datetime.datetime.now()) + tle = self.db._get_tle(SATE_ID, datetime.datetime.now()) + self.assertEqual(tle, SAMPLE_TLE) + + def test_empty(self): + with self.assertRaises(LookupError): + self.db._get_tle(SATE_ID, datetime.datetime.now()) + + # this methods are from TLESource() + def test_get(self): + date = datetime.datetime.now() + self.db.add_tle(SATE_ID, SAMPLE_TLE, date) + tle = self.db.get_tle(SATE_ID, date) + self.assertEqual(tle.lines, SAMPLE_TLE) + self.assertEqual(tle.sate_id, SATE_ID) + self.assertEqual(tle.date, date) + + def test_get_predictor(self): + predictor = self.db.get_predictor(SATE_ID) + + self.assertIsInstance(predictor, TLEPredictor) + self.assertEqual(predictor.sate_id, SATE_ID) + self.assertEqual(predictor.source, self.db) + + def test_get_predictor_precise(self): + predictor = self.db.get_predictor(SATE_ID, precise=True) + self.assertIsInstance(predictor, HighAccuracyTLEPredictor) + self.assertEqual(predictor.sate_id, SATE_ID) + self.assertEqual(predictor.source, self.db) + + +class TestEtcTLESource(unittest.TestCase): + def setUp(self): + self.dirname = tempfile.mkdtemp() + self.filename = os.path.join(self.dirname, "tle_file") + + with open(self.filename, "w") as fd: + fd.write(SATE_ID + "\n") + for l in SAMPLE_TLE: + fd.write(l + "\n") + + def test_add_tle(self): + db = sources.EtcTLESource(self.filename) + + db.add_tle(SATE_ID, SAMPLE_TLE2, datetime.datetime.now()) + tle = db._get_tle(SATE_ID, datetime.datetime.now()) + self.assertEqual(tle, SAMPLE_TLE2) + + def test_read_tle(self): + db = sources.EtcTLESource(self.filename) + + tle = db._get_tle(SATE_ID, datetime.datetime.now()) + self.assertEqual(tle, SAMPLE_TLE) + + def test_wrong_sate(self): + db = sources.EtcTLESource(self.filename) + + with self.assertRaises(LookupError): + db._get_tle("fake_id", datetime.datetime.now()) + + def tearDown(self): + shutil.rmtree(self.dirname) + + +class TestWSTLESource(unittest.TestCase): + def setUp(self): + self.mock_json = { + "date": "2015-01-15T08:56:33Z", + "lines": [ + "1 40014U 14033E 15014.37260739 .00003376 00000-0 33051-3 0 6461", + "2 40014 97.9772 276.7067 0034296 17.9121 342.3159 14.92570778 31092" + ] + } + + self.expected_lines = ( + "1 40014U 14033E 15014.37260739 .00003376 00000-0 33051-3 0 6461", + "2 40014 97.9772 276.7067 0034296 17.9121 342.3159 14.92570778 31092") + + self.headers = {'user-agent': 'orbit-predictor', 'Accept': 'application/json'} + + @patch("requests.get") + def test_get_tle(self, mocked_requests): + expected_url = urlparse.urlparse( + "http://test.none/api/tle/closest/?date=2015-01-01&satellite_number=40014U") + expected_qs = urlparse.parse_qs(expected_url.query) + + mocked_response = Mock() + mocked_response.ok = True + mocked_response.json.return_value = self.mock_json + mocked_requests.return_value = mocked_response + + source = sources.WSTLESource(url="http://test.none/") + tle = source._get_tle('40014U', datetime.datetime(2015, 1, 1)) + + call_args = mocked_requests.call_args + url = urlparse.urlparse(call_args[0][0]) + url_qs = urlparse.parse_qs(url.query) + + self.assertEqual(url.path, expected_url.path) + self.assertEqual(url_qs, expected_qs) + self.assertEqual(call_args[1], {'headers': self.headers}) + self.assertEqual(tle, self.expected_lines) + + @patch("requests.get") + def test_get_last_update(self, mocked_requests): + url = "http://test.none/api/tle/last/?satellite_number=40014U" + mocked_response = Mock() + mocked_response.ok = True + mocked_response.json.return_value = self.mock_json + mocked_requests.return_value = mocked_response + + source = sources.WSTLESource(url="http://test.none/") + tle = source.get_last_update('40014U') + + mocked_requests.assert_called_with(url, headers=self.headers) + self.assertEqual(tle, self.expected_lines) diff --git a/tests/test_sun_elevation.py b/tests/test_sun_elevation.py new file mode 100644 index 0000000..6c07289 --- /dev/null +++ b/tests/test_sun_elevation.py @@ -0,0 +1,215 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import datetime +import time +import unittest + +from orbit_predictor.utils import sun_azimuth_elevation + +buenos_aires_llh = (-34.628284, -58.436045) +buenos_aires_timezone = -3 +north_pole_llh = (84.0, 0.0) +greenwich_timezone = 0 +tokio_llh = (35.6833333, 139.7666666) +tokio_timezone = 9 + + +class ElevationTestCase(unittest.TestCase): + + def assertElevation(self, date, expected): + lat, lon = self.location + _, actual = sun_azimuth_elevation(lat, lon, date) + self.assertAlmostEqual(actual, expected, delta=4) + + +class TestOnBuenosAires(ElevationTestCase): + + def setUp(self): + super(self.__class__, self).setUp() + + self.january_10th = Date.from_utc('10/1/2008', buenos_aires_timezone) + self.march_10th = Date.from_utc('10/3/2008', buenos_aires_timezone) + self.june_10th = Date.from_utc('10/6/2008', buenos_aires_timezone) + self.september_10th = Date.from_utc('10/9/2008', buenos_aires_timezone) + + self.location = buenos_aires_llh + + def test_q1(self): + self.assertElevation(self.january_10th.noon(), 71.58688) + self.assertElevation(self.january_10th.dawn(), 12.10156) + self.assertElevation(self.january_10th.dust(), 12.10156) + + def test_q2(self): + self.assertElevation(self.march_10th.noon(), 55.78998) + self.assertElevation(self.march_10th.dawn(), 1.39351) + self.assertElevation(self.march_10th.dust(), 2.87432) + + def test_q3(self): + self.assertElevation(self.june_10th.noon(), 30.93056) + self.assertElevation(self.june_10th.dawn(), -11.55022) + self.assertElevation(self.june_10th.dust(), -14.17221) + + def test_q4(self): + self.assertElevation(self.september_10th.noon(), 48.96952) + self.assertElevation(self.september_10th.dawn(), -0.75196) + self.assertElevation(self.september_10th.dust(), -4.5362) + + +class TestOnNorthPole(ElevationTestCase): + + def setUp(self): + super(self.__class__, self).setUp() + + self.january_10th = Date.from_utc('10/1/2008', greenwich_timezone) + self.march_10th = Date.from_utc('10/3/2008', greenwich_timezone) + self.june_10th = Date.from_utc('10/6/2008', greenwich_timezone) + self.september_10th = Date.from_utc('10/9/2008', greenwich_timezone) + + self.location = north_pole_llh + + def test_q1(self): + self.assertElevation(self.january_10th.noon(), -16.00948) + self.assertElevation(self.january_10th.dawn(), -20.55059) + self.assertElevation(self.january_10th.dust(), -23.20414) + + def test_q2(self): + self.assertElevation(self.march_10th.noon(), 2.14842) + self.assertElevation(self.march_10th.dawn(), -2.6164) + self.assertElevation(self.march_10th.dust(), -5.00546) + + def test_q3(self): + self.assertElevation(self.june_10th.noon(), 29.06153) + self.assertElevation(self.june_10th.dawn(), 24.4841) + self.assertElevation(self.june_10th.dust(), 21.39807) + + def test_q4(self): + self.assertElevation(self.september_10th.noon(), 10.68454) + self.assertElevation(self.september_10th.dawn(), 6.36817) + self.assertElevation(self.september_10th.dust(), 2.91842) + + +class TestOnTokio(ElevationTestCase): + + def setUp(self): + super(self.__class__, self).setUp() + + self.january_10th = Date.from_utc('10/1/2008', tokio_timezone) + self.march_10th = Date.from_utc('10/3/2008', tokio_timezone) + self.june_10th = Date.from_utc('10/6/2008', tokio_timezone) + self.september_10th = Date.from_utc('10/9/2008', tokio_timezone) + + self.location = tokio_llh + + def test_q1(self): + self.assertElevation(self.january_10th.noon(), 32.18769) + self.assertElevation(self.january_10th.dawn(), 0.75589) + self.assertElevation(self.january_10th.dust(), -26.7901) + + def test_q2(self): + self.assertElevation(self.march_10th.noon(), 50.27033) + self.assertElevation(self.march_10th.dawn(), 11.41627) + self.assertElevation(self.march_10th.dust(), -16.22019) + + def test_q3(self): + self.assertElevation(self.june_10th.noon(), 76.65314) + self.assertElevation(self.june_10th.dawn(), 28.86539) + self.assertElevation(self.june_10th.dust(), -1.48132) + + def test_q4(self): + self.assertElevation(self.september_10th.noon(), 58.7265) + self.assertElevation(self.september_10th.dawn(), 19.47538) + self.assertElevation(self.september_10th.dust(), -13.66267) + + +class Date(object): + def __init__(self, dt, tz): + self.date = dt - datetime.timedelta(hours=tz) + + def noon(self): + return self.date + datetime.timedelta(hours=12) + + def dawn(self): + return self.date + datetime.timedelta(hours=7) + + def dust(self): + return self.date + datetime.timedelta(hours=19) + + @classmethod + def from_utc(cls, time_string, tz): + t = time.strptime(time_string, '%d/%m/%Y') + return cls(datetime.datetime(t.tm_year, t.tm_mon, t.tm_mday, tzinfo=UTC), tz) + + +# NOTE: These classes seem necessary to allow the dates to be introduced as UTC and then +# localized so that when input to get_sun_azimuth_elevation they are interpreted +# correctly as the original UTC... +# It's horrible, and if you know of a better way to do this please let me know. +class UTCTimezone(datetime.tzinfo): + ZERO = datetime.timedelta(0) + + def utcoffset(self, dt): + return UTCTimezone.ZERO + + def dst(self, dt): + return UTCTimezone.ZERO + + def tzname(self, dt): + return 'UTC' + + +UTC = UTCTimezone() + + +class LocalTimezone(datetime.tzinfo): + STDOFFSET = datetime.timedelta(seconds=-time.timezone) + if time.daylight: + DSTOFFSET = datetime.timedelta(seconds=-time.altzone) + else: + DSTOFFSET = STDOFFSET + DSTDIFF = DSTOFFSET - STDOFFSET + + def utcoffset(self, dt): + if self._isdst(dt): + return LocalTimezone.DSTOFFSET + else: + return LocalTimezone.STDOFFSET + + def dst(self, dt): + if self._isdst(dt): + return LocalTimezone.DSTDIFF + else: + return UTCTimezone.ZERO + + def tzname(self, dt): + return time.tzname[self._isdst(dt)] + + def _isdst(self, dt): + tt = (dt.year, dt.month, dt.day, + dt.hour, dt.minute, dt.second, + dt.weekday(), 0, 0) + stamp = time.mktime(tt) + tt = time.localtime(stamp) + return tt.tm_isdst > 0 + + +Local = LocalTimezone() diff --git a/tests/test_tle_predictor.py b/tests/test_tle_predictor.py new file mode 100644 index 0000000..3364c13 --- /dev/null +++ b/tests/test_tle_predictor.py @@ -0,0 +1,279 @@ +# MIT License +# +# Copyright (c) 2017 Satellogic SA +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import datetime +import unittest + +from orbit_predictor.coordinate_systems import llh_to_ecef +from orbit_predictor.locations import Location, tortu1 +from orbit_predictor.predictors import ( + NotReachable, + Position, + PredictedPass, + TLEPredictor +) +from orbit_predictor.sources import MemoryTLESource + +try: + from unittest.mock import patch +except ImportError: + from mock import patch # Python2 + + +SATE_ID = 'BUGSAT-1' +BUGSAT1_TLE_LINES = ( + "1 40014U 14033E 14294.41438078 .00003468 00000-0 34565-3 0 3930", + "2 40014 97.9781 190.6418 0032692 299.0467 60.7524 14.91878099 18425") + + +class TLEPredictorTestCase(unittest.TestCase): + + @classmethod + def setUpClass(cls): + # Source + cls.db = MemoryTLESource() + cls.db.add_tle(SATE_ID, BUGSAT1_TLE_LINES, datetime.datetime.now()) + # Predictor + cls.predictor = TLEPredictor(SATE_ID, cls.db) + + def test_predicted_pass_eq(self): + aos = datetime.datetime.utcnow() + max_elevation_date = datetime.datetime.utcnow() + datetime.timedelta(minutes=5) + los = datetime.datetime.utcnow() + datetime.timedelta(minutes=10) + max_elevation_position = Position( + when_utc=max_elevation_date, + position_ecef=(1, 1, 1), + velocity_ecef=(1, 1, 1), + error_estimate=0) + + p1 = PredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=10) + p2 = PredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=10) + + self.assertEqual(p1, p2) + self.assertEqual(p2, p1) + + def test_predicted_pass_no_eq(self): + aos = datetime.datetime.utcnow() + max_elevation_date = datetime.datetime.utcnow() + datetime.timedelta(minutes=5) + los = datetime.datetime.utcnow() + datetime.timedelta(minutes=10) + max_elevation_position = Position( + when_utc=max_elevation_date, + position_ecef=(1, 1, 1), + velocity_ecef=(1, 1, 1), + error_estimate=0) + + p1 = PredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=10) + p2 = PredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=50) + + self.assertNotEqual(p1, p2) + self.assertNotEqual(p2, p1) + + def test_predicted_pass_eq_subclass(self): + + class SubPredictedPass(PredictedPass): + pass + + aos = datetime.datetime.utcnow() + max_elevation_date = datetime.datetime.utcnow() + datetime.timedelta(minutes=5) + los = datetime.datetime.utcnow() + datetime.timedelta(minutes=10) + max_elevation_position = Position( + when_utc=max_elevation_date, + position_ecef=(1, 1, 1), + velocity_ecef=(1, 1, 1), + error_estimate=0) + + p1 = PredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=10) + p2 = SubPredictedPass( + sate_id=1, location=tortu1, + aos=aos, los=los, duration_s=600, + max_elevation_date=max_elevation_date, + max_elevation_position=max_elevation_position, + max_elevation_deg=10) + + self.assertEqual(p1, p2) + self.assertEqual(p2, p1) + + def test_get_next_pass_with_gpredict_data(self): + GPREDICT_DATA = """ + ------------------------------------------------------------------------------------------------- + AOS TCA LOS Duration Max El AOS Az LOS Az + ------------------------------------------------------------------------------------------------- + 2014/10/23 01:27:09 2014/10/23 01:33:03 2014/10/23 01:38:57 00:11:47 25.85 40.28 177.59 + 2014/10/23 03:02:44 2014/10/23 03:08:31 2014/10/23 03:14:17 00:11:32 20.55 341.35 209.65 + 2014/10/23 14:48:23 2014/10/23 14:54:39 2014/10/23 15:00:55 00:12:31 75.31 166.30 350.27 + 2014/10/23 16:25:19 2014/10/23 16:29:32 2014/10/23 16:33:46 00:08:27 7.14 200.60 287.00 + 2014/10/24 01:35:34 2014/10/24 01:41:37 2014/10/24 01:47:39 00:12:05 32.20 34.97 180.38 + 2014/10/24 03:11:40 2014/10/24 03:17:11 2014/10/24 03:22:42 00:11:02 16.30 335.44 213.21 + 2014/10/24 14:57:00 2014/10/24 15:03:16 2014/10/24 15:09:32 00:12:32 84.30 169.06 345.11 + 2014/10/24 16:34:18 2014/10/24 16:38:02 2014/10/24 16:41:45 00:07:27 5.09 205.18 279.57 + 2014/10/25 01:44:01 2014/10/25 01:50:11 2014/10/25 01:56:20 00:12:19 40.61 29.75 183.12 + 2014/10/25 03:20:39 2014/10/25 03:25:51 2014/10/25 03:31:04 00:10:25 12.78 329.21 217.10""" # NOQA + + for line in GPREDICT_DATA.splitlines()[4:]: + line_parts = line.split() + aos = datetime.datetime.strptime(" ".join(line_parts[:2]), '%Y/%m/%d %H:%M:%S') + max_elevation_date = datetime.datetime.strptime(" ".join(line_parts[2:4]), + '%Y/%m/%d %H:%M:%S') + los = datetime.datetime.strptime(" ".join(line_parts[4:6]), '%Y/%m/%d %H:%M:%S') + duration = datetime.datetime.strptime(line_parts[6], '%H:%M:%S') + duration_s = datetime.timedelta( + minutes=duration.minute, seconds=duration.second).total_seconds() + max_elev_deg = float(line_parts[7]) + + try: + date = pass_.los # NOQA + except UnboundLocalError: + date = datetime.datetime.strptime( + "2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + pass_ = self.predictor.get_next_pass(tortu1, date) + self.assertAlmostEqual((pass_.aos - aos).total_seconds(), 0, delta=25) + self.assertAlmostEqual((pass_.max_elevation_date - max_elevation_date).total_seconds(), + 0, delta=25) + self.assertAlmostEqual((pass_.los - los).total_seconds(), 0, delta=25) + self.assertAlmostEqual(pass_.max_elevation_deg, max_elev_deg, delta=1) + self.assertAlmostEqual(pass_.duration_s, duration_s, delta=10) + + def test_get_next_pass(self): + date = datetime.datetime.strptime("2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + pass_ = self.predictor.get_next_pass(tortu1, date, max_elevation_gt=15) + for i in range(20): + pass_ = self.predictor.get_next_pass(tortu1, pass_.los, max_elevation_gt=15) + self.assertLess(self.predictor._iterations, 200) + self.assertGreater(self.predictor._iterations, 30) + self.assertGreaterEqual(pass_.max_elevation_deg, 15) + + def test_get_next_pass_with_limit_exception(self): + date = datetime.datetime.strptime("2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + pass_ = self.predictor.get_next_pass(tortu1, date, max_elevation_gt=15) + with self.assertRaises(NotReachable): + self.predictor.get_next_pass(tortu1, date, max_elevation_gt=15, + limit_date=pass_.aos - datetime.timedelta(minutes=1)) + + def test_get_next_pass_with_limit(self): + date = datetime.datetime.strptime("2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + pass_ = self.predictor.get_next_pass(tortu1, date, max_elevation_gt=15) + new_pass = self.predictor.get_next_pass( + tortu1, date, max_elevation_gt=15, + limit_date=pass_.los + datetime.timedelta(seconds=1)) + self.assertEqual(pass_, new_pass) + + def test_get_next_pass_while_passing(self): + date = datetime.datetime.strptime("2014/10/23 01:32:09", '%Y/%m/%d %H:%M:%S') + pass_ = self.predictor.get_next_pass(tortu1, date) + self.assertEqual(pass_.aos, date) + self.assertTrue(date < pass_.los) + + position = self.predictor.get_position(date) + self.assertTrue(tortu1.is_visible(position)) + + def test_grater_than_deg(self): + date = datetime.datetime.strptime("2014/10/23 01:25:09", '%Y/%m/%d %H:%M:%S') + pass5 = self.predictor.get_next_pass(tortu1, date, aos_at_dg=5) + pass10 = self.predictor.get_next_pass(tortu1, date, aos_at_dg=10) + pass15 = self.predictor.get_next_pass(tortu1, date, aos_at_dg=15) + self.assertTrue(pass5.aos < pass10.aos < pass15.aos) + self.assertTrue(pass5.los > pass10.los > pass15.los) + self.assertTrue(pass5.max_elevation_deg == pass10.max_elevation_deg) + + @patch("orbit_predictor.predictors.TLEPredictor._propagate_ecef") + def test_get_position(self, mocked_propagate): + mocked_propagate.return_value = ('foo', 'bar') + when_utc = datetime.datetime.utcnow() + position = self.predictor.get_position(when_utc) + + self.assertIsInstance(position, Position) + self.assertEqual(when_utc, position.when_utc) + self.assertIsNone(position.error_estimate) + self.assertEqual(position.position_ecef, 'foo') + self.assertEqual(position.velocity_ecef, 'bar') + + def test_off_nadir_computable_and_reasonable(self): + date = datetime.datetime.strptime("2014-10-22 20:18:11.921921", '%Y-%m-%d %H:%M:%S.%f') + pass_ = self.predictor.get_next_pass(tortu1, date) + self.assertLessEqual(abs(pass_.off_nadir_deg), 90) + + +class OffNadirAngleTests(unittest.TestCase): + def setUp(self): + self.location = Location("A random location", latitude_deg=0, longitude_deg=0, + elevation_m=0) + + def test_off_nadir_satellite_exactly_over(self): + position_ecef = llh_to_ecef(0, 0, 500 * 1000) # A satellite exactyl over the point + velocity_ecef = (0, 0, 1) + max_elevation_position = Position(None, position_ecef, velocity_ecef, None) + + pass_ = PredictedPass(sate_id=1, location=self.location, + aos=None, los=None, duration_s=None, + max_elevation_position=max_elevation_position, + max_elevation_deg=None) + + self.assertAlmostEqual(pass_.off_nadir_deg, 0, delta=0.02) + + def test_off_nadir_satellite_passing_left_means_positive_sign(self): + position_ecef = llh_to_ecef(0, -10, 500 * 1000) # A satellite exactyl over the point + velocity_ecef = (0, 0, -1) + max_elevation_position = Position(None, position_ecef, velocity_ecef, None) + + pass_ = PredictedPass(sate_id=1, location=self.location, + aos=None, los=None, duration_s=None, + max_elevation_position=max_elevation_position, + max_elevation_deg=None) + + self.assertGreaterEqual(pass_.off_nadir_deg, 0) + + def test_off_nadir_satellite_passing_right_means_negative_sign(self): + position_ecef = llh_to_ecef(0, 10, 500 * 1000) # A satellite exactyl over the point + velocity_ecef = (0, 0, -1) + max_elevation_position = Position(None, position_ecef, velocity_ecef, None) + + pass_ = PredictedPass(sate_id=1, location=self.location, + aos=None, los=None, duration_s=None, + max_elevation_position=max_elevation_position, + max_elevation_deg=None) + + self.assertLessEqual(pass_.off_nadir_deg, 0)