Merge pull request #175 from aidotse/allow_using_TLE

Allow using TLEs

README.md

@@ -244,7 +244,11 @@ sim.add_known_actor(grndStation)
 
 #### Set an orbit for a PASEOS SpacecraftActor
 
-Once you have defined a [SpacecraftActor](#spacecraftactor), you can assign a [Keplerian orbit](https://en.wikipedia.org/wiki/Kepler_orbit) to it. To this aim, you need to define the central body the [SpacecraftActor](#spacecraftactor) is orbiting around and specify its position and velocity (in the central body's [inertial frame](https://en.wikipedia.org/wiki/Inertial_frame_of_reference)) and an epoch. In this case, we will use `Earth` as a central body.
+Once you have defined a [SpacecraftActor](#spacecraftactor), you can assign a [Keplerian orbit](https://en.wikipedia.org/wiki/Kepler_orbit) or use [SGP4 (Earth orbit only)](https://en.wikipedia.org/wiki/Simplified_perturbations_models). 
+
+##### Keplerian Orbit
+
+To this aim, you need to define the central body the [SpacecraftActor](#spacecraftactor) is orbiting around and specify its position and velocity (in the central body's [inertial frame](https://en.wikipedia.org/wiki/Inertial_frame_of_reference)) and an epoch. In this case, we will use `Earth` as a central body.
 
 ```py
 import pykep as pk
@@ -264,6 +268,37 @@ ActorBuilder.set_orbit(actor=sat_actor,
                        epoch=pk.epoch(0), central_body=earth)
 ```
 
+##### SGP4 / Two-line element (TLE) 
+
+For using SGP4 / [Two-line element (TLE)](https://en.wikipedia.org/wiki/Two-line_element_set) you need to specify the TLE of the [SpacecraftActor](#spacecraftactor). In this case, we will use the TLE of the [Sentinel-2A](https://en.wikipedia.org/wiki/Sentinel-2) satellite from [celestrak](https://celestrak.com/).
+
+```py
+from paseos import ActorBuilder, SpacecraftActor
+# Define an actor of type SpacecraftActor
+sat_actor = ActorBuilder.get_actor_scaffold(name="Sentinel-2A",
+                                       actor_type=SpacecraftActor,
+                                       epoch=pk.epoch(0))
+
+# Specify your TLE
+line1 = "1 40697U 15028A   23188.15862373  .00000171  00000+0  81941-4 0  9994"
+line2 = "2 40697  98.5695 262.3977 0001349  91.8221 268.3116 14.30817084419867"
+
+# Set the orbit of the actor
+ActorBuilder.set_TLE(sat_actor, line1, line2)
+```
+
+##### Accessing the orbit
+You can access the orbit of a [SpacecraftActor](#spacecraftactor) with
+
+```py
+# Position, velocity and altitude can be accessed like this
+t0 = pk.epoch("2022-06-16 00:00:00.000") # Define the time (epoch)
+print(sat_actor.get_position(t0)) 
+print(sat_actor.get_position_velocity(t0))
+print(sat_actor.get_altitude(t0))
+```
+
+
 #### How to add a communication device
 
 The following code snippet shows how to add a communication device to a [SpacecraftActors] (#spacecraftactor). A communication device is needed to model the communication between [SpacecraftActors] (#spacecraftactor) or a [SpacecraftActor](#spacecraftactor) and [GroundstationActor](#ground-stationactor). Currently, given the maximum transmission data rate of a communication device, PASEOS calculates the maximum data that can be transmitted by multiplying the transmission data rate by the length of the communication window. The latter is calculated by taking the period for which two actors are in line-of-sight into account.

examples/Sentinel_2_example_notebook/Sentinel2_example_notebook.ipynb

@@ -25,7 +25,6 @@
     "import os\n",
     "sys.path.insert(1, os.path.join(\"..\",\"..\"))\n",
     "import pykep as pk\n",
-    "import numpy as np\n",
     "import paseos\n",
     "from paseos import ActorBuilder, SpacecraftActor\n",
     "from utils import s2pix_detector, acquire_data\n",
@@ -55,8 +54,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "#Define today as pykep epoch (27-10-22)\n",
+    "#please, refer to https://esa.github.io/pykep/documentation/core.html#pykep.epoch\n",
+    "today = pk.epoch_from_string('2022-10-27 12:00:00.000')\n",
+    "\n",
     "# Define local actor\n",
-    "S2B = ActorBuilder.get_actor_scaffold(name=\"Sentinel2-B\", actor_type=SpacecraftActor, epoch=pk.epoch(0))"
+    "S2B = ActorBuilder.get_actor_scaffold(name=\"Sentinel2-B\", actor_type=SpacecraftActor, epoch=today)"
    ]
   },
   {
@@ -65,30 +68,7 @@
    "metadata": {},
    "source": [
     "#### 1.a) - Add an orbit for S2B\n",
-    "\n",
-    "Since **S2B** is orbiting around Earth, let's define `earth` as `pykep.planet` object. \n",
-    "\n",
-    "Internally, paseos uses pykep.planet objects to model gravitational acceleration and Keplerian orbits."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "30691b25",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Define central body\n",
-    "earth = pk.planet.jpl_lp(\"earth\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "013bf2e2",
-   "metadata": {},
-   "source": [
-    "To find realistic orbits for **S2B**, we can exploit `Two Line Elements (TLEs)` (Downloaded on 27-10-2022). This would allow finding their ephemerides at time = 27-10-2022 12:00:00.\n",
-    "Please, refer to [Two-line_element_set](https://en.wikipedia.org/wiki/Two-line_element_set) to know more about TLEs."
+    "We can make use of the [two-line element](https://en.wikipedia.org/wiki/Two-line_element_set) actor creation [inside PASEOS](https://github.com/aidotse/PASEOS/tree/allow_using_TLE#sgp4--two-line-element-tle) for set the orbit of the PASEOS actor (TLE Downloaded on 27-10-2022)."
    ]
   },
   {
@@ -98,38 +78,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "#Define today as pykep epoch (27-10-22)\n",
-    "#please, refer to https://esa.github.io/pykep/documentation/core.html#pykep.epoch\n",
-    "today = pk.epoch_from_string('2022-10-27 12:00:00.000')\n",
-    "\n",
     "sentinel2B_line1 = \"1 42063U 17013A   22300.18652110  .00000099  00000+0  54271-4 0  9998\"\n",
     "sentinel2B_line2 = \"2 42063  98.5693  13.0364 0001083 104.3232 255.8080 14.30819357294601\"\n",
-    "sentinel2B = pk.planet.tle(sentinel2B_line1, sentinel2B_line2)\n",
     "\n",
-    "#Calculating S2B ephemerides.\n",
-    "sentinel2B_eph=sentinel2B.eph(today)\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "8eab58f6",
-   "metadata": {},
-   "source": [
-    "Now we define the actor's orbit for **S2B**."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "3567ebb8",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "#Adding orbit around Earth based on previously calculated ephemerides\n",
-    "ActorBuilder.set_orbit(actor=S2B, \n",
-    "                       position=sentinel2B_eph[0], \n",
-    "                       velocity=sentinel2B_eph[1], \n",
-    "                       epoch=today, central_body=earth)"
+    "ActorBuilder.set_TLE(S2B, sentinel2B_line1, sentinel2B_line2)"
    ]
   },
   {
@@ -463,7 +415,9 @@
    "source": [
     "### 3.d) - Showing detected volcanic eruptions\n",
     "\n",
-    "The next plot will show an example of onboard coarse volcanic eruptions detection on some Sentinel-2 L1C tiles. The different eruptions will be surrounded a bounding box, and their coordinates will be printed to raise an alert."
+    "The next plot will show an example of onboard coarse volcanic eruptions detection on some Sentinel-2 L1C tiles. The different eruptions will be surrounded a bounding box, and their coordinates will be printed to raise an alert.\n",
+    "\n",
+    "The execution and rendering of the images may take a few minutes."
    ]
   },
   {
@@ -480,7 +434,9 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "85bbd3ac",
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
     "fig, ax=plt.subplots(1,3,figsize=(10,4))\n",
@@ -502,7 +458,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3.10.6 ('paseos')",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },

paseos/actors/actor_builder.py

@@ -16,18 +16,19 @@
 class ActorBuilder:
     """This class is used to construct actors."""
 
-    def __new__(self):
-        if not hasattr(self, "instance"):
-            self.instance = super(ActorBuilder, self).__new__(self)
+    def __new__(cls):
+        if not hasattr(cls, "instance"):
+            cls.instance = super(ActorBuilder, cls).__new__(cls)
         else:
             logger.debug(
                 "Tried to create another instance of ActorBuilder. Keeping original one..."
             )
-        return self.instance
+        return cls.instance
 
     def __init__(self):
         logger.trace("Initializing ActorBuilder")
 
+    @staticmethod
     def get_actor_scaffold(name: str, actor_type: object, epoch: pk.epoch):
         """Initiates an actor with minimal properties.
 
@@ -50,6 +51,7 @@ def get_actor_scaffold(name: str, actor_type: object, epoch: pk.epoch):
 
         return actor_type(name, epoch)
 
+    @staticmethod
     def set_ground_station_location(
         actor: GroundstationActor,
         latitude: float,
@@ -81,6 +83,36 @@ def set_ground_station_location(
         )
         actor._minimum_altitude_angle = minimum_altitude_angle
 
+    @staticmethod
+    def set_TLE(
+        actor: SpacecraftActor,
+        line1: str,
+        line2: str,
+    ):
+        """Define the orbit of the actor using a TLE. For more information on TLEs see
+        https://en.wikipedia.org/wiki/Two-line_element_set .
+
+        TLEs can be obtained from https://www.space-track.org/ or https://celestrak.com/NORAD/elements/
+
+        Args:
+            actor (SpacecraftActor): Actor to update.
+            line1 (str): First line of the TLE.
+            line2 (str): Second line of the TLE.
+
+        Raises:
+            RuntimeError: If the TLE could not be read.
+        """
+        try:
+            actor._orbital_parameters = pk.planet.tle(line1, line2)
+            # TLE only works around Earth
+            actor._central_body = pk.planet.jpl_lp("earth")
+        except RuntimeError:
+            logger.error("Error reading TLE \n", line1, "\n", line2)
+            raise RuntimeError("Error reading TLE")
+
+        logger.debug(f"Added TLE to actor {actor}")
+
+    @staticmethod
     def set_orbit(
         actor: SpacecraftActor,
         position,
@@ -115,6 +147,7 @@ def set_orbit(
 
         logger.debug(f"Added orbit to actor {actor}")
 
+    @staticmethod
     def set_position(actor: BaseActor, position: list):
         """Sets the actors position. Use this if you do not want the actor to have a keplerian orbit around a central body.
 
@@ -133,6 +166,7 @@ def set_position(actor: BaseActor, position: list):
         actor._position = position
         logger.debug(f"Setting position {position} on actor {actor}")
 
+    @staticmethod
     def set_power_devices(
         actor: SpacecraftActor,
         battery_level_in_Ws: float,
@@ -182,6 +216,7 @@ def set_power_devices(
             + f"ChargingRate={charging_rate_in_W}W to actor {actor}"
         )
 
+    @staticmethod
     def set_radiation_model(
         actor: SpacecraftActor,
         data_corruption_events_per_s: float,
@@ -215,6 +250,7 @@ def set_radiation_model(
         )
         logger.debug(f"Added radiation model to actor {actor}.")
 
+    @staticmethod
     def set_thermal_model(
         actor: SpacecraftActor,
         actor_mass: float,
@@ -310,6 +346,7 @@ def set_thermal_model(
             power_consumption_to_heat_ratio=power_consumption_to_heat_ratio,
         )
 
+    @staticmethod
     def add_comm_device(actor: BaseActor, device_name: str, bandwidth_in_kbps: float):
         """Creates a communication device.
 
@@ -327,6 +364,7 @@ def add_comm_device(actor: BaseActor, device_name: str, bandwidth_in_kbps: float
 
         logger.debug(f"Added comm device with bandwith={bandwidth_in_kbps} kbps to actor {actor}.")
 
+    @staticmethod
     def add_custom_property(
         actor: BaseActor, property_name: str, initial_value: Any, update_function: Callable
     ):

paseos/tests/actor_builder_test.py

@@ -5,11 +5,53 @@
 
 sys.path.append("../..")
 
-from paseos import ActorBuilder
+from paseos import ActorBuilder, SpacecraftActor
 
 from test_utils import get_default_instance
 
 
+def test_set_TLE():
+    """Check if we can set a TLE correctly"""
+
+    _, sentinel2a, earth = get_default_instance()
+    # Set the TLE
+    line1 = "1 40697U 15028A   23188.15862373  .00000171  00000+0  81941-4 0  9994"
+    line2 = "2 40697  98.5695 262.3977 0001349  91.8221 268.3116 14.30817084419867"
+    ActorBuilder.set_TLE(sentinel2a, line1, line2)
+
+    # Check that get_altitude returns a sensible value
+    earth_radius = 6371000
+    assert sentinel2a.get_altitude() > earth_radius + 780000
+    assert sentinel2a.get_altitude() < earth_radius + 820000
+
+    # Check that get_position_velocity returns sensible values
+    position, velocity = sentinel2a.get_position_velocity(sentinel2a.local_time)
+    assert position is not None
+    assert velocity is not None
+
+    # Create an actor with a keplerian orbit and check that the position and velocity
+    # diverge over time
+    s2a_kep = ActorBuilder.get_actor_scaffold("s2a_kep", SpacecraftActor, sentinel2a.local_time)
+    ActorBuilder.set_orbit(s2a_kep, position, velocity, sentinel2a.local_time, earth)
+
+    # After some orbits the differences should be significant
+    # since the TLE uses SGP4 and the other actor uses Keplerian elements
+    t0_later = pk.epoch(sentinel2a.local_time.mjd2000 + 1)
+    r, v = sentinel2a.get_position_velocity(t0_later)
+    r_kep, v_kep = s2a_kep.get_position_velocity(t0_later)
+    print("r,v SGP4 after  1 day")
+    print(r)
+    print(v)
+    print("r,v Kep  after  1 day")
+    print(r_kep)
+    print(v_kep)
+    print("Differences in r and v")
+    print(np.linalg.norm(np.array(r) - np.array(r_kep)))
+    print(np.linalg.norm(np.array(v) - np.array(v_kep)))
+    assert np.linalg.norm(np.array(r) - np.array(r_kep)) > 100000
+    assert np.linalg.norm(np.array(v) - np.array(v_kep)) > 400
+
+
 def test_set_orbit():
     """Check if we can specify an orbit correctly"""
     _, sat1, earth = get_default_instance()