Merge pull request #4 from Juanlu001/keplerian-j2-predictor

Add Keplerian and J2 predictor

orbit_predictor/angles.py

@@ -0,0 +1,178 @@
+# -*- 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.
+
+"""
+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

orbit_predictor/keplerian.py

@@ -0,0 +1,127 @@
+# -*- 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
+# 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

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