ACROSS API: Add SAA and EarthLimb constraint classes

python/across_api/base/constraints.py

@@ -0,0 +1,167 @@
+# Copyright © 2023 United States Government as represented by the
+# Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+
+from typing import Union
+
+import astropy.units as u  # type: ignore[import]
+import numpy as np
+from astropy.coordinates import Angle, SkyCoord  # type: ignore[import]
+from astropy.time import Time  # type: ignore[import]
+from shapely import Polygon, points  # type: ignore[import]
+
+from .ephem import EphemBase
+
+
+def get_slice(time: Time, ephem: EphemBase) -> slice:
+    """
+    Return a slice for what the part of the ephemeris that we're using.
+
+    Arguments
+    ---------
+    time
+        The time to calculate the slice for
+    ephem
+        The spacecraft ephemeris
+
+    Returns
+    -------
+        The slice for the ephemeris
+    """
+    # If we're just passing a single time, we can just find the index for that
+    if time.isscalar:
+        # Check that the time is within the ephemeris range, this should never
+        # happen as the ephemeris is generated based on this range.
+        assert (
+            time >= ephem.begin and time <= ephem.end
+        ), "Time outside of ephemeris of range"
+        # Find the index for the time and return a slice for that index
+        index = ephem.ephindex(time)
+        return slice(index, index + 1)
+    else:
+        # Check that the time range is within the ephemeris range, as above.
+        assert (
+            time[0] >= ephem.begin and time[-1] <= ephem.end
+        ), "Time outside of ephemeris of range"
+
+        # Find the indices for the start and end of the time range and return a
+        # slice for that range
+        return slice(ephem.ephindex(time[0]), ephem.ephindex(time[-1]) + 1)
+
+
+class SAAPolygonConstraint:
+    """
+    Polygon based SAA constraint. The SAA is defined by a Shapely Polygon, and
+    this constraint will calculate for a given set of times and a given
+    ephemeris whether the spacecraft is in that SAA polygon.
+
+    Attributes
+    ----------
+    polygon
+        Shapely Polygon object defining the SAA polygon.
+    """
+
+    polygon: Polygon
+
+    def __init__(self, polygon: list):
+        self.polygon = Polygon(polygon)
+
+    def __call__(self, time: Time, ephem: EphemBase) -> Union[bool, np.ndarray]:
+        """
+        Return a bool array indicating whether the spacecraft is in constraint
+        for a given ephemeris.
+
+        Arguments
+        ---------
+        time
+            The time(s) to calculate the constraint for.
+        ephem
+            The spacecraft ephemeris, must be precalculated. NOTE: The
+            ephemeris can be calculated for a longer time range than the `time`
+            argument, but it must contain the time(s) in the `time` argument.
+
+        Returns
+        -------
+            Array of booleans for every value in `time` returning True if the
+            spacecraft is in the SAA polygon at that time. If `time` is a
+            scalar then a single boolean is returned.
+        """
+
+        # Find a slice what the part of the ephemeris that we're using
+        i = get_slice(time, ephem)
+
+        in_constraint = self.polygon.contains(
+            points(ephem.longitude[i], ephem.latitude[i])
+        )
+        # Return the result as True or False, or an array of True/False
+        return in_constraint[0] if time.isscalar else in_constraint
+
+
+class EarthLimbConstraint:
+    """
+    For a given Earth limb avoidance angle, is a given coordinate inside this
+    constraint?
+
+    Parameters
+    ----------
+    min_angle
+        The minimum angle from the Earth limb that the spacecraft can point.
+
+    Methods
+    -------
+    __call__(coord, ephem, earthsize=None)
+        Checks if a given coordinate is inside the constraint.
+
+    """
+
+    min_angle: u.Quantity
+
+    def __init__(self, min_angle: u.Quantity):
+        self.min_angle = Angle(min_angle)
+
+    def __call__(
+        self,
+        time: Time,
+        ephem: EphemBase,
+        skycoord: SkyCoord,
+    ) -> Union[bool, np.ndarray]:
+        """
+        Check for a given time, ephemeris and coordinate if positions given are
+        inside the Earth limb constraint. This is done by checking if the
+        separation between the Earth and the spacecraft is less than the
+        Earth's angular radius plus the minimum angle.
+
+        NOTE: Assumes a circular approximation for Earth.
+
+        Parameters
+        ----------
+        coord
+            The coordinate to check.
+        time
+            The time to check.
+        ephem
+            The ephemeris object.
+
+        Returns
+        -------
+        bool
+            `True` if the coordinate is inside the constraint, `False`
+            otherwise.
+
+        """
+        # Find a slice what the part of the ephemeris that we're using
+        i = get_slice(time, ephem)
+
+        # Calculate the angular distance between the center of the Earth and
+        # the object. Note that creating the SkyCoord here from ra/dec stored
+        # in the ephemeris `earth` is 3x faster than just doing the separation
+        # directly with `earth`.
+        in_constraint = (
+            SkyCoord(ephem.earth[i].ra, ephem.earth[i].dec).separation(skycoord)
+            < ephem.earthsize[i] + self.min_angle
+        )
+
+        # Return the result as True or False, or an array of True/False
+        return (
+            in_constraint[0] if time.isscalar and skycoord.isscalar else in_constraint
+        )

