Comms expansion

.gitignore

@@ -128,6 +128,9 @@ dmypy.json
 # Pyre type checker
 .pyre/
 
+.idea
+pylintrc
+
 my_notebooks
 .vscode
 examples/Sentinel_2_example_notebook/Etna_00.tif

README.md

@@ -328,6 +328,10 @@ print(sat_actor.get_altitude(t0))
 
 #### How to add a communication device
 
+PASEOS supports two types of communication modelling.
+1. A simplified model where a constant bitrate is set between two actors. See the first example below on how to add such a communication device.
+2. A link budget based model where the bitrate is adjusted based on distance between actors and the transmitter, receiver and link characteristics. This other examples in this section will further explore this model.
+
 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.
 
 ```py
@@ -345,6 +349,95 @@ ActorBuilder.add_comm_device(actor=sat_actor,
                              bandwidth_in_kbps=100000)
 ```
 
+Link budget communication modelling is possible between satellites or between a satellite and ground stations. The example below explores an optical connection between sat_actor and comms_sat and a radio connection between comms_sat and maspalomas_groundstation.
+
+The first step is to add an optical transmitter to a spacecraft:
+
+```py
+    import pykep as pk
+    from paseos import SpacecraftActor, ActorBuilder
+    from paseos.communication.device_type import DeviceType
+    sat_actor = ActorBuilder.get_actor_scaffold(name="mySat",
+                                       actor_type=SpacecraftActor,
+                                       epoch=pk.epoch(0))
+
+    ActorBuilder.set_orbit(
+    actor=sat_actor,
+    position=[10000000, 0, 0],
+    velocity=[0, 8000.0, 0],
+    epoch=pk.epoch(0),
+    central_body=pk.planet.jpl_lp("earth"),  # use Earth from pykep
+    )
+    optical_transmitter_name = "sat_optical_transmitter_1"
+    ActorBuilder.add_comm_device(actor=sat_actor,
+                                 device_name=optical_transmitter_name, input_power=1, 
+                                 power_efficiency=1, antenna_efficiency=1, line_losses=1, 
+                                 point_losses=3, fwhm=1E-3, device_type=DeviceType.OPTICAL_TRANSMITTER)
+```
+The full width at half maximum (fwhm) is used to determine the transmitter gain. This can alternatively be determined by setting antenna_gain, or by setting antenna_diameter.
+
+The second step is to add an optical receiver and a radio transmitter to the comms satellite.
+
+```py
+import pykep as pk
+from paseos import SpacecraftActor, ActorBuilder
+from paseos.communication.device_type import DeviceType
+comms_sat:SpacecraftActor = ActorBuilder.get_actor_scaffold(name="comms_1",actor_type=SpacecraftActor, epoch=pk.epoch(0))
+ActorBuilder.set_orbit(
+    actor=comms_sat,
+    position=[10000000, 0, 0],
+    velocity=[0, 8000.0, 0],
+    epoch=pk.epoch(0),
+    central_body=pk.planet.jpl_lp("earth"),  # use Earth from pykep
+    )
+
+optical_receiver_name = "optical_receiver_1"
+ActorBuilder.add_comm_device(actor=comms_sat,device_name=optical_receiver_name, line_losses=4.1, antenna_gain=114.2, device_type=DeviceType.OPTICAL_RECEIVER)
+
+radio_name = "sat_radio_transmitter_1"
+ActorBuilder.add_comm_device(actor=comms_sat,device_name=radio_name, input_power=2, power_efficiency=0.5,
+                                antenna_efficiency=0.5, line_losses=1, point_losses=5, antenna_diameter=0.3, device_type=DeviceType.RADIO_TRANSMITTER)
+```
+Next, add a radio receiver the ground station.
+
+```py
+from paseos import SpacecraftActor, ActorBuilder, GroundstationActor
+maspalomas_groundstation = ActorBuilder.get_actor_scaffold(
+        name="maspalomas_groundstation", actor_type=GroundstationActor, epoch=pk.epoch(0)
+    )
+ActorBuilder.set_ground_station_location(maspalomas_groundstation,latitude=27.7629, longitude=-15.6338, elevation=205.1, minimum_altitude_angle=5)
+receiver_name = "maspalomas_radio_receiver_1"
+ActorBuilder.add_comm_device(actor=maspalomas_groundstation, device_name=receiver_name, noise_temperature=135, 
+                             line_losses=1, polarization_losses=3, antenna_gain=62.6, device_type=DeviceType.RADIO_RECEIVER)
+```
+
+The final step is to add links to the transmitting actors.
+
+```py
+from paseos import SpacecraftActor, ActorBuilder
+optical_link_name = "optical_link_1"
+ActorBuilder.add_comm_link(sat_actor, optical_transmitter_name, comms_sat, optical_receiver_name, optical_link_name)
+
+frequency = 8475E6 #Hz
+radio_link_name = "radio_link_1"
+ActorBuilder.add_comm_link(comms_sat, radio_link_name, maspalomas_groundstation, receiver_name, radio_link_name, frequency=frequency)
+```
+
+The link and bitrate can then be evaluated during the simulation to determine the bitrate:
+
+```py
+from paseos.communication.link_budget_calc import calc_dist_and_alt_angle
+
+while t <= simulation_time:
+    for instance in paseos_instances:
+        local_t = instance.local_actor.local_time
+        for link in instance.local_actor.communication_links:
+            distance, elevation_angle = calc_dist_and_alt_angle(instance.local_actor,
+                                                                link.receiver_actor, local_t)
+            if instance.local_actor.is_in_line_of_sight(link.receiver_actor, epoch=local_t):
+                bitrate = link.get_bitrate(distance, elevation_angle)
+```
+
 #### How to add a power device
 
 The following code snippet shows how to add a power device to a [SpacecraftActor](#spacecraftactor).

examples/Communication_example/communication_example_utils.py

@@ -0,0 +1,123 @@
+import numpy as np
+import pandas as pd
+
+
+def get_known_actor_comms_status(values):
+    """Helper function to track comms status
+    Args:
+        values (list): a list with known actors histories.
+    Returns:
+        A list with status values.
+    """
+    conv_values = []
+    for val in values:
+        status = ["No signal", "Ground only", "CommSat only", "Ground + Sat"]
+        idx = ("comms_1" in val) * 2 + 1 * ("gs_1" in val)
+        conv_values.append(status[idx])
+    return conv_values
+
+
+def get_bitrate_status(link):
+    """Helper function to get bitrate status history
+    Args:
+        link (LinkModel): the link that contains the bitrate history.
+    Returns:
+            The bitrate history, a list with bitrates in bps.
+    """
+    return link.bitrate_history
+
+
+def get_line_of_sight_status(link):
+    """Helper function to get line of sight status history
+    Args:
+        link (LinkModel): the link that contains the line of sight history.
+    Returns:
+            The line of sight history, a list with booleans.
+    """
+    return link.line_of_sight_history
+
+
+def get_distance_status(link):
+    """Helper function to distance status history
+    Args:
+        link (LinkModel): the link that contains the distance history.
+    Returns:
+            The distance history, a list with distances in metres.
+    """
+    return link.distance_history
+
+
+def get_elevation_angle_status(link):
+    """Helper function to get elevation angle status history
+    Args:
+        link (LinkModel): the link that contains the elevation angle history.
+    Returns:
+            The elevation angle history, a list with angles in degrees.
+    """
+    return link.elevation_angle_history
+
+
+def get_closest_entry(df, t, id):
+    """Helper function to the closest entry in the dataframe of a particular time.
+    Args:
+        df (DataFrame): the dataframe.
+        t (df.Time): the timestamp that is being searched for.
+        id (int): the id of the satellite in the dataframe.
+    Returns:
+            The elevation angle history, a list with angles in degrees.
+    """
+    df_id = df[df.ID == id]
+    return df_id.iloc[(df_id["Time"] - t).abs().argsort().iloc[:1]]
+
+
+def get_analysis_df(df, timestep=60, orbital_period=1):
+    """Helper function to get elevation angle status history
+    Args:
+        df (Dataframe): the dataframe to analyze
+        timestep (int): the timestep to construct the analysis
+        orbital_period (int): the orbital period
+    Returns:
+            The dataframe constructed with the analysis outputs.
+    """
+    t = np.round(np.linspace(0, df.Time.max(), int(df.Time.max() // timestep)))
+    sats = df.ID.unique()
+    df["known_actors"] = pd.Categorical(df.known_actors)
+    df["comm_cat"] = df.known_actors.cat.codes
+    standby = []
+    processing = []
+    is_in_eclipse = []
+    comm_stat = [[], [], [], []]
+
+    for idx, t_cur in enumerate(t):
+        n_c = 0
+        n_ec = 0
+        for c in comm_stat:
+            c.append(0)
+        for sat in sats:
+            vals = get_closest_entry(df, t_cur, sat)
+            n_c += vals.current_activity.values[0] == "Standby"
+            n_ec += vals.is_in_eclipse.values[0]
+            c_idx = vals.comm_cat.values[0]
+            comm_stat[c_idx][idx] += 1
+        standby.append(n_c)
+        processing.append(len(sats) - n_c)
+        is_in_eclipse.append(n_ec)
+
+    ana_df = pd.DataFrame(
+        {
+            "Time[s]": t,
+            "# of Standby": standby,
+            "# of Processing": processing,
+            "# in Eclipse": is_in_eclipse,
+            "# of " + df.known_actors.cat.categories[0]: comm_stat[0],
+            "# of " + df.known_actors.cat.categories[1]: comm_stat[1],
+            "# of " + df.known_actors.cat.categories[2]: comm_stat[2],
+            # "# of " + df.known_actors.cat.categories[3]: comm_stat[3],
+        }
+    )
+    ana_df["Completed orbits"] = ana_df["Time[s]"] / orbital_period
+    ana_df = ana_df.round({"Completed orbits": 2})
+    ana_df["Share Processing"] = ana_df["# of Processing"] / len(sats)
+    ana_df["Share in Eclipse"] = ana_df["# in Eclipse"] / len(sats)
+
+    return ana_df