From bced15602cf083ddd50b05f238731a5515b2c32b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Tue, 14 Nov 2017 15:11:56 +0100 Subject: [PATCH 1/6] Move predictors to their own package --- orbit_predictor/predictors.py | 246 ------------------------- orbit_predictor/predictors/__init__.py | 3 + orbit_predictor/predictors/base.py | 135 ++++++++++++++ orbit_predictor/predictors/tle.py | 91 +++++++++ orbit_predictor/utils.py | 8 + 5 files changed, 237 insertions(+), 246 deletions(-) delete mode 100644 orbit_predictor/predictors.py create mode 100644 orbit_predictor/predictors/__init__.py create mode 100644 orbit_predictor/predictors/base.py create mode 100644 orbit_predictor/predictors/tle.py diff --git a/orbit_predictor/predictors.py b/orbit_predictor/predictors.py deleted file mode 100644 index 7e872d9..0000000 --- a/orbit_predictor/predictors.py +++ /dev/null @@ -1,246 +0,0 @@ -# -*- 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/predictors/__init__.py b/orbit_predictor/predictors/__init__.py new file mode 100644 index 0000000..572ba15 --- /dev/null +++ b/orbit_predictor/predictors/__init__.py @@ -0,0 +1,3 @@ +from orbit_predictor.predictors.base import Position, PredictedPass # noqa +from orbit_predictor.exceptions import NotReachable # noqa +from orbit_predictor.predictors.tle import TLEPredictor # noqa diff --git a/orbit_predictor/predictors/base.py b/orbit_predictor/predictors/base.py new file mode 100644 index 0000000..818cd38 --- /dev/null +++ b/orbit_predictor/predictors/base.py @@ -0,0 +1,135 @@ +# -*- coding: utf8 -*- + +import datetime +import logging +import warnings +from collections import namedtuple +from math import acos, degrees + +from orbit_predictor import coordinate_systems +from orbit_predictor.utils import ( + cross_product, + dot_product, + reify, + vector_diff, + vector_norm, + gstime_from_datetime) + +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 CartesianPredictor(Predictor): + + 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 = 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) + + +class GPSPredictor(Predictor): + pass diff --git a/orbit_predictor/predictors/tle.py b/orbit_predictor/predictors/tle.py new file mode 100644 index 0000000..6868d96 --- /dev/null +++ b/orbit_predictor/predictors/tle.py @@ -0,0 +1,91 @@ +import datetime +from math import degrees + +from sgp4.earth_gravity import wgs84 +from sgp4.io import twoline2rv + +from orbit_predictor.exceptions import NotReachable +from orbit_predictor.locations import Location +from orbit_predictor.predictors import PredictedPass +from orbit_predictor.predictors.base import CartesianPredictor, logger + + +class TLEPredictor(CartesianPredictor): + + 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 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_ diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index 142edd2..c2fc9f6 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -25,6 +25,9 @@ from datetime import datetime from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt +from sgp4.ext import jday +from sgp4.propagation import _gstime + try: from functools import lru_cache except ImportError: @@ -185,6 +188,11 @@ def sidereal_time(utc_tuple, local_lon, sun_lon): return GMST0 + UTH + local_lon / 15 +def gstime_from_datetime(when_utc): + timetuple = when_utc.timetuple()[:6] + return _gstime(jday(*timetuple)) + + class reify(object): """ Use as a class method decorator. It operates almost exactly like the From d93ff1a6fd48d53341721abdd72d8d211df25089 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Mon, 13 Nov 2017 17:15:47 +0100 Subject: [PATCH 2/6] Add Keplerian propagator --- orbit_predictor/angles.py | 156 ++++++++++++++++++++++++ orbit_predictor/keplerian.py | 105 ++++++++++++++++ orbit_predictor/predictors/keplerian.py | 121 ++++++++++++++++++ orbit_predictor/utils.py | 61 +++++++++ setup.py | 1 + tests/GMAT_SimplePropagation.script | 49 ++++++++ tests/test_angles.py | 70 +++++++++++ tests/test_keplerian.py | 54 ++++++++ tests/test_keplerian_predictor.py | 117 ++++++++++++++++++ 9 files changed, 734 insertions(+) create mode 100644 orbit_predictor/angles.py create mode 100644 orbit_predictor/keplerian.py create mode 100644 orbit_predictor/predictors/keplerian.py create mode 100644 tests/GMAT_SimplePropagation.script create mode 100644 tests/test_angles.py create mode 100644 tests/test_keplerian.py create mode 100644 tests/test_keplerian_predictor.py diff --git a/orbit_predictor/angles.py b/orbit_predictor/angles.py new file mode 100644 index 0000000..d56fad9 --- /dev/null +++ b/orbit_predictor/angles.py @@ -0,0 +1,156 @@ +# coding: utf-8 +# Inspired by https://github.com/poliastro/poliastro/blob/86f971c/src/poliastro/twobody/angles.py +# Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license +"""Angles and anomalies. + +""" +from math import sin, cos, tan, atan, sqrt + +from numpy import isclose + + +def _kepler_equation(E, M, ecc): + return E - ecc * sin(E) - M + + +def _kepler_equation_prime(E, _, ecc): + return 1 - ecc * cos(E) + + +def ta_to_E(ta, ecc): + """Eccentric anomaly from true anomaly. + + Parameters + ---------- + ta : float + True anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + E : float + Eccentric anomaly. + + """ + E = 2 * atan(sqrt((1 - ecc) / (1 + ecc)) * tan(ta / 2)) + return E + + +def E_to_ta(E, ecc): + """True anomaly from eccentric anomaly. + + Parameters + ---------- + E : float + Eccentric anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + ta : float + True anomaly (rad). + + """ + ta = 2 * atan(sqrt((1 + ecc) / (1 - ecc)) * tan(E / 2)) + return ta + + +def M_to_E(M, ecc): + """Eccentric anomaly from mean anomaly. + + Parameters + ---------- + M : float + Mean anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + E : float + Eccentric anomaly. + + Note + ----- + Algorithm taken from Vallado 2007, pp. 73. + + """ + E = M + while True: + E_new = E + (M - E + ecc * sin(E)) / (1 - ecc * cos(E)) + if isclose(E_new, E, rtol=1e-15): + break + else: + E = E_new + + return E_new + + +def E_to_M(E, ecc): + """Mean anomaly from eccentric anomaly. + + Parameters + ---------- + E : float + Eccentric anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + M : float + Mean anomaly (rad). + + """ + M = _kepler_equation(E, 0.0, ecc) + return M + + +def M_to_ta(M, ecc): + """True anomaly from mean anomaly. + + Parameters + ---------- + M : float + Mean anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + ta : float + True anomaly (rad). + + Examples + -------- + >>> ta = M_to_ta(radians(30.0), 0.06) + >>> rad2deg(ta) + 33.673284930211658 + + """ + E = M_to_E(M, ecc) + ta = E_to_ta(E, ecc) + return ta + + +def ta_to_M(ta, ecc): + """Mean anomaly from true anomaly. + + Parameters + ---------- + ta : float + True anomaly (rad). + ecc : float + Eccentricity. + + Returns + ------- + M : float + Mean anomaly (rad). + + """ + E = ta_to_E(ta, ecc) + M = E_to_M(E, ecc) + return M diff --git a/orbit_predictor/keplerian.py b/orbit_predictor/keplerian.py new file mode 100644 index 0000000..f3e8a3f --- /dev/null +++ b/orbit_predictor/keplerian.py @@ -0,0 +1,105 @@ +# coding: utf-8 +# Inspired by +# https://github.com/poliastro/poliastro/blob/1d2f3ca/src/poliastro/twobody/classical.py +# https://github.com/poliastro/poliastro/blob/1d2f3ca/src/poliastro/twobody/rv.py +# Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license +from math import cos, sin, sqrt + +import numpy as np +from numpy.linalg import norm + +from orbit_predictor.utils import transform + + +def rv_pqw(k, p, ecc, nu): + """Returns r and v vectors in perifocal frame. + + """ + position_pqw = np.array([cos(nu), sin(nu), 0]) * p / (1 + ecc * cos(nu)) + velocity_pqw = np.array([-sin(nu), (ecc + cos(nu)), 0]) * sqrt(k / p) + + return position_pqw, velocity_pqw + + +def coe2rv(k, p, ecc, inc, raan, argp, ta): + """Converts from classical orbital elements to vectors. + + Parameters + ---------- + k : float + Standard gravitational parameter (km^3 / s^2). + p : float + Semi-latus rectum or parameter (km). + ecc : float + Eccentricity. + inc : float + Inclination (rad). + raan : float + Longitude of ascending node (rad). + argp : float + Argument of perigee (rad). + ta : float + True anomaly (rad). + + """ + position_pqw, velocity_pqw = rv_pqw(k, p, ecc, ta) + + position_eci = transform(np.array(position_pqw), 'z', -argp) + position_eci = transform(position_eci, 'x', -inc) + position_eci = transform(position_eci, 'z', -raan) + + velocity_eci = transform(np.array(velocity_pqw), 'z', -argp) + velocity_eci = transform(velocity_eci, 'x', -inc) + velocity_eci = transform(velocity_eci, 'z', -raan) + + return position_eci, velocity_eci + + +def rv2coe(k, r, v, tol=1e-8): + """Converts from vectors to classical orbital elements. + + Parameters + ---------- + k : float + Standard gravitational parameter (km^3 / s^2). + r : ndarray + Position vector (km). + v : ndarray + Velocity vector (km / s). + tol : float, optional + Tolerance for eccentricity and inclination checks, default to 1e-8. + + """ + h = np.cross(r, v) + n = np.cross([0, 0, 1], h) / norm(h) + e = ((v.dot(v) - k / (norm(r))) * r - r.dot(v) * v) / k + ecc = norm(e) + p = h.dot(h) / k + inc = np.arccos(h[2] / norm(h)) + + circular = ecc < tol + equatorial = abs(inc) < tol + + if equatorial and not circular: + raan = 0 + argp = np.arctan2(e[1], e[0]) % (2 * np.pi) # Longitude of periapsis + ta = (np.arctan2(h.dot(np.cross(e, r)) / norm(h), r.dot(e)) % + (2 * np.pi)) + elif not equatorial and circular: + raan = np.arctan2(n[1], n[0]) % (2 * np.pi) + argp = 0 + # Argument of latitude + ta = (np.arctan2(r.dot(np.cross(h, n)) / norm(h), r.dot(n)) % + (2 * np.pi)) + elif equatorial and circular: + raan = 0 + argp = 0 + ta = np.arctan2(r[1], r[0]) % (2 * np.pi) # True longitude + else: + raan = np.arctan2(n[1], n[0]) % (2 * np.pi) + argp = (np.arctan2(e.dot(np.cross(h, n)) / norm(h), e.dot(n)) % + (2 * np.pi)) + ta = (np.arctan2(r.dot(np.cross(h, e)) / norm(h), r.dot(e)) + % (2 * np.pi)) + + return p, ecc, inc, raan, argp, ta diff --git a/orbit_predictor/predictors/keplerian.py b/orbit_predictor/predictors/keplerian.py new file mode 100644 index 0000000..da72394 --- /dev/null +++ b/orbit_predictor/predictors/keplerian.py @@ -0,0 +1,121 @@ +# coding: utf-8 +import datetime +from math import degrees, radians, sqrt + +import numpy as np +from sgp4.earth_gravity import wgs84 +from sgp4.io import twoline2rv + +from orbit_predictor import coordinate_systems +from orbit_predictor.angles import ta_to_M, M_to_ta +from orbit_predictor.keplerian import rv2coe, coe2rv +from orbit_predictor.predictors import TLEPredictor +from orbit_predictor.predictors.base import CartesianPredictor +from orbit_predictor.utils import gstime_from_datetime + +MU_E = wgs84.mu + + +def kepler(argp, delta_t_sec, ecc, inc, p, raan, sma, ta): + # Initial mean anomaly + M_0 = ta_to_M(ta, ecc) + + # Mean motion + n = sqrt(wgs84.mu / sma ** 3) + + # Propagation + M = M_0 + n * delta_t_sec + + # New true anomaly + ta = M_to_ta(M, ecc) + + # Position and velocity vectors + position_eci, velocity_eci = coe2rv(MU_E, p, ecc, inc, raan, argp, ta) + + return position_eci, velocity_eci + + +class KeplerianPredictor(CartesianPredictor): + """Propagator that uses the Keplerian osculating orbital elements. + + We use a naïve propagation algorithm that advances the anomaly the + corresponding amount depending on the time difference and keeps all + the rest of the osculating elements. It's robust against singularities + as long as the starting elements are well specified but only works for + elliptical orbits (ecc < 1). This limitation is not a problem since the object + of study are artificial satellites orbiting the Earth. + + """ + def __init__(self, sma, ecc, inc, raan, argp, ta, epoch): + """Initializes predictor. + + :param sma: Semimajor axis, km + :param ecc: Eccentricity + :param inc: Inclination, deg + :param raan: Right ascension of the ascending node, deg + :param argp: Argument of perigee, deg + :param ta: True anomaly, deg + :param epoch: Epoch, datetime + """ + if ecc >= 1.0: + raise NotImplementedError("Parabolic and elliptic orbits" + "are not implemented") + self._sma = sma + self._ecc = ecc + self._inc = inc + self._raan = raan + self._argp = argp + self._ta = ta + self._epoch = epoch + + @classmethod + def from_tle(cls, sate_id, source, date=None): + """Returns approximate keplerian elements from TLE. + + The conversion between mean elements in the TEME reference + frame to osculating elements in any standard reference frame + is not well defined in literature (see Vallado 3rd edition, pp 236 to 240) + + """ + # Get latest TLE, or the one corresponding to a specified date + if date is None: + date = datetime.datetime.utcnow() + + tle = source.get_tle(sate_id, date) + + # Retrieve TLE epoch and corresponding position + epoch = twoline2rv(tle.lines[0], tle.lines[1], wgs84).epoch + pos = TLEPredictor(sate_id, source).get_position(epoch) + + # Convert position from ECEF to ECI + gmst = gstime_from_datetime(epoch) + position_eci = coordinate_systems.ecef_to_eci(pos.position_ecef, gmst) + velocity_eci = coordinate_systems.ecef_to_eci(pos.velocity_ecef, gmst) + + # Convert position to Keplerian osculating elements + p, ecc, inc, raan, argp, ta = rv2coe( + wgs84.mu, np.array(position_eci), np.array(velocity_eci)) + sma = p / (1 - ecc ** 2) + + return cls(sma, ecc, degrees(inc), degrees(raan), degrees(argp), degrees(ta), epoch) + + def _propagate_eci(self, when_utc=None): + """Return position and velocity in the given date using ECI coordinate system. + + """ + # Orbit parameters + sma = self._sma + ecc = self._ecc + p = sma * (1 - ecc ** 2) + inc = radians(self._inc) + raan = radians(self._raan) + argp = radians(self._argp) + ta = radians(self._ta) + + # Time increment + delta_t_sec = (when_utc - self._epoch).total_seconds() + + # Propagate + position_eci, velocity_eci = kepler(argp, delta_t_sec, ecc, inc, p, raan, sma, ta) + + return tuple(position_eci), tuple(velocity_eci) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index c2fc9f6..017c6aa 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -25,6 +25,7 @@ from datetime import datetime from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt +import numpy as np from sgp4.ext import jday from sgp4.propagation import _gstime @@ -90,6 +91,66 @@ def vector_norm(a): return euclidean_distance(*a) +# Inspired by https://github.com/poliastro/poliastro/blob/aaa1bb2/poliastro/util.py +# and https://github.com/poliastro/poliastro/blob/06ef6ba/poliastro/util.py +# Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license +def rotate(vec, ax, angle): + """Rotates the coordinate system around axis x, y or z a CCW angle. + + Parameters + ---------- + vec : array + Dimension 3 vector. + ax : int + Axis to be rotated. + angle : float + Angle of rotation (rad). + + Notes + ----- + This performs a so-called active or alibi transformation: rotates the + vector while the coordinate system remains unchanged. To do the opposite + operation (passive or alias transformation) call the function as + `rotate(vec, ax, -angle)` or use the convenience function `transform`, + see `[1]_`. + + References + ---------- + .. [1] http://en.wikipedia.org/wiki/Rotation_matrix#Ambiguities + + """ + assert vec.shape == (3,) + + rot = np.eye(3) + if ax == 'x': + sl = slice(1, 3) + elif ax == 'y': + sl = slice(0, 3, 2) + elif ax == 'z': + sl = slice(0, 2) + else: + raise ValueError("Invalid axis: must be one of 'x', 'y' or 'z'") + + rot[sl, sl] = np.array([ + [np.cos(angle), -np.sin(angle)], + [np.sin(angle), np.cos(angle)] + ]) + + return np.dot(rot, vec) + + +def transform(vec, ax, angle): + """Rotates a coordinate system around axis a positive right-handed angle. + + Notes + ----- + This is a convenience function, equivalent to `rotate(vec, ax, -angle)`. + Refer to the documentation of that function for further information. + + """ + return rotate(vec, ax, -angle) + + def sun_azimuth_elevation(latitude_deg, longitude_deg, when=None): """ Return (azimuth, elevation) of the Sun at ground point diff --git a/setup.py b/setup.py index 5bfef57..c58659f 100755 --- a/setup.py +++ b/setup.py @@ -35,6 +35,7 @@ 'Programming Language :: Python :: 3', ], install_requires=[ + 'numpy>=1.8.2', 'sgp4', 'requests', ], diff --git a/tests/GMAT_SimplePropagation.script b/tests/GMAT_SimplePropagation.script new file mode 100644 index 0000000..93c1b89 --- /dev/null +++ b/tests/GMAT_SimplePropagation.script @@ -0,0 +1,49 @@ +%General Mission Analysis Tool(GMAT) Script +%Author: Juan Luis Cano + +%---------- Spacecraft + +Create Spacecraft geoSat; +GMAT geoSat.DateFormat = TAIModJulian; +GMAT geoSat.Epoch = '21545'; +GMAT geoSat.CoordinateSystem = EarthMJ2000Eq; +GMAT geoSat.DisplayStateType = Keplerian; +GMAT geoSat.SMA = 6780; +GMAT geoSat.ECC = 0.001; +GMAT geoSat.INC = 28.5; +GMAT geoSat.RAAN = 67.0; +GMAT geoSat.AOP = 355.0; +GMAT geoSat.TA = 250.0; +GMAT geoSat.Id = 'SatId'; + +%---------- ForceModels + +Create ForceModel AllForces; +GMAT AllForces.CentralBody = Earth; +GMAT AllForces.PrimaryBodies = {Earth}; +GMAT AllForces.Drag = None; +GMAT AllForces.SRP = Off; +GMAT AllForces.RelativisticCorrection = Off; +GMAT AllForces.ErrorControl = None; +GMAT AllForces.GravityField.Earth.Degree = 2; +GMAT AllForces.GravityField.Earth.Order = 0; +GMAT AllForces.GravityField.Earth.PotentialFile = 'EGM96.cof'; +GMAT AllForces.GravityField.Earth.EarthTideModel = 'None'; + +%---------- Propagators + +Create Propagator DefaultProp; +GMAT DefaultProp.FM = AllForces; +GMAT DefaultProp.Type = PrinceDormand853; +GMAT DefaultProp.InitialStepSize = 60; +GMAT DefaultProp.Accuracy = 1.0e-011; +GMAT DefaultProp.MinStep = 0.0; +GMAT DefaultProp.MaxStep = 86400; +GMAT DefaultProp.MaxStepAttempts = 1000; +GMAT DefaultProp.StopIfAccuracyIsViolated = true; + +%---------- Mission Sequence + +BeginMissionSequence; + +Propagate 'Prop 3 hours' DefaultProp(geoSat) {geoSat.ElapsedDays = 0.125}; diff --git a/tests/test_angles.py b/tests/test_angles.py new file mode 100644 index 0000000..c4ff806 --- /dev/null +++ b/tests/test_angles.py @@ -0,0 +1,70 @@ +# coding: utf-8 +# Inspired by +# https://github.com/poliastro/poliastro/blob/86f971c/src/poliastro/twobody/tests/test_angles.py +# Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license +from unittest import TestCase + +from math import radians, degrees + +from orbit_predictor import angles + + +class AnglesTests(TestCase): + def test_true_to_eccentric(self): + # Data from NASA-TR-R-158 + data = [ + # ecc, E (deg), ta(deg) + (0.0, 0.0, 0.0), + (0.05, 10.52321, 11.05994), + (0.10, 54.67466, 59.49810), + (0.35, 142.27123, 153.32411), + (0.61, 161.87359, 171.02189) + ] + for row in data: + ecc, expected_E, ta = row + + E = angles.ta_to_E(radians(ta), ecc) + + self.assertAlmostEqual(degrees(E), expected_E, places=4) + + def test_mean_to_true(self): + # Data from Schlesinger & Udick, 1912 + data = [ + # ecc, M (deg), ta (deg) + (0.0, 0.0, 0.0), + (0.05, 10.0, 11.06), + (0.06, 30.0, 33.67), + (0.04, 120.0, 123.87), + (0.14, 65.0, 80.50), + (0.19, 21.0, 30.94), + (0.35, 65.0, 105.71), + (0.48, 180.0, 180.0), + (0.75, 125.0, 167.57) + ] + for row in data: + ecc, M, expected_ta = row + + ta = angles.M_to_ta(radians(M), ecc) + + self.assertAlmostEqual(degrees(ta), expected_ta, places=2) + + def test_true_to_mean(self): + # Data from Schlesinger & Udick, 1912 + data = [ + # ecc, M (deg), ta (deg) + (0.0, 0.0, 0.0), + (0.05, 10.0, 11.06), + (0.06, 30.0, 33.67), + (0.04, 120.0, 123.87), + (0.14, 65.0, 80.50), + (0.19, 21.0, 30.94), + (0.35, 65.0, 105.71), + (0.48, 180.0, 180.0), + (0.75, 125.0, 167.57) + ] + for row in data: + ecc, expected_M, ta = row + + M = angles.ta_to_M(radians(ta), ecc) + + self.assertAlmostEqual(degrees(M), expected_M, places=1) diff --git a/tests/test_keplerian.py b/tests/test_keplerian.py new file mode 100644 index 0000000..f67953f --- /dev/null +++ b/tests/test_keplerian.py @@ -0,0 +1,54 @@ +# coding: utf-8 +# Inspired by +# Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license +from unittest import TestCase + +from math import radians +import numpy as np +from numpy.testing import assert_allclose + +from sgp4.earth_gravity import wgs84 + +from orbit_predictor.keplerian import coe2rv, rv2coe + + +class COE2RVTests(TestCase): + def test_convert_coe_to_rv(self): + # Data from Vallado, example 2.6 + p = 11067.790 + ecc = 0.83285 + inc = radians(87.87) + raan = radians(227.89) + argp = radians(53.38) + ta = radians(92.335) + + expected_position = [6525.344, 6861.535, 6449.125] + expected_velocity = [4.902276, 5.533124, -1.975709] + + position, velocity = coe2rv(wgs84.mu, p, ecc, inc, raan, argp, ta) + + assert_allclose(position, expected_position, rtol=1e-5) + assert_allclose(velocity, expected_velocity, rtol=1e-5) + + +class RV2COETests(TestCase): + def test_convert_rv_to_coe(self): + # Data from Vallado, example 2.5 + position = np.array([6524.384, 6862.875, 6448.296]) + velocity = np.array([4.901327, 5.533756, -1.976341]) + + expected_p = 11067.79 + expected_ecc = 0.832853 + expected_inc = radians(87.870) + expected_raan = radians(227.89) + expected_argp = radians(53.38) + expected_ta = radians(92.335) + + p, ecc, inc, raan, argp, ta = rv2coe(wgs84.mu, position, velocity) + + self.assertAlmostEqual(p, expected_p, places=0) + self.assertAlmostEqual(ecc, expected_ecc, places=4) + self.assertAlmostEqual(inc, expected_inc, places=4) + self.assertAlmostEqual(raan, expected_raan, places=3) + self.assertAlmostEqual(argp, expected_argp, places=3) + self.assertAlmostEqual(ta, expected_ta, places=5) diff --git a/tests/test_keplerian_predictor.py b/tests/test_keplerian_predictor.py new file mode 100644 index 0000000..964b530 --- /dev/null +++ b/tests/test_keplerian_predictor.py @@ -0,0 +1,117 @@ +from unittest import TestCase + +from datetime import datetime, timedelta + +import numpy as np +from numpy.testing import assert_allclose + +from orbit_predictor.predictors import TLEPredictor +from orbit_predictor.predictors.keplerian import KeplerianPredictor +from orbit_predictor.sources import MemoryTLESource + + +class KeplerianPredictorTests(TestCase): + def setUp(self): + # Data from Vallado, example 2.4 + # Converted to classical orbital elements + sma = 7200.478692389954 + ecc = 0.00810123424807035 + inc = 98.59998936154028 + raan = 319.7043176816153 + argp = 70.87958362589532 + ta = 0.004121614893481961 + + self.epoch = datetime(2000, 1, 1, 12, 0) + + self.predictor = KeplerianPredictor(sma, ecc, inc, raan, argp, ta, self.epoch) + + def test_initial_position(self): + expected_position = np.array([1131.340, -2282.343, 6672.423]) + expected_velocity = np.array([-5.64305, 4.30333, 2.42879]) + + position_eci, velocity_eci = self.predictor._propagate_eci(self.epoch) + + assert_allclose(position_eci, expected_position, rtol=1e-7) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-6) + + def test_propagate_eci(self): + expected_position = np.array([-4219.7527, 4363.0292, -3958.7666]) + expected_velocity = np.array([3.689866, -1.916735, -6.112511]) + + when_utc = self.epoch + timedelta(minutes=40) + + position_eci, velocity_eci = self.predictor._propagate_eci(when_utc) + + assert_allclose(position_eci, expected_position, rtol=1e-5) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-5) + + +class KeplerianPredictorOneDayTests(TestCase): + def setUp(self): + # Converted to classical orbital elements + sma = 6780 + ecc = 0.001 + inc = 28.5 + raan = 67.0 + argp = 355.0 + ta = 250.0 + + self.epoch = datetime(2000, 1, 1, 12, 0) + + self.predictor = KeplerianPredictor(sma, ecc, inc, raan, argp, ta, self.epoch) + + def test_initial_position(self): + expected_position = np.array([3852.57404763, -4749.1872318, -2933.02952967]) + expected_velocity = np.array([5.33068317, 5.28723659, -1.54255441]) + + position_eci, velocity_eci = self.predictor._propagate_eci(self.epoch) + + assert_allclose(position_eci, expected_position, rtol=1e-7) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-6) + + def test_propagate_eci(self): + expected_position = np.array([-5154.02724044, 3011.19175291, 3214.77198183]) + expected_velocity = np.array([-3.67652279, -6.71613987, 0.41267465]) + + when_utc = self.epoch + timedelta(hours=23.999999999) + + position_eci, velocity_eci = self.predictor._propagate_eci(when_utc) + + assert_allclose(position_eci, expected_position, rtol=1e-2) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-2) + + def test_propagate_one_day(self): + expected_position = np.array([-5154.02724044, 3011.19175291, 3214.77198183]) + expected_velocity = np.array([-3.67652279, -6.71613987, 0.41267465]) + + when_utc = self.epoch + timedelta(hours=24) + + position_eci, velocity_eci = self.predictor._propagate_eci(when_utc) + + assert_allclose(position_eci, expected_position, rtol=1e-2) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-2) + + +class TLEConversionTests(TestCase): + + 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', + ) + + def test_from_tle_returns_same_initial_conditions_on_epoch(self): + start = datetime(2017, 3, 6, 7, 51) + db = MemoryTLESource() + db.add_tle(self.SATE_ID, self.LINES, start) + + keplerian_predictor = KeplerianPredictor.from_tle(self.SATE_ID, db, start) + tle_predictor = TLEPredictor(self.SATE_ID, db) + + epoch = keplerian_predictor._epoch + + pos_keplerian = keplerian_predictor.get_position(epoch) + pos_tle = tle_predictor.get_position(epoch) + + assert_allclose(pos_keplerian.position_ecef, pos_tle.position_ecef, rtol=1e-11) + assert_allclose(pos_keplerian.velocity_ecef, pos_tle.velocity_ecef, rtol=1e-13) From 5bebefca9894c1d5215210ddd64100f8497c6bd4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Tue, 14 Nov 2017 15:36:01 +0100 Subject: [PATCH 3/6] Add numerical predictor, only J2 for the moment --- orbit_predictor/predictors/numerical.py | 83 +++++++++++++++++++++++++ tests/test_numerical_predictor.py | 35 +++++++++++ 2 files changed, 118 insertions(+) create mode 100644 orbit_predictor/predictors/numerical.py create mode 100644 tests/test_numerical_predictor.py diff --git a/orbit_predictor/predictors/numerical.py b/orbit_predictor/predictors/numerical.py new file mode 100644 index 0000000..40bd0b3 --- /dev/null +++ b/orbit_predictor/predictors/numerical.py @@ -0,0 +1,83 @@ +from math import radians, sqrt, cos, sin + +from sgp4.earth_gravity import wgs84 + +from orbit_predictor.predictors.keplerian import KeplerianPredictor +from orbit_predictor.angles import ta_to_M, M_to_ta +from orbit_predictor.keplerian import coe2rv + +MU_E = wgs84.mu +R_E_KM = wgs84.radiusearthkm +J2 = wgs84.j2 + + +def pkepler(argp, delta_t_sec, ecc, inc, p, raan, sma, ta): + """Perturbed Kepler problem (only J2) + + Notes + ----- + Based on algorithm 64 of Vallado 3rd edition + + """ + # Mean motion + n = sqrt(MU_E / sma ** 3) + + # Initial mean anomaly + M_0 = ta_to_M(ta, ecc) + + # Update for perturbations + delta_raan = ( + - (3 * n * R_E_KM ** 2 * J2) / (2 * p ** 2) * + cos(inc) * delta_t_sec + ) + raan = raan + delta_raan + + delta_argp = ( + (3 * n * R_E_KM ** 2 * J2) / (4 * p ** 2) * + (4 - 5 * sin(inc) ** 2) * delta_t_sec + ) + argp = argp + delta_argp + + M0_dot = ( + (3 * n * R_E_KM ** 2 * J2) / (4 * p ** 2) * + (2 - 3 * sin(inc) ** 2) * sqrt(1 - ecc ** 2) + ) + M_dot = n + M0_dot + + # Propagation + M = M_0 + M_dot * delta_t_sec + + # New true anomaly + ta = M_to_ta(M, ecc) + + # Position and velocity vectors + position_eci, velocity_eci = coe2rv(MU_E, p, ecc, inc, raan, argp, ta) + + return position_eci, velocity_eci + + +class J2Predictor(KeplerianPredictor): + """Propagator that uses secular variations due to J2. + + """ + def _propagate_eci(self, when_utc=None): + """Return position and velocity in the given date using ECI coordinate system. + + """ + # TODO: Remove duplicated code + # Orbit parameters + sma = self._sma + ecc = self._ecc + p = sma * (1 - ecc ** 2) + inc = radians(self._inc) + raan = radians(self._raan) + argp = radians(self._argp) + ta = radians(self._ta) + + # Time increment + delta_t_sec = (when_utc - self._epoch).total_seconds() + + # Propagate + position_eci, velocity_eci = pkepler(argp, delta_t_sec, ecc, inc, p, raan, sma, ta) + + return tuple(position_eci), tuple(velocity_eci) diff --git a/tests/test_numerical_predictor.py b/tests/test_numerical_predictor.py new file mode 100644 index 0000000..a5e6f58 --- /dev/null +++ b/tests/test_numerical_predictor.py @@ -0,0 +1,35 @@ +from datetime import datetime, timedelta + +from unittest import TestCase + +import numpy as np +from numpy.testing import assert_allclose + +from orbit_predictor.predictors.numerical import J2Predictor + + +class J2PredictorTests(TestCase): + def setUp(self): + # Converted to classical orbital elements + sma = 6780 + ecc = 0.001 + inc = 28.5 + raan = 67.0 + argp = 355.0 + ta = 250.0 + + self.epoch = datetime(2000, 1, 1, 12, 0) + + self.predictor = J2Predictor(sma, ecc, inc, raan, argp, ta, self.epoch) + + def test_propagate_eci(self): + # Data from GMAT + expected_position = np.array([2085.9287615146, -6009.5713894563, -2357.3802307070]) + expected_velocity = np.array([6.4787522759177, 3.2366136616580, -2.5063420188165]) + + when_utc = self.epoch + timedelta(hours=3) + + position_eci, velocity_eci = self.predictor._propagate_eci(when_utc) + + assert_allclose(position_eci, expected_position, rtol=1e-2) + assert_allclose(velocity_eci, expected_velocity, rtol=1e-2) From e9bff5bb73dd8c30e28f2e3e889206df59741ab3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Tue, 14 Nov 2017 15:58:47 +0100 Subject: [PATCH 4/6] Add license and encoding where missing --- orbit_predictor/angles.py | 24 +++++++++++++++++++++++- orbit_predictor/keplerian.py | 24 +++++++++++++++++++++++- orbit_predictor/predictors/base.py | 23 ++++++++++++++++++++++- orbit_predictor/predictors/keplerian.py | 24 +++++++++++++++++++++++- orbit_predictor/predictors/numerical.py | 23 +++++++++++++++++++++++ orbit_predictor/predictors/tle.py | 23 +++++++++++++++++++++++ orbit_predictor/utils.py | 1 + 7 files changed, 138 insertions(+), 4 deletions(-) diff --git a/orbit_predictor/angles.py b/orbit_predictor/angles.py index d56fad9..0583ff3 100644 --- a/orbit_predictor/angles.py +++ b/orbit_predictor/angles.py @@ -1,4 +1,26 @@ -# coding: utf-8 +# -*- 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. +# # Inspired by https://github.com/poliastro/poliastro/blob/86f971c/src/poliastro/twobody/angles.py # Copyright (c) 2012-2017 Juan Luis Cano Rodríguez, MIT license """Angles and anomalies. diff --git a/orbit_predictor/keplerian.py b/orbit_predictor/keplerian.py index f3e8a3f..977e800 100644 --- a/orbit_predictor/keplerian.py +++ b/orbit_predictor/keplerian.py @@ -1,4 +1,26 @@ -# coding: utf-8 +# -*- 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. +# # Inspired by # https://github.com/poliastro/poliastro/blob/1d2f3ca/src/poliastro/twobody/classical.py # https://github.com/poliastro/poliastro/blob/1d2f3ca/src/poliastro/twobody/rv.py diff --git a/orbit_predictor/predictors/base.py b/orbit_predictor/predictors/base.py index 818cd38..9470bf0 100644 --- a/orbit_predictor/predictors/base.py +++ b/orbit_predictor/predictors/base.py @@ -1,4 +1,25 @@ -# -*- coding: utf8 -*- +# -*- 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 diff --git a/orbit_predictor/predictors/keplerian.py b/orbit_predictor/predictors/keplerian.py index da72394..eb24189 100644 --- a/orbit_predictor/predictors/keplerian.py +++ b/orbit_predictor/predictors/keplerian.py @@ -1,4 +1,26 @@ -# coding: utf-8 +# -*- 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 from math import degrees, radians, sqrt diff --git a/orbit_predictor/predictors/numerical.py b/orbit_predictor/predictors/numerical.py index 40bd0b3..38dc1d2 100644 --- a/orbit_predictor/predictors/numerical.py +++ b/orbit_predictor/predictors/numerical.py @@ -1,3 +1,26 @@ +# -*- 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. + from math import radians, sqrt, cos, sin from sgp4.earth_gravity import wgs84 diff --git a/orbit_predictor/predictors/tle.py b/orbit_predictor/predictors/tle.py index 6868d96..44f7bd4 100644 --- a/orbit_predictor/predictors/tle.py +++ b/orbit_predictor/predictors/tle.py @@ -1,3 +1,26 @@ +# -*- 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 from math import degrees diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index 017c6aa..44655f6 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- # MIT License # # Copyright (c) 2017 Satellogic SA From 7ce6f63cdbb884ebb0b63842905f39abe4ac508f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Tue, 14 Nov 2017 21:45:15 +0100 Subject: [PATCH 5/6] Backport lru_cache for older Python versions --- orbit_predictor/utils.py | 8 +------- setup.py | 3 ++- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index 44655f6..e26cb71 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -33,13 +33,7 @@ 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 + from backports.functools_lru_cache import lru_cache # This function was ported from its Matlab equivalent here: # http://www.mathworks.com/matlabcentral/fileexchange/23051-vectorized-solar-azimuth-and-elevation-estimation diff --git a/setup.py b/setup.py index c58659f..71d6b39 100755 --- a/setup.py +++ b/setup.py @@ -49,6 +49,7 @@ "pytest", "pytest-cov", "pytz", - ] + ], + ':python_version=="2.7"': ["backports.functools_lru_cache"], } ) From 0a7aa273cf56569b47b035a3d3c41574c12a81a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Juan=20Luis=20Cano=20Rodr=C3=ADguez?= Date: Tue, 14 Nov 2017 22:11:35 +0100 Subject: [PATCH 6/6] Add rotation tests and some fixes --- orbit_predictor/utils.py | 6 ++++-- tests/test_angles.py | 17 +++++++++++++++++ 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/orbit_predictor/utils.py b/orbit_predictor/utils.py index e26cb71..b960e68 100644 --- a/orbit_predictor/utils.py +++ b/orbit_predictor/utils.py @@ -94,9 +94,9 @@ def rotate(vec, ax, angle): Parameters ---------- - vec : array + vec : ndarray Dimension 3 vector. - ax : int + ax : str Axis to be rotated. angle : float Angle of rotation (rad). @@ -130,6 +130,8 @@ def rotate(vec, ax, angle): [np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)] ]) + if ax == 'y': + rot = rot.T return np.dot(rot, vec) diff --git a/tests/test_angles.py b/tests/test_angles.py index c4ff806..ecdea6c 100644 --- a/tests/test_angles.py +++ b/tests/test_angles.py @@ -6,7 +6,11 @@ from math import radians, degrees +import numpy as np +from numpy.testing import assert_allclose + from orbit_predictor import angles +from orbit_predictor.utils import rotate class AnglesTests(TestCase): @@ -68,3 +72,16 @@ def test_true_to_mean(self): M = angles.ta_to_M(radians(ta), ecc) self.assertAlmostEqual(degrees(M), expected_M, places=1) + + +class RotateTests(TestCase): + def test_rotate_simple(self): + vec = np.array([1, 0, 0]) + + assert_allclose(rotate(vec, 'x', np.radians(90)), np.array([1, 0, 0]), atol=1e-16) + assert_allclose(rotate(vec, 'y', np.radians(90)), np.array([0, 0, -1]), atol=1e-16) + assert_allclose(rotate(vec, 'z', np.radians(90)), np.array([0, 1, 0]), atol=1e-16) + + def test_rotate_raises_error(self): + vec_unused = np.ones(3) + self.assertRaises(ValueError, rotate, vec_unused, "q", 0)