python/across_api/burstcube/constraints.py

@@ -0,0 +1,36 @@
+# Copyright © 2023 United States Government as represented by the
+# Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+
+from ..base.constraints import EarthLimbConstraint, SAAPolygonConstraint
+import astropy.units as u  # type: ignore
+
+"""Define constraints for the BurstCube Mission."""
+
+burstcube_saa_constraint = SAAPolygonConstraint(
+    polygon=[
+        (33.900000, -30.0),
+        (12.398, -19.876),
+        (-9.103, -9.733),
+        (-30.605, 0.4),
+        (-38.4, 2.0),
+        (-45.0, 2.0),
+        (-65.0, -1.0),
+        (-84.0, -6.155),
+        (-89.2, -8.880),
+        (-94.3, -14.220),
+        (-94.3, -18.404),
+        (-84.48631, -31.84889),
+        (-86.100000, -30.0),
+        (-72.34921, -43.98599),
+        (-54.5587, -52.5815),
+        (-28.1917, -53.6258),
+        (-0.2095279, -46.88834),
+        (28.8026, -34.0359),
+        (33.900000, -30.0),
+    ]
+)
+
+
+# EarthLimbConstraint
+burstcube_earth_constraint = EarthLimbConstraint(min_angle=0 * u.deg)

python/across_api/swift/constraints.py

@@ -0,0 +1,31 @@
+# Copyright © 2023 United States Government as represented by the
+# Administrator of the National Aeronautics and Space Administration.
+# All Rights Reserved.
+
+import astropy.units as u  # type: ignore
+
+from ..base.constraints import EarthLimbConstraint, SAAPolygonConstraint
+
+"""Define constraints for the BurstCube Mission."""
+
+# Swift Spacecraft SAA polygon as supplied by Swift team.
+swift_saa_constraint = SAAPolygonConstraint(
+    polygon=[
+        (39.0, -30.0),
+        (36.0, -26.0),
+        (28.0, -21.0),
+        (6.0, -12.0),
+        (-5.0, -6.0),
+        (-21.0, 2.0),
+        (-30.0, 3.0),
+        (-45.0, 2.0),
+        (-60.0, -2.0),
+        (-75.0, -7.0),
+        (-83.0, -10.0),
+        (-87.0, -16.0),
+        (-86.0, -23.0),
+        (-83.0, -30.0),
+    ]
+)
+
+swift_earth_constraint = EarthLimbConstraint(min_angle=33 * u.deg)

python/requirements.in

@@ -5,5 +5,6 @@ email-validator
 fastapi
 mangum
 requests
+shapely
 sgp4
 spacetrack

python/requirements.txt

@@ -59,6 +59,7 @@ numpy==1.26.2
     # via
     #   astropy
     #   pyerfa
+    #   shapely
 packaging==23.2
     # via astropy
 pyasn1==0.5.1
@@ -87,6 +88,8 @@ rush==2021.4.0
     # via spacetrack
 sgp4==2.23
     # via -r requirements.in
+shapely==2.0.2
+    # via -r requirements.in
 simplejson==3.19.2
     # via architect-functions
 six==1.16.0