Merge a4c8c62 into 04c38ae

orbit_predictor/constants.py

@@ -0,0 +1,12 @@
+from math import pi
+
+from sgp4.earth_gravity import wgs84
+
+
+AU = 149597870.700  # km
+
+OMEGA = 2 * pi / (86400 * 365.2421897)  # rad / s
+MU_E = wgs84.mu  # km3 / s2
+R_E_KM = wgs84.radiusearthkm  # km
+J2 = wgs84.j2
+OMEGA_E = 7.292115e-5  # rad / s

orbit_predictor/coordinate_systems.py

@@ -24,6 +24,11 @@
 from math import asin, atan, atan2, cos, degrees, pi, radians, sin, sqrt
 
 
+def _euclidean_distance(*components):
+    # TODO: Remove code duplication with utils
+    return sqrt(sum(c**2 for c in components))
+
+
 def llh_to_ecef(lat_deg, lon_deg, h_km):
     """
     Latitude is geodetic, height is above ellipsoid. Output is in km.
@@ -113,6 +118,21 @@ def ecef_to_eci(eci_coords, gmst):
     return x, y, z
 
 
+def eci_to_radec(eci_coords):
+    xequat, yequat, zequat = eci_coords
+
+    # convert equatorial rectangular coordinates to RA and Decl:
+    r = _euclidean_distance(xequat, yequat, zequat)
+    RA = atan2(yequat, xequat)
+    DEC = asin(zequat/r)
+
+    return RA, DEC, r
+
+
+def radec_to_eci(radec_coords):
+    raise NotImplementedError
+
+
 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)

orbit_predictor/predictors/numerical.py

@@ -24,21 +24,14 @@
 import datetime as dt
 
 import numpy as np
-from sgp4.earth_gravity import wgs84
 
+from orbit_predictor.constants import OMEGA, MU_E, R_E_KM, J2, OMEGA_E
 from orbit_predictor.predictors.keplerian import KeplerianPredictor
 from orbit_predictor.angles import ta_to_M, M_to_ta
 from orbit_predictor.keplerian import coe2rv
 from orbit_predictor.utils import njit, raan_from_ltan, float_to_hms
 
 
-OMEGA = 2 * np.pi / (86400 * 365.2421897)  # rad / s
-MU_E = wgs84.mu  # km3 / s2
-R_E_KM = wgs84.radiusearthkm  # km
-J2 = wgs84.j2
-OMEGA_E = 7.292115e-5  # rad / s
-
-
 def sun_sync_plane_constellation(num_satellites, *,
                                  alt_km=None, ecc=None, inc_deg=None, ltan_h=12, date=None):
     """Creates num_satellites in the same Sun-synchronous plane, uniformly spaced.

orbit_predictor/utils.py

@@ -23,13 +23,16 @@
 import functools
 from collections import namedtuple
 import datetime as dt
-from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt, modf
+from math import asin, atan2, cos, degrees, floor, radians, sin, sqrt, tan, modf
 
 import numpy as np
 from sgp4.earth_gravity import wgs84
 from sgp4.ext import jday
 from sgp4.propagation import _gstime
 
+from .constants import AU
+from .coordinate_systems import eci_to_radec, ecef_to_eci
+
 # Inspired in https://github.com/poliastro/poliastro/blob/88edda8/src/poliastro/jit.py
 try:
     from numba import njit
@@ -73,6 +76,18 @@ def euclidean_distance(*components):
     return sqrt(sum(c**2 for c in components))
 
 
+def angle_between(a, b):
+    """
+    Computes angle between two vectors in degrees.
+
+    Notes
+    -----
+    Naïve algorithm, see https://scicomp.stackexchange.com/q/27689/782.
+
+    """
+    return degrees(np.arccos(dot_product(a, b) / (vector_norm(a) * vector_norm(b))))
+
+
 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))
@@ -175,41 +190,11 @@ def transform(vec, ax, angle):
 
 
 def raan_from_ltan(when, ltan=12.0):
-    # TODO: Avoid code duplication
-    # compute apparent right ascension of the sun (radians)
-    jd = juliandate(timetuple_from_dt(when))
-    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
-    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 equatorial rectangular coordinates
-    xequat = xeclip
-    yequat = yeclip * cos_d(oblecl) + zeclip * sin_d(oblecl)
+    sun_eci = get_sun(when)
 
     # convert equatorial rectangular coordinates to RA and Decl:
-    RA = atan2_d(yequat, xequat)  # degrees
+    RA, _, _ = eci_to_radec(sun_eci)
+    RA = degrees(RA)
 
     # Idea from
     # https://www.mathworks.com/matlabcentral/fileexchange/39085-mean-local-time-of-the-ascending-node
@@ -234,64 +219,147 @@ def sun_azimuth_elevation(latitude_deg, longitude_deg, when=None):
     jd = juliandate(timetuple_from_dt(when))
     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
+    w, M, L, eccentricity, oblecl = _sun_mean_ecliptic_elements(date)
+    sun_eci = _sun_eci(w, M, L, eccentricity, oblecl)
+
+    # convert equatorial rectangular coordinates to RA and Decl:
+    RA, DEC, r = eci_to_radec(sun_eci)
+
+    RA = degrees(RA)
+    DEC = degrees(DEC)
+
+    # 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(DEC)
+    y = sin_d(HA) * cos_d(DEC)
+    z = sin_d(DEC)
+
+    # 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 _sun_mean_ecliptic_elements(t_ut1):
+    w = 282.9404 + 4.70935e-5 * t_ut1    # longitude of perihelion degrees
+    eccentricity = 0.016709 - 1.151e-9 * t_ut1      # eccentricity
+    M = (356.0470 + 0.9856002585 * t_ut1) % 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
+    oblecl = 23.4393 - 3.563e-7 * t_ut1  # Sun's obliquity of the ecliptic
+
+    return w, M, L, eccentricity, oblecl
 
+
+def _sun_eci(w, M, L, eccentricity, oblecl):
     # 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)
+    # rectangular coordinates in the plane of the ecliptic (x axis toward perihelion)
     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
+    # find the true 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
+    # rotate these coordinates to equatorial 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)
+    return [xequat, yequat, zequat]
 
-    # 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)
+def get_sun(when):
+    """
+    Returns inertial position of the Sun, in au.
+    """
+    jd = juliandate(timetuple_from_dt(when))
+    date = jd - DECEMBER_31TH_1999_MIDNIGHT_JD
 
-    # Replace RA with hour angle HA
-    HA = sidereal * 15 - RA
+    w, M, L, eccentricity, oblecl = _sun_mean_ecliptic_elements(date)
+    sun_eci = _sun_eci(w, M, L, eccentricity, oblecl)
 
-    # convert to rectangular coordinate system
-    x = cos_d(HA) * cos_d(delta)
-    y = sin_d(HA) * cos_d(delta)
-    z = sin_d(delta)
+    return np.array(sun_eci)
 
-    # 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)
+def get_shadow(r, when_utc):
+    """
+    Gives illumination of Earth satellite (2 for illuminated, 1 for penumbra, 0 for umbra).
 
-    return AzimuthElevation(azimuth, elevation)
+    Parameters
+    ----------
+    r : numpy.ndarray or list
+        ECEF vector pointing to the satellite in km.
+    when_utc : datetime.datetime
+        Time of calculation.
+
+    """
+    gmst = gstime_from_datetime(when_utc)
+    r_sun = get_sun(when_utc) * AU
+
+    return shadow(r_sun, ecef_to_eci(r, gmst))
+
+
+def shadow(r_sun, r, r_p=wgs84.radiusearthkm):
+    """
+    Gives illumination of Earth satellite (2 for illuminated, 1 for penumbra, 0 for umbra).
+
+    Parameters
+    ----------
+    r_sun : numpy.ndarray or list
+        Vector pointing to the Sun in km.
+    r : numpy.ndarray or list
+        Vector pointing to the satellite in km.
+    r_p : float, optional
+        Radius of the planet, default to Earth WGS84.
+
+    Notes
+    -----
+    Algorithm 34 from Vallado, section 5.3.
+
+    """
+    alpha_umb = radians(0.264121687)
+    alpha_pen = radians(0.269007205)
+
+    if dot_product(r_sun, r) < 0:
+        angle = angle_between(-r_sun, r)
+        sat_horiz = vector_norm(r) * cos_d(angle)
+        sat_vert = vector_norm(r) * sin_d(angle)
+        x = r_p / sin(alpha_pen)
+        pen_vert = tan(alpha_pen) * (x + sat_horiz)
+
+        if sat_vert <= pen_vert:
+            y = r_p / sin(alpha_umb)
+            umb_vert = tan(alpha_umb) * (y - sat_horiz)
+
+            if sat_vert <= umb_vert:
+                return 0
+            else:
+                return 1
+    else:
+        return 2
 
 
 def juliandate(utc_tuple):
@@ -310,7 +378,7 @@ def sidereal_time(utc_tuple, local_lon, sun_lon):
     # 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
+    # Calculate local sidereal time
     GMST0 = ((sun_lon + 180) % 360) / 15
     return GMST0 + UTH + local_lon / 15
 

tests/test_sun.py

@@ -0,0 +1,21 @@
+import datetime as dt
+
+import numpy as np
+import pytest
+
+from orbit_predictor.utils import angle_between, get_sun
+
+
+# Data obtained from Astropy using the JPL ephemerides
+# coords = get_body("sun", Time(when_utc)).represent_as(CartesianRepresentation).xyz.to("au").T.value
+@pytest.mark.parametrize("when_utc,expected_eci", [
+    [dt.datetime(2000, 1, 1, 12), np.array([0.17705013, -0.88744275, -0.38474906])],
+    [dt.datetime(2009, 6, 1, 18, 30), np.array([0.32589889, 0.88109849, 0.38197646])],
+    [dt.datetime(2019, 11, 25, 18, 46, 0), np.array([-0.449363, -0.80638653, -0.34956405])],
+    [dt.datetime(2025, 12, 1, 12), np.array([-0.35042293, -0.84565374, -0.36657211])],
+])
+def test_get_sun_matches_expected_result_within_precision(when_utc, expected_eci):
+    eci = get_sun(when_utc)
+
+    assert angle_between(eci, expected_eci) < 1.0  # Claimed precision
+    assert angle_between(eci, expected_eci) < 0.5  # Actual precision