From 67ac410ce84854b9395913e4170eea1e531450c3 Mon Sep 17 00:00:00 2001
From: artistrea <apiravesi@gmail.com>
Date: Mon, 3 Nov 2025 23:42:16 -0300
Subject: [PATCH 1/9] feat: minimum elevation for service grid point to
 satellite link

---
 sharc/parameters/imt/parameters_grid.py      |  6 +++
 sharc/topology/topology_imt_mss_dc.py        | 12 +++--
 tests/parameters/parameters_for_testing.yaml |  2 +
 tests/parameters/test_parameters.py          |  4 ++
 tests/test_topology_imt_mss_dc.py            | 46 ++++++++++++++++++++
 5 files changed, 66 insertions(+), 4 deletions(-)

diff --git a/sharc/parameters/imt/parameters_grid.py b/sharc/parameters/imt/parameters_grid.py
index bd6e0144..551a5d2d 100644
--- a/sharc/parameters/imt/parameters_grid.py
+++ b/sharc/parameters/imt/parameters_grid.py
@@ -339,6 +339,9 @@ class ParametersSatelliteWithServiceGrid(ParametersTerrestrialGrid):
 
     beam_radius: float = None
 
+    # [deg]
+    minimum_service_angle: float = 5.0
+
     def validate(self, ctx):
         """Validates instance parameters.
         Parameters
@@ -353,6 +356,9 @@ def validate(self, ctx):
             )
         self.cell_radius = self.beam_radius
 
+        if self.minimum_service_angle < 0. or self.minimum_service_angle > 90:
+            raise ValueError(f"{ctx}.minimum_service_angle should be in [0, 90]")
+
         if not isinstance(
                 self.eligible_sats_margin_from_border,
                 float) and not isinstance(
diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py
index 01debf44..8c22dc7c 100644
--- a/sharc/topology/topology_imt_mss_dc.py
+++ b/sharc/topology/topology_imt_mss_dc.py
@@ -408,7 +408,7 @@ def get_satellite_pointing(
             eligible_sats_msk &= polygon_mask
             eligible_sats_idx = np.where(eligible_sats_msk)[0]
 
-            elev = calc_elevation(
+            all_elevations = calc_elevation(
                 grid_lat[:, np.newaxis],
                 all_sat_lat[eligible_sats_msk][np.newaxis, :],
                 grid_lon[:, np.newaxis],
@@ -417,7 +417,7 @@ def get_satellite_pointing(
                 es_height=0,
             )
 
-            best_sats = elev.argmax(axis=-1)
+            best_sats = all_elevations.argmax(axis=-1)
 
             best_sats_true = eligible_sats_idx[best_sats]
 
@@ -450,8 +450,12 @@ def get_satellite_pointing(
 
             sat_points_towards = defaultdict(list)
 
-            for i, sat in enumerate(best_sats_true):
-                sat_points_towards[sat].append(i)
+            min_service_angle = orbit_params.beam_positioning.service_grid.minimum_service_angle
+            for i, sat in enumerate(best_sats):
+                if all_elevations[i][sat] < min_service_angle:
+                    continue
+                actual_best_sat_idx = eligible_sats_idx[sat]
+                sat_points_towards[actual_best_sat_idx].append(i)
 
             # now only return the angles that
             # the caller asked with the active_sat_idxs parameter
diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml
index 1cd0d9fc..208116b6 100644
--- a/tests/parameters/parameters_for_testing.yaml
+++ b/tests/parameters/parameters_for_testing.yaml
@@ -264,6 +264,7 @@ imt:
                     # if positive, makes border smaller by x km
                     # if negative, makes border bigger by x km
                     eligible_sats_margin_from_border: -2.1
+                    minimum_service_angle: 10.1
             sat_is_active_if:
                 # for a satellite to be active, it needs to respect ALL conditions
                 conditions:
@@ -1096,6 +1097,7 @@ mss_d2d:
             # if positive, makes border smaller by x km
             # if negative, makes border bigger by x km
             eligible_sats_margin_from_border: -2.1
+            minimum_service_angle: 12.1
     sat_is_active_if:
         # for a satellite to be active, it needs to respect ALL conditions
         conditions:
diff --git a/tests/parameters/test_parameters.py b/tests/parameters/test_parameters.py
index e174b53d..f87b5eee 100644
--- a/tests/parameters/test_parameters.py
+++ b/tests/parameters/test_parameters.py
@@ -313,6 +313,8 @@ def test_parameters_imt(self):
             0.11)
         self.assertEqual(
             self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.eligible_sats_margin_from_border, -2.1)
+        self.assertEqual(
+            self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.minimum_service_angle, 10.1)
         self.assertEqual(len(
             self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.grid_in_zone.from_countries.country_names), 2)
         self.assertEqual(
@@ -709,6 +711,8 @@ def test_parametes_mss_d2d(self):
             0.11)
         self.assertEqual(
             self.parameters.mss_d2d.beam_positioning.service_grid.eligible_sats_margin_from_border, -2.1)
+        self.assertEqual(
+            self.parameters.mss_d2d.beam_positioning.service_grid.minimum_service_angle, 12.1)
         self.assertEqual(
             len(self.parameters.mss_d2d.beam_positioning.service_grid.grid_in_zone.from_countries.country_names), 2)
         self.assertEqual(
diff --git a/tests/test_topology_imt_mss_dc.py b/tests/test_topology_imt_mss_dc.py
index 481e7ac8..cfa669d9 100644
--- a/tests/test_topology_imt_mss_dc.py
+++ b/tests/test_topology_imt_mss_dc.py
@@ -6,6 +6,7 @@
 from sharc.station_manager import StationManager
 from sharc.parameters.parameters_orbit import ParametersOrbit
 from sharc.support.sharc_geom import CoordinateSystem, lla2ecef
+from sharc.satellite.utils.sat_utils import calc_elevation
 
 
 class TestTopologyImtMssDc(unittest.TestCase):
@@ -156,6 +157,51 @@ def test_visible_satellites(self):
         # oblateness
         npt.assert_array_less(min_elevation_angle, xy_plane_elevations)
 
+    def test_minimum_service_angle(self):
+        """Test minimum visilibity angle for service grid service."""
+        orbit = ParametersOrbit(
+            n_planes=1,
+            sats_per_plane=1,
+            phasing_deg=3.9,
+            long_asc_deg=18.0,
+            inclination_deg=54.5,
+            perigee_alt_km=525,
+            apogee_alt_km=525,
+        )
+        self.coordinate_system.set_reference(0.0, 0.0, 0)
+
+        self.params.orbits = [orbit]
+        self.params.beam_positioning.type = "SERVICE_GRID"
+        self.params.beam_positioning.service_grid.grid_in_zone.type = "CIRCLE"
+        self.params.beam_positioning.service_grid.grid_in_zone.circle.center_lat = 0.0
+        self.params.beam_positioning.service_grid.grid_in_zone.circle.center_lon = 0.0
+        self.params.beam_positioning.service_grid.grid_in_zone.circle.radius_km = 10 * 111.0
+        self.params.validate("")
+
+        self.imt_mss_dc_topology = TopologyImtMssDc(
+            self.params, self.coordinate_system)
+
+        n_previous_selected = np.inf
+        for a in [5, 50, 80]:
+            self.imt_mss_dc_topology.orbit_params.beam_positioning.service_grid.minimum_service_angle = a
+            self.imt_mss_dc_topology.calculate_coordinates(
+                random_number_gen=np.random.RandomState(8))
+
+            lon_lat_grid = self.params.beam_positioning.service_grid.lon_lat_grid
+            elev_from_bs = calc_elevation(
+                lon_lat_grid[1],
+                self.imt_mss_dc_topology.lat[0],
+                lon_lat_grid[0],
+                self.imt_mss_dc_topology.lon[0],
+                sat_height=self.imt_mss_dc_topology.height[0],
+                es_height=0.0,
+            )
+            n_selected = np.sum(elev_from_bs >= a)
+            self.assertLess(n_selected, n_previous_selected)
+            self.assertLess(n_selected, len(elev_from_bs))
+            self.assertEqual(n_selected, self.imt_mss_dc_topology.num_base_stations)
+            n_previous_selected = n_selected
+
 
 if __name__ == '__main__':
     unittest.main()

From 34f455d31ab1b1eaab6b240a4e356f3cbf5aa13e Mon Sep 17 00:00:00 2001
From: artistrea <apiravesi@gmail.com>
Date: Wed, 5 Nov 2025 15:24:02 -0300
Subject: [PATCH 2/9] fix(topology): imt mss dc topology could create
 num_stations = 0 fixed

---
 sharc/topology/topology_imt_mss_dc.py | 447 +++++++++++++-------------
 1 file changed, 225 insertions(+), 222 deletions(-)

diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py
index 8c22dc7c..e79e1dcb 100644
--- a/sharc/topology/topology_imt_mss_dc.py
+++ b/sharc/topology/topology_imt_mss_dc.py
@@ -75,117 +75,101 @@ def get_coordinates(
     ):
         """Compute the coordinates of the visible space stations."""
         orbit_params.sat_is_active_if.validate("orbit_params.sat_is_active_if")
-        # Calculate the total number of satellites across all orbits
-        total_satellites = sum(
-            orbit.n_planes *
-            orbit.sats_per_plane for orbit in orbit_params.orbits)
-        if any([
-            not hasattr(orbit_params, attr)
-            for attr in ["sat_is_active_if", "orbits", "beam_radius", "num_beams", "beam_positioning"]
-        ]):
-            raise ValueError(
-                "Parameter passed to TopologyImtMssDc needs to contain all of the attributes:\n"
-                '["sat_is_active_if", "orbits", "beam_radius", "num_beams", "beam_positioning"]')
-
-        idx_orbit = np.zeros(
-            total_satellites,
-            dtype=int)  # Add orbit index array
-
-        # List to store indices of active satellites
-        active_satellite_idxs = []
+        num_base_stations = 0
 
         MAX_ITER = 10000  # Maximum iterations to find at least one visible satellite
-        i = 0  # Iteration counter for ensuring satellite visibility
-        while len(active_satellite_idxs) == 0:
-            # Initialize arrays to store satellite positions, angles and
-            # distance from center of earth
-            all_positions = {
-                "R": [],
-                "lat": [],
-                "lon": [],
-                "sx": [],
-                "sy": [],
-                "sz": [],
-                "alt": []}
-            all_elevations = []  # Store satellite elevations
-            all_azimuths = []  # Store satellite azimuths
-
-            current_sat_idx = 0  # Index tracker for satellites across all orbits
-
-            # Iterate through each orbit defined in the parameters
-            for orbit_idx, param in enumerate(orbit_params.orbits):
-                orbit = OrbitModel(
-                    Nsp=param.sats_per_plane,  # Satellites per plane
-                    Np=param.n_planes,  # Number of orbital planes
-                    phasing=param.phasing_deg,  # Phasing angle in degrees
-                    long_asc=param.long_asc_deg,  # Longitude of ascending node in degrees
-                    omega=param.omega_deg,  # Argument of perigee in degrees
-                    delta=param.inclination_deg,  # Orbital inclination in degrees
-                    hp=param.perigee_alt_km,  # Perigee altitude in kilometers
-                    ha=param.apogee_alt_km,  # Apogee altitude in kilometers
-                    Mo=param.initial_mean_anomaly,  # Initial mean anomaly in degrees
-                    # whether to use only time as random variable
-                    model_time_as_random_variable=param.model_time_as_random_variable,
-                    t_min=param.t_min,
-                    t_max=param.t_max,
-                )
-                # Generate random positions for satellites in this orbit
-                pos_vec = orbit.get_orbit_positions_random(
-                    rng=random_number_gen)
-
-                # Determine the number of satellites in this orbit
-                num_satellites = len(pos_vec["sx"])
-
-                # Assign orbit index to satellites
-                idx_orbit[current_sat_idx:current_sat_idx +
-                          num_satellites] = orbit_idx
-
-                # Extract satellite positions and calculate distances
-                sx, sy, sz = pos_vec['sx'], pos_vec['sy'], pos_vec['sz']
-                # Distance from Earth's center
-                r = np.sqrt(sx**2 + sy**2 + sz**2)
-
-                # When getting azimuth and elevation, we need to consider sx, sy and sz points
-                # from the center of earth to the satellite, and we need to point the satellite
-                # towards the center of earth
-                # Calculate elevation angles
-                elevations = np.degrees(np.arcsin(-sz / r))
-                # Calculate azimuth angles
-                azimuths = np.degrees(np.arctan2(-sy, -sx))
-
-                # Append satellite positions and angles to global lists
-                all_positions['lat'].extend(pos_vec['lat'])  # Latitudes
-                all_positions['lon'].extend(pos_vec['lon'])  # Longitudes
-                all_positions['sx'].extend(sx)  # X-coordinates
-                all_positions['sy'].extend(sy)  # Y-coordinates
-                all_positions['sz'].extend(sz)  # Z-coordinates
-                all_positions["R"].extend(r)
-                all_positions["alt"].extend(pos_vec['alt'])
-                all_elevations.extend(elevations)  # Elevation angles
-                all_azimuths.extend(azimuths)  # Azimuth angles
-
-                active_sats_mask = np.ones(len(pos_vec['lat']), dtype=bool)
-
-                if "MINIMUM_ELEVATION_FROM_ES" in orbit_params.sat_is_active_if.conditions:
-                    # Calculate satellite visibility from base stations
-                    elev_from_bs = calc_elevation(
-                        coordinate_system.ref_lat,  # Latitude of base station
-                        pos_vec['lat'],  # Latitude of satellites
-                        coordinate_system.ref_long,  # Longitude of base station
-                        pos_vec['lon'],  # Longitude of satellites
-                        # Perigee altitude in kilometers
-                        sat_height=pos_vec['alt'] * 1e3,
-                        es_height=coordinate_system.ref_alt,
+        while num_base_stations == 0:
+            # Calculate the total number of satellites across all orbits
+            total_satellites = sum(
+                orbit.n_planes *
+                orbit.sats_per_plane for orbit in orbit_params.orbits)
+            if any([
+                not hasattr(orbit_params, attr)
+                for attr in ["sat_is_active_if", "orbits", "beam_radius", "num_beams", "beam_positioning"]
+            ]):
+                raise ValueError(
+                    "Parameter passed to TopologyImtMssDc needs to contain all of the attributes:\n"
+                    '["sat_is_active_if", "orbits", "beam_radius", "num_beams", "beam_positioning"]')
+
+            idx_orbit = np.zeros(
+                total_satellites,
+                dtype=int)  # Add orbit index array
+
+            # List to store indices of active satellites
+            active_satellite_idxs = []
+
+            i = 0  # Iteration counter for ensuring satellite visibility
+            while len(active_satellite_idxs) == 0:
+                # Initialize arrays to store satellite positions, angles and
+                # distance from center of earth
+                all_positions = {
+                    "R": [],
+                    "lat": [],
+                    "lon": [],
+                    "sx": [],
+                    "sy": [],
+                    "sz": [],
+                    "alt": []}
+                all_elevations = []  # Store satellite elevations
+                all_azimuths = []  # Store satellite azimuths
+
+                current_sat_idx = 0  # Index tracker for satellites across all orbits
+
+                # Iterate through each orbit defined in the parameters
+                for orbit_idx, param in enumerate(orbit_params.orbits):
+                    orbit = OrbitModel(
+                        Nsp=param.sats_per_plane,  # Satellites per plane
+                        Np=param.n_planes,  # Number of orbital planes
+                        phasing=param.phasing_deg,  # Phasing angle in degrees
+                        long_asc=param.long_asc_deg,  # Longitude of ascending node in degrees
+                        omega=param.omega_deg,  # Argument of perigee in degrees
+                        delta=param.inclination_deg,  # Orbital inclination in degrees
+                        hp=param.perigee_alt_km,  # Perigee altitude in kilometers
+                        ha=param.apogee_alt_km,  # Apogee altitude in kilometers
+                        Mo=param.initial_mean_anomaly,  # Initial mean anomaly in degrees
+                        # whether to use only time as random variable
+                        model_time_as_random_variable=param.model_time_as_random_variable,
+                        t_min=param.t_min,
+                        t_max=param.t_max,
                     )
-
-                    # Determine visible satellites based on minimum elevation
-                    # angle
-                    active_sats_mask = active_sats_mask & (elev_from_bs.flatten(
-                    ) >= orbit_params.sat_is_active_if.minimum_elevation_from_es)
-
-                if "MAXIMUM_ELEVATION_FROM_ES" in orbit_params.sat_is_active_if.conditions:
-                    # no need to recalculate if already calculated above
-                    if "MINIMUM_ELEVATION_FROM_ES" not in orbit_params.sat_is_active_if.conditions:
+                    # Generate random positions for satellites in this orbit
+                    pos_vec = orbit.get_orbit_positions_random(
+                        rng=random_number_gen)
+
+                    # Determine the number of satellites in this orbit
+                    num_satellites = len(pos_vec["sx"])
+
+                    # Assign orbit index to satellites
+                    idx_orbit[current_sat_idx:current_sat_idx +
+                              num_satellites] = orbit_idx
+
+                    # Extract satellite positions and calculate distances
+                    sx, sy, sz = pos_vec['sx'], pos_vec['sy'], pos_vec['sz']
+                    # Distance from Earth's center
+                    r = np.sqrt(sx**2 + sy**2 + sz**2)
+
+                    # When getting azimuth and elevation, we need to consider sx, sy and sz points
+                    # from the center of earth to the satellite, and we need to point the satellite
+                    # towards the center of earth
+                    # Calculate elevation angles
+                    elevations = np.degrees(np.arcsin(-sz / r))
+                    # Calculate azimuth angles
+                    azimuths = np.degrees(np.arctan2(-sy, -sx))
+
+                    # Append satellite positions and angles to global lists
+                    all_positions['lat'].extend(pos_vec['lat'])  # Latitudes
+                    all_positions['lon'].extend(pos_vec['lon'])  # Longitudes
+                    all_positions['sx'].extend(sx)  # X-coordinates
+                    all_positions['sy'].extend(sy)  # Y-coordinates
+                    all_positions['sz'].extend(sz)  # Z-coordinates
+                    all_positions["R"].extend(r)
+                    all_positions["alt"].extend(pos_vec['alt'])
+                    all_elevations.extend(elevations)  # Elevation angles
+                    all_azimuths.extend(azimuths)  # Azimuth angles
+
+                    active_sats_mask = np.ones(len(pos_vec['lat']), dtype=bool)
+
+                    if "MINIMUM_ELEVATION_FROM_ES" in orbit_params.sat_is_active_if.conditions:
                         # Calculate satellite visibility from base stations
                         elev_from_bs = calc_elevation(
                             coordinate_system.ref_lat,  # Latitude of base station
@@ -197,122 +181,141 @@ def get_coordinates(
                             es_height=coordinate_system.ref_alt,
                         )
 
-                    # Determine visible satellites based on minimum elevation
-                    # angle
-                    active_sats_mask = active_sats_mask & (elev_from_bs.flatten(
-                    ) <= orbit_params.sat_is_active_if.maximum_elevation_from_es)
-
-                if "LAT_LONG_INSIDE_COUNTRY" in orbit_params.sat_is_active_if.conditions:
-                    flat_active_lon = pos_vec["lon"].flatten()[
-                        active_sats_mask]
-                    flat_active_lat = pos_vec["lat"].flatten()[
-                        active_sats_mask]
-
-                    # create points(lon, lat) to compare to country
-                    sats_points = gpd.points_from_xy(
-                        flat_active_lon, flat_active_lat, crs=EARTH_DEFAULT_CRS)
-
-                    # Check if the satellite is inside the country polygon
-                    polygon_mask = np.zeros_like(active_sats_mask)
-                    polygon_mask[active_sats_mask] = sats_points.within(
-                        orbit_params.sat_is_active_if.lat_long_inside_country.filter_polygon)
-
-                    active_sats_mask = active_sats_mask & polygon_mask
-
-                visible_sat_idxs = np.arange(
-                    current_sat_idx, current_sat_idx + len(pos_vec['lat']), dtype=int
-                )[active_sats_mask]
-                active_satellite_idxs.extend(visible_sat_idxs)
-
-                # Update the index tracker for the next orbit
-                current_sat_idx += len(sx)
-
-            i += 1  # Increment iteration counter
-            if i >= MAX_ITER:  # Check if maximum iterations reached
-                raise RuntimeError(
-                    "Maximum iterations reached, and no satellite was selected within the minimum elevation criteria."
-                )
-        # We have the list of visible satellites, now create a Topolgy of this subset and move the coordinate system
-        # reference.
-        # Convert X-coordinates to meters
-        all_space_station_x = np.ravel(np.array(all_positions['sx'])) * 1e3
-        # Convert Y-coordinates to meters
-        all_space_station_y = np.ravel(np.array(all_positions['sy'])) * 1e3
-        # Convert Z-coordinates to meters
-        all_space_station_z = np.ravel(np.array(all_positions['sz'])) * 1e3
-        all_elevation = np.ravel(np.array(all_elevations))  # Elevation angles
-        all_azimuth = np.ravel(np.array(all_azimuths))  # Azimuth angles
-        all_lat = np.ravel(np.array(all_positions['lat']))
-        all_lon = np.ravel(np.array(all_positions['lon']))
-        all_sat_altitude = np.ravel(np.array(all_positions['alt'])) * 1e3
-
-        total_active_satellites = len(active_satellite_idxs)
-        space_station_x = all_space_station_x[active_satellite_idxs]
-        space_station_y = all_space_station_y[active_satellite_idxs]
-        space_station_z = all_space_station_z[active_satellite_idxs]
-        elevation = all_elevation[active_satellite_idxs]
-        azimuth = all_azimuth[active_satellite_idxs]
-        lat = all_lat[active_satellite_idxs]
-        lon = all_lon[active_satellite_idxs]
-        sat_altitude = all_sat_altitude[active_satellite_idxs]
-        sat_altitude = all_sat_altitude[active_satellite_idxs]
+                        # Determine visible satellites based on minimum elevation
+                        # angle
+                        active_sats_mask = active_sats_mask & (elev_from_bs.flatten(
+                        ) >= orbit_params.sat_is_active_if.minimum_elevation_from_es)
+
+                    if "MAXIMUM_ELEVATION_FROM_ES" in orbit_params.sat_is_active_if.conditions:
+                        # no need to recalculate if already calculated above
+                        if "MINIMUM_ELEVATION_FROM_ES" not in orbit_params.sat_is_active_if.conditions:
+                            # Calculate satellite visibility from base stations
+                            elev_from_bs = calc_elevation(
+                                coordinate_system.ref_lat,  # Latitude of base station
+                                pos_vec['lat'],  # Latitude of satellites
+                                coordinate_system.ref_long,  # Longitude of base station
+                                pos_vec['lon'],  # Longitude of satellites
+                                # Perigee altitude in kilometers
+                                sat_height=pos_vec['alt'] * 1e3,
+                                es_height=coordinate_system.ref_alt,
+                            )
+
+                        # Determine visible satellites based on minimum elevation
+                        # angle
+                        active_sats_mask = active_sats_mask & (elev_from_bs.flatten(
+                        ) <= orbit_params.sat_is_active_if.maximum_elevation_from_es)
+
+                    if "LAT_LONG_INSIDE_COUNTRY" in orbit_params.sat_is_active_if.conditions:
+                        flat_active_lon = pos_vec["lon"].flatten()[
+                            active_sats_mask]
+                        flat_active_lat = pos_vec["lat"].flatten()[
+                            active_sats_mask]
+
+                        # create points(lon, lat) to compare to country
+                        sats_points = gpd.points_from_xy(
+                            flat_active_lon, flat_active_lat, crs=EARTH_DEFAULT_CRS)
+
+                        # Check if the satellite is inside the country polygon
+                        polygon_mask = np.zeros_like(active_sats_mask)
+                        polygon_mask[active_sats_mask] = sats_points.within(
+                            orbit_params.sat_is_active_if.lat_long_inside_country.filter_polygon)
+
+                        active_sats_mask = active_sats_mask & polygon_mask
+
+                    visible_sat_idxs = np.arange(
+                        current_sat_idx, current_sat_idx + len(pos_vec['lat']), dtype=int
+                    )[active_sats_mask]
+                    active_satellite_idxs.extend(visible_sat_idxs)
+
+                    # Update the index tracker for the next orbit
+                    current_sat_idx += len(sx)
+
+                i += 1  # Increment iteration counter
+                if i >= MAX_ITER:  # Check if maximum iterations reached
+                    raise RuntimeError(
+                        "Maximum iterations reached, and no satellite was selected within the minimum elevation criteria."
+                    )
+            # We have the list of visible satellites, now create a Topolgy of this subset and move the coordinate system
+            # reference.
+            # Convert X-coordinates to meters
+            all_space_station_x = np.ravel(np.array(all_positions['sx'])) * 1e3
+            # Convert Y-coordinates to meters
+            all_space_station_y = np.ravel(np.array(all_positions['sy'])) * 1e3
+            # Convert Z-coordinates to meters
+            all_space_station_z = np.ravel(np.array(all_positions['sz'])) * 1e3
+            all_elevation = np.ravel(np.array(all_elevations))  # Elevation angles
+            all_azimuth = np.ravel(np.array(all_azimuths))  # Azimuth angles
+            all_lat = np.ravel(np.array(all_positions['lat']))
+            all_lon = np.ravel(np.array(all_positions['lon']))
+            all_sat_altitude = np.ravel(np.array(all_positions['alt'])) * 1e3
+
+            total_active_satellites = len(active_satellite_idxs)
+            space_station_x = all_space_station_x[active_satellite_idxs]
+            space_station_y = all_space_station_y[active_satellite_idxs]
+            space_station_z = all_space_station_z[active_satellite_idxs]
+            elevation = all_elevation[active_satellite_idxs]
+            azimuth = all_azimuth[active_satellite_idxs]
+            lat = all_lat[active_satellite_idxs]
+            lon = all_lon[active_satellite_idxs]
+            sat_altitude = all_sat_altitude[active_satellite_idxs]
+            sat_altitude = all_sat_altitude[active_satellite_idxs]
+
+            # Convert the ECEF coordinates to the transformed cartesian coordinates and set the Space Station positions
+            # used to generetate the IMT Base Stations
+            space_station_x, space_station_y, space_station_z = \
+                coordinate_system.ecef2enu(space_station_x, space_station_y, space_station_z)
+
+            # Rotate the azimuth and elevation angles off the center beam the new
+            # transformed cartesian coordinates
+            r = 1
+            # transform pointing vectors, without considering geodesical earth
+            # coord system
+            pointing_vec_x, pointing_vec_y, pointing_vec_z = polar_to_cartesian(
+                r, all_azimuth, all_elevation)
+            pointing_vec_x, pointing_vec_y, pointing_vec_z = \
+                coordinate_system.ecef2enu(
+                    pointing_vec_x, pointing_vec_y, pointing_vec_z, translate=0)
+            _, all_azimuth, all_elevation = cartesian_to_polar(
+                pointing_vec_x, pointing_vec_y, pointing_vec_z)
 
-        # Convert the ECEF coordinates to the transformed cartesian coordinates and set the Space Station positions
-        # used to generetate the IMT Base Stations
-        space_station_x, space_station_y, space_station_z = \
-            coordinate_system.ecef2enu(space_station_x, space_station_y, space_station_z)
-
-        # Rotate the azimuth and elevation angles off the center beam the new
-        # transformed cartesian coordinates
-        r = 1
-        # transform pointing vectors, without considering geodesical earth
-        # coord system
-        pointing_vec_x, pointing_vec_y, pointing_vec_z = polar_to_cartesian(
-            r, all_azimuth, all_elevation)
-        pointing_vec_x, pointing_vec_y, pointing_vec_z = \
-            coordinate_system.ecef2enu(
-                pointing_vec_x, pointing_vec_y, pointing_vec_z, translate=0)
-        _, all_azimuth, all_elevation = cartesian_to_polar(
-            pointing_vec_x, pointing_vec_y, pointing_vec_z)
-
-        beams_elev, beams_azim, sx, sy = TopologyImtMssDc.get_satellite_pointing(
-            random_number_gen,
-            coordinate_system,
-            orbit_params,
-            total_active_satellites,
-            all_space_station_x, all_space_station_y, all_space_station_z,
-            all_azimuth,
-            all_elevation,
-            all_lat, all_lon, all_sat_altitude,
-            active_satellite_idxs
-        )
+            beams_elev, beams_azim, sx, sy = TopologyImtMssDc.get_satellite_pointing(
+                random_number_gen,
+                coordinate_system,
+                orbit_params,
+                total_active_satellites,
+                all_space_station_x, all_space_station_y, all_space_station_z,
+                all_azimuth,
+                all_elevation,
+                all_lat, all_lon, all_sat_altitude,
+                active_satellite_idxs
+            )
 
-        # In SHARC each sector is treated as a separate base station, so we need to repeat the satellite positions
-        # for each sector.
-        sat_ocurr = [len(x) for x in beams_elev]
-
-        elevation = np.array(
-            functools.reduce(
-                lambda x,
-                y: list(x) +
-                list(y),
-                beams_elev))
-        azimuth = np.array(
-            functools.reduce(
-                lambda x,
-                y: list(x) +
-                list(y),
-                beams_azim))
-
-        space_station_x = np.repeat(space_station_x, sat_ocurr)
-        space_station_y = np.repeat(space_station_y, sat_ocurr)
-        space_station_z = np.repeat(space_station_z, sat_ocurr)
-
-        num_base_stations = np.sum(sat_ocurr)
-        lat = np.repeat(lat, sat_ocurr)
-        lon = np.repeat(lon, sat_ocurr)
-
-        altitudes = np.repeat(sat_altitude, sat_ocurr)
+            # In SHARC each sector is treated as a separate base station, so we need to repeat the satellite positions
+            # for each sector.
+            sat_ocurr = [len(x) for x in beams_elev]
+
+            elevation = np.array(
+                functools.reduce(
+                    lambda x,
+                    y: list(x) +
+                    list(y),
+                    beams_elev))
+            azimuth = np.array(
+                functools.reduce(
+                    lambda x,
+                    y: list(x) +
+                    list(y),
+                    beams_azim))
+
+            space_station_x = np.repeat(space_station_x, sat_ocurr)
+            space_station_y = np.repeat(space_station_y, sat_ocurr)
+            space_station_z = np.repeat(space_station_z, sat_ocurr)
+
+            num_base_stations = np.sum(sat_ocurr)
+            lat = np.repeat(lat, sat_ocurr)
+            lon = np.repeat(lon, sat_ocurr)
+
+            altitudes = np.repeat(sat_altitude, sat_ocurr)
 
         assert (space_station_x.shape == (num_base_stations,))
         assert (space_station_y.shape == (num_base_stations,))

From 9cd3acd8ef11e97b754fe663a3c61a4817f67ec7 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Mon, 10 Nov 2025 17:29:46 -0300
Subject: [PATCH 3/9] update(parameters): Added z and Lf as parameters for
 ITU-R S.1528 Rec 1.2 antenna pattern.

---
 sharc/antenna/antenna_s1528.py                | 41 +++++++++++--------
 .../antenna/parameters_antenna_s1528.py       |  8 ++++
 tests/parameters/parameters_for_testing.yaml  |  3 ++
 tests/parameters/test_parameters.py           | 19 +++++++++
 tests/test_antenna_s1528.py                   |  2 +
 5 files changed, 57 insertions(+), 16 deletions(-)

diff --git a/sharc/antenna/antenna_s1528.py b/sharc/antenna/antenna_s1528.py
index b8ee9076..de9a966b 100644
--- a/sharc/antenna/antenna_s1528.py
+++ b/sharc/antenna/antenna_s1528.py
@@ -175,12 +175,18 @@ def __init__(self, param: ParametersAntennaS1528):
         # the system design
         self.l_s = param.antenna_l_s
 
+        # far-out side-lobe level [dBi]
+        if param.far_out_side_lobe is None:
+            self.l_f = 0
+        else:
+            self.l_f = param.far_out_side_lobe
+
         # for elliptical antennas, this is the ratio major axis/minor axis
         # we assume circular antennas, so z = 1
-        self.z = 1
-
-        # far-out side-lobe level [dBi]
-        self.l_f = 0
+        if param.major_minor_axis_ratio is None:
+            self.z = 1
+        else:
+            self.z = param.major_minor_axis_ratio
 
         # back-lobe level
         self.l_b = np.maximum(
@@ -190,19 +196,22 @@ def __init__(self, param: ParametersAntennaS1528):
         # one-half the 3 dB beamwidth in the plane of interest
         self.psi_b = param.antenna_3_dB_bw / 2
 
-        if self.l_s == -15:
-            self.a = 2.58 * math.sqrt(1 - 1.4 * math.log10(self.z))
-        elif self.l_s == -20:
-            self.a = 2.58 * math.sqrt(1 - 1.0 * math.log10(self.z))
-        elif self.l_s == -25:
-            self.a = 2.58 * math.sqrt(1 - 0.6 * math.log10(self.z))
-        elif self.l_s == -30:
-            self.a = 2.58 * math.sqrt(1 - 0.4 * math.log10(self.z))
+        if math.isclose(self.z, 1.0):
+            self.a = 2.58
         else:
-            sys.stderr.write(
-                "ERROR\nInvalid AntennaS1528 L_s parameter: " + str(self.l_s),
-            )
-            sys.exit(1)
+            if self.l_s == -15:
+                self.a = 2.58 * math.sqrt(1 - 1.4 * math.log10(self.z))
+            elif self.l_s == -20:
+                self.a = 2.58 * math.sqrt(1 - 1.0 * math.log10(self.z))
+            elif self.l_s == -25:
+                self.a = 2.58 * math.sqrt(1 - 0.6 * math.log10(self.z))
+            elif self.l_s == -30:
+                self.a = 2.58 * math.sqrt(1 - 0.4 * math.log10(self.z))
+            else:
+                sys.stderr.write(
+                    f"ERROR\nInvalid AntennaS1528 L_s parameter {self.l_s} for z={self.z}"
+                )
+                sys.exit(1)
 
         self.b = 6.32
         self.alpha = 1.5
diff --git a/sharc/parameters/antenna/parameters_antenna_s1528.py b/sharc/parameters/antenna/parameters_antenna_s1528.py
index 4cd5ea35..ed2f204f 100644
--- a/sharc/parameters/antenna/parameters_antenna_s1528.py
+++ b/sharc/parameters/antenna/parameters_antenna_s1528.py
@@ -37,6 +37,14 @@ class ParametersAntennaS1528(ParametersBase):
     l_r: float = None
     l_t: float = None
 
+    # Recommends 1.2 only
+    # For elliptical antennas, this is the ratio major axis/minor axis
+    # we assume circular antennas, so z = 1
+    major_minor_axis_ratio: float = None
+
+    # Far-out side-lobe level
+    far_out_side_lobe: float = None
+
     def load_parameters_from_file(self, config_file: str):
         """Load the parameters from file an run a sanity check.
 
diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml
index 208116b6..86631172 100644
--- a/tests/parameters/parameters_for_testing.yaml
+++ b/tests/parameters/parameters_for_testing.yaml
@@ -999,6 +999,9 @@ single_space_station:
             # Antenna diameter [m]
             # if no diameter is passed, a diameter will be assumed according to document
             diameter: 2.12
+        itu_r_s_1528:
+            major_minor_axis_ratio: 1.0
+            far_out_side_lobe: -25
     ###########################################################################
     # Selected channel model
     channel_model: "P619"
diff --git a/tests/parameters/test_parameters.py b/tests/parameters/test_parameters.py
index f87b5eee..b959acfb 100644
--- a/tests/parameters/test_parameters.py
+++ b/tests/parameters/test_parameters.py
@@ -875,6 +875,25 @@ def test_parameters_single_space_station(self):
             self.parameters.single_space_station.antenna.gain,
         )
 
+        self.assertEqual(
+            self.parameters.single_space_station.antenna.itu_r_s_1528.antenna_gain,
+            self.parameters.single_space_station.antenna.gain,
+        )
+        self.assertEqual(
+            self.parameters.single_space_station.antenna.itu_r_s_1528.major_minor_axis_ratio,
+            1.0
+        )
+
+        self.assertEqual(
+            self.parameters.single_space_station.antenna.itu_r_s_1528.major_minor_axis_ratio,
+            1.0
+        )
+
+        self.assertEqual(
+            self.parameters.single_space_station.antenna.itu_r_s_1528.far_out_side_lobe,
+            -25
+        )
+
         self.assertEqual(
             self.parameters.single_space_station.param_p619.earth_station_alt_m,
             self.parameters.single_space_station.geometry.es_altitude,
diff --git a/tests/test_antenna_s1528.py b/tests/test_antenna_s1528.py
index 0d2aae6c..b025506b 100644
--- a/tests/test_antenna_s1528.py
+++ b/tests/test_antenna_s1528.py
@@ -23,6 +23,8 @@ def setUp(self):
         param.antenna_gain = 39
         param.antenna_pattern = "ITU-R S.1528-0"
         param.antenna_3_dB_bw = 2
+        param.major_minor_axis_ratio = 1
+        param.far_out_side_lobe = 0
 
         param.antenna_l_s = -20
         self.antenna20 = AntennaS1528(param)

From 741f854ea0bf43a62988639be7ce77b0313c8a34 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Tue, 2 Dec 2025 13:26:07 -0300
Subject: [PATCH 4/9] fix(propagation): Fixed the season parameter check in
 Atmosfere P.835 to

---
 sharc/propagation/atmosphere.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/sharc/propagation/atmosphere.py b/sharc/propagation/atmosphere.py
index 846ce151..d4d382fb 100644
--- a/sharc/propagation/atmosphere.py
+++ b/sharc/propagation/atmosphere.py
@@ -385,6 +385,9 @@ def get_reference_atmosphere_p835(
         """
 
         h_km = altitude / 1000
+        season = season.lower()
+        if season not in ['winter', 'summer']:
+            raise ValueError(f"Invalid season name f{season}.")
 
         if latitude <= 22:
             # low latitude

From 595bddff4de2567190c8ff6f0b751eda147479bc Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Thu, 4 Dec 2025 10:52:40 -0300
Subject: [PATCH 5/9] fix(simulation): Fixed simulator log facility setup *
 Also added a CLI option to stream simulation log to a file

---
 sharc/main_cli.py             | 35 +++++++++-----------
 sharc/support/logging.yaml    | 27 ----------------
 sharc/support/sharc_logger.py | 61 ++++++++++++++++++++---------------
 3 files changed, 51 insertions(+), 72 deletions(-)
 delete mode 100644 sharc/support/logging.yaml

diff --git a/sharc/main_cli.py b/sharc/main_cli.py
index 6d13173e..a9ffa536 100644
--- a/sharc/main_cli.py
+++ b/sharc/main_cli.py
@@ -7,9 +7,9 @@
 
 import os
 import sys
-import getopt
+import argparse
 
-from sharc.support.sharc_logger import Logging, SimulationLogger
+from sharc.support.sharc_logger import setup_logging, SimulationLogger
 from sharc.controller import Controller
 from sharc.gui.view_cli import ViewCli
 from sharc.model import Model
@@ -34,28 +34,25 @@ def main(argv):
 
     param_file = ""
 
-    try:
-        opts, _ = getopt.getopt(argv, "hp:")
-    except getopt.GetoptError:
-        print("usage: main_cli.py -p <param_file>")
-        sys.exit(2)
-
-    if not opts:
-        param_file = os.path.join(os.getcwd(), "input", "parameters.yaml")
-    else:
-        for opt, arg in opts:
-            if opt == "-h":
-                print("usage: main_cli.py -p <param_file>")
-                sys.exit()
-            elif opt == "-p":
-                param_file = os.path.join(os.getcwd(), arg)
+    parser = argparse.ArgumentParser(description="SHARC - Radio Sharing and Compatiblity Monte Carlo Simulator")
+    parser.add_argument("-p", "--param-file", default=os.path.join(os.getcwd(), "input", "parameters.yaml"),
+                        help="Path to parameter file (default: input/parameters.yaml)")
+    parser.add_argument("-l", "--log-file", default=None,
+                        help="Path to output log file (optional)")
+
+    args = parser.parse_args(argv)
+    param_file = os.path.join(os.getcwd(), args.param_file) if not os.path.isabs(args.param_file) else args.param_file
+
+    log_file = None
+    if args.log_file is not None:
+        log_file = os.path.join(os.getcwd(), args.log_file) if not os.path.isabs(args.log_file) else args.log_file
+
+    setup_logging(log_file=log_file)
 
     # Logger setup start
     sim_logger = SimulationLogger(param_file)
     sim_logger.start()
 
-    Logging.setup_logging()
-
     model = Model()
     view_cli = ViewCli()
     controller = Controller()
diff --git a/sharc/support/logging.yaml b/sharc/support/logging.yaml
deleted file mode 100644
index 0e04fc41..00000000
--- a/sharc/support/logging.yaml
+++ /dev/null
@@ -1,27 +0,0 @@
----
-version: 1
-disable_existing_loggers: False
-formatters:
-    simple:
-        format: "%(message)s"
-#        format: "%(asctime)s - %(name)s - %(message)s"
-
-handlers:
-    console:
-        class: logging.StreamHandler
-        level: INFO
-        formatter: simple
-        stream: ext://sys.stdout
-
-    file:
-        class: logging.handlers.RotatingFileHandler
-        formatter: simple
-        filename: output/logfile.log
-        maxBytes: 10485760 # 10MB
-        backupCount: 20
-        encoding: utf8
-        
-
-root:
-    level: INFO
-    handlers: [console, file]
diff --git a/sharc/support/sharc_logger.py b/sharc/support/sharc_logger.py
index c4771ec5..5d7cc5ac 100644
--- a/sharc/support/sharc_logger.py
+++ b/sharc/support/sharc_logger.py
@@ -2,33 +2,42 @@
 import os
 import sys
 import yaml
-import logging
+import logging.config
 import subprocess
 from pathlib import Path
 from datetime import datetime
-from typing import Optional
-
-
-class Logging:
-    """Logging utility class for configuring application logging."""
-
-    @staticmethod
-    def setup_logging(
-        default_path="support/logging.yaml",
-        default_level=logging.INFO,
-        env_key="LOG_CFG",
-    ):
-        """Set up logging configuration for the application."""
-        path = default_path
-        value = os.getenv(env_key, None)
-        if value:
-            path = value
-        if os.path.exists(path):
-            with open(path, "rt") as f:
-                config = yaml.safe_load(f.read())
-            logging.config.dictConfig(config)
-        else:
-            logging.basicConfig(level=default_level)
+from typing import Optional, List
+
+level_mapping = logging.getLevelNamesMapping()
+
+
+def setup_logging(log_file=None, default_level="INFO"):
+    """Setup logging configuration for the root logger.
+
+    Run this function in the beginning of the simulation to setup the root logger.
+    """
+
+    try:
+        level = level_mapping[default_level]
+    except KeyError:
+        raise ValueError("Invalid log level option {}".format(default_level))
+
+    root_logger = logging.getLogger()
+    root_logger.setLevel(level)
+    formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
+    root_logger.handlers = []
+
+    # Stream to stdout
+    handler = logging.StreamHandler(sys.stdout)
+    handler.setFormatter(formatter)
+    root_logger.addHandler(handler)
+
+    # Stream to file if specified
+    if log_file is not None:
+        file_handler = logging.FileHandler(log_file)
+        file_handler.setFormatter(formatter)
+        root_logger.addHandler(file_handler)
+
 
 
 class SimulationLogger:
@@ -102,7 +111,7 @@ def _find_root_dir(self, folder_name: str) -> Optional[Path]:
                 return parent
         return None
 
-    def _run_git_cmd(self, args: list[str]) -> Optional[str]:
+    def _run_git_cmd(self, args: List[str]) -> Optional[str]:
         try:
             return (
                 subprocess.check_output(["git"] + args, stderr=subprocess.DEVNULL)
@@ -132,7 +141,7 @@ def _get_invocation_command(self) -> str:
     def _get_python_version(self) -> str:
         return sys.version.replace("\n", " ")
 
-    def _get_installed_packages(self) -> list[str]:
+    def _get_installed_packages(self) -> List[str]:
         try:
             output = subprocess.check_output(
                 [sys.executable, "-m", "pip", "freeze"], stderr=subprocess.DEVNULL

From 4cea239eaab1eef555c916c3ed2cd8cfd26a7e96 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Wed, 17 Dec 2025 13:46:58 -0300
Subject: [PATCH 6/9] Adding a modification that uses 1528 antenna + elevation
 angle

---
 sharc/station_factory.py              | 18 ++++++++++++++++++
 sharc/topology/topology_imt_mss_dc.py | 17 ++++++++++++++---
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/sharc/station_factory.py b/sharc/station_factory.py
index 02384d4f..1f25d035 100644
--- a/sharc/station_factory.py
+++ b/sharc/station_factory.py
@@ -1752,6 +1752,7 @@ def generate_mss_d2d(
         z = mss_d2d_values["sat_z"]
         elev = mss_d2d_values["sat_antenna_elev"]
         azim = mss_d2d_values["sat_antenna_azim"]
+        beams_ground_elev = mss_d2d_values["beams_ground_elev"]
         mss_d2d.geom.set_global_coords(
             x, y, z,
             azim, elev,
@@ -1788,6 +1789,23 @@ def generate_mss_d2d(
                 f"generate_mss_ss: Invalid antenna type: {params.antenna.pattern}")
 
         for i in range(mss_d2d.num_stations):
+            if params.antenna.pattern == "Satellite Beamforming":
+                if beams_ground_elev[i] >= 50:
+                    params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
+                    params.antenna.itu_r_s_1528.antenna_gain = 34.6
+                    params.antenna.itu_r_s_1528.antenna_3_dB_bw = 3.06
+                    params.antenna.itu_r_s_1528.antenna_l_s = -35.0
+                    params.antenna.itu_r_s_1528.far_out_side_lobe = -25
+
+                else:
+                    params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
+                    params.antenna.itu_r_s_1528.antenna_gain = 40.0
+                    params.antenna.itu_r_s_1528.antenna_3_dB_bw = 1.68
+                    params.antenna.itu_r_s_1528.antenna_l_s = -20.0
+                    params.antenna.itu_r_s_1528.far_out_side_lobe = -15
+
+                antenna_pattern = AntennaS1528(params.antenna.itu_r_s_1528)
+
             mss_d2d.antenna[i] = antenna_pattern
 
         return mss_d2d  # Return the configured StationManager
diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py
index 337cbe43..88917282 100644
--- a/sharc/topology/topology_imt_mss_dc.py
+++ b/sharc/topology/topology_imt_mss_dc.py
@@ -278,7 +278,7 @@ def get_coordinates(
             _, all_azimuth, all_elevation = cartesian_to_polar(
                 pointing_vec_x, pointing_vec_y, pointing_vec_z)
 
-            beams_elev, beams_azim, sx, sy = TopologyImtMssDc.get_satellite_pointing(
+            beams_elev, beams_azim, beams_ground_elev, sx, sy = TopologyImtMssDc.get_satellite_pointing(
                 random_number_gen,
                 coordinate_system,
                 orbit_params,
@@ -306,6 +306,12 @@ def get_coordinates(
                     y: list(x) +
                     list(y),
                     beams_azim))
+            ground_elev = np.array(
+                functools.reduce(
+                    lambda x,
+                    y: list(x) +
+                    list(y),
+                    beams_ground_elev))
 
             space_station_x = np.repeat(space_station_x, sat_ocurr)
             space_station_y = np.repeat(space_station_y, sat_ocurr)
@@ -349,6 +355,7 @@ def get_coordinates(
             "sat_alt": altitudes,
             "sat_antenna_elev": elevation,
             "sat_antenna_azim": azimuth,
+            "beams_ground_elev": ground_elev,
             "sectors_x": sx,
             "sectors_y": sy,
             "sectors_z": np.zeros_like(sx)
@@ -464,6 +471,7 @@ def get_satellite_pointing(
             # the caller asked with the active_sat_idxs parameter
             beams_azim = []
             beams_elev = []
+            beams_ground_elev = []
             n = 0
 
             for act_sat in active_satellite_idxs:
@@ -471,6 +479,8 @@ def get_satellite_pointing(
                     n += len(sat_points_towards[act_sat])
                     beams_azim.append(azim[sat_points_towards[act_sat]])
                     beams_elev.append(elev[sat_points_towards[act_sat]])
+                    beams_ground_elev.append(all_elevations[
+                        sat_points_towards[act_sat], eligible_sats_idx == act_sat])
                 else:
                     beams_azim.append([])
                     beams_elev.append([])
@@ -480,7 +490,7 @@ def get_satellite_pointing(
             sx = np.zeros(n)
             sy = np.zeros(n)
 
-            return beams_elev, beams_azim, sx, sy
+            return beams_elev, beams_azim, beams_ground_elev, sx, sy
         # We borrow the TopologyNTN method to calculate the sectors azimuth and elevation angles from their
         # respective x and y boresight coordinates
         sx, sy = TopologyNTN.get_sectors_xy(
@@ -572,6 +582,7 @@ def get_satellite_pointing(
                 np.repeat(
                     sat_altitude,
                     orbit_params.num_beams))) - 90
+        beams_ground_elev = 90 - beams_elev
 
         beams_azim = beams_azim.reshape(
             (total_active_satellites, orbit_params.num_beams)
@@ -596,7 +607,7 @@ def get_satellite_pointing(
                 nadir_azim[i]
             )
 
-        return beams_elev, beams_azim, sx, sy
+        return beams_elev, beams_azim, beams_ground_elev, sx, sy
 
     @staticmethod
     def get_distr(

From a0f0b615d23dbe8e997be6d6fafa0942713892e8 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Fri, 19 Dec 2025 12:00:35 -0300
Subject: [PATCH 7/9] Prevented multiple antenna instances to be created for
 Satellite Beamforming

---
 sharc/station_factory.py | 31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/sharc/station_factory.py b/sharc/station_factory.py
index 1f25d035..18bf86a3 100644
--- a/sharc/station_factory.py
+++ b/sharc/station_factory.py
@@ -1784,6 +1784,22 @@ def generate_mss_d2d(
             antenna_pattern = AntennaS1528Taylor(params.antenna.itu_r_s_1528)
         elif params.antenna.pattern == "MSS Adjacent":
             antenna_pattern = AntennaMSSAdjacent(params.frequency)
+        elif params.antenna.pattern == "Satellite Beamforming":
+            params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
+            params.antenna.itu_r_s_1528.antenna_gain = 34.6
+            params.antenna.itu_r_s_1528.antenna_3_dB_bw = 3.06
+            params.antenna.itu_r_s_1528.antenna_l_s = -35.0
+            params.antenna.itu_r_s_1528.far_out_side_lobe = -25
+            antenna_pattern_high = AntennaS1528(params.antenna.itu_r_s_1528)
+
+            # Different antenna parameters for low elevation beams
+            params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
+            params.antenna.itu_r_s_1528.antenna_gain = 40.0
+            params.antenna.itu_r_s_1528.antenna_3_dB_bw = 1.68
+            params.antenna.itu_r_s_1528.antenna_l_s = -20.0
+            params.antenna.itu_r_s_1528.far_out_side_lobe = -15
+            antenna_pattern_low = AntennaS1528(params.antenna.itu_r_s_1528)
+
         else:
             raise ValueError(
                 f"generate_mss_ss: Invalid antenna type: {params.antenna.pattern}")
@@ -1791,20 +1807,9 @@ def generate_mss_d2d(
         for i in range(mss_d2d.num_stations):
             if params.antenna.pattern == "Satellite Beamforming":
                 if beams_ground_elev[i] >= 50:
-                    params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
-                    params.antenna.itu_r_s_1528.antenna_gain = 34.6
-                    params.antenna.itu_r_s_1528.antenna_3_dB_bw = 3.06
-                    params.antenna.itu_r_s_1528.antenna_l_s = -35.0
-                    params.antenna.itu_r_s_1528.far_out_side_lobe = -25
-
+                    antenna_pattern = antenna_pattern_high
                 else:
-                    params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
-                    params.antenna.itu_r_s_1528.antenna_gain = 40.0
-                    params.antenna.itu_r_s_1528.antenna_3_dB_bw = 1.68
-                    params.antenna.itu_r_s_1528.antenna_l_s = -20.0
-                    params.antenna.itu_r_s_1528.far_out_side_lobe = -15
-
-                antenna_pattern = AntennaS1528(params.antenna.itu_r_s_1528)
+                    antenna_pattern = antenna_pattern_low
 
             mss_d2d.antenna[i] = antenna_pattern
 

From 425d52a2a44c084f79968e0b44a82eb1fe304f21 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Mon, 22 Dec 2025 10:37:10 -0300
Subject: [PATCH 8/9] update(antenna): Adding DC-MSS System4 antenna based on
 high and low elevation parameters.

---
 sharc/antenna/antenna_factory.py              |  6 +++
 .../antenna/parameters_antenna_system4.py     | 48 +++++++++++++++++++
 sharc/parameters/parameters_antenna.py        | 13 ++++-
 sharc/station_factory.py                      | 18 ++-----
 4 files changed, 68 insertions(+), 17 deletions(-)
 create mode 100644 sharc/parameters/antenna/parameters_antenna_system4.py

diff --git a/sharc/antenna/antenna_factory.py b/sharc/antenna/antenna_factory.py
index 051e1f33..5df02cb9 100644
--- a/sharc/antenna/antenna_factory.py
+++ b/sharc/antenna/antenna_factory.py
@@ -60,6 +60,12 @@ def create_antenna(
                     azimuth,
                     elevation
                 )
+            case "Antenna System 4":
+                # Handled in station_factory.py since it requires two antennas
+                raise NotImplementedError(
+                    "Antenna System 4 requires two antenna instances and "
+                    "should be created in station_factory.py"
+                )
             case _:
                 raise ValueError(
                     f"Antenna factory does not support pattern {
diff --git a/sharc/parameters/antenna/parameters_antenna_system4.py b/sharc/parameters/antenna/parameters_antenna_system4.py
new file mode 100644
index 00000000..b98109b9
--- /dev/null
+++ b/sharc/parameters/antenna/parameters_antenna_system4.py
@@ -0,0 +1,48 @@
+# Implementation of ParametersAntennaSystem4 class
+# The System 4 antenna was defined in WP4C Working Document 4C/356 from October 2025
+# The antenna is based on S.1528 recommends 1.2
+# It defines two sets of parameters for high elevation and low elevation beams
+from dataclasses import dataclass, field
+
+from sharc.parameters.parameters_base import ParametersBase
+from sharc.parameters.antenna.parameters_antenna_s1528 import ParametersAntennaS1528
+
+
+@dataclass
+class ParametersAntennaSystem4(ParametersBase):
+    """Dataclass containing the Antenna System 4 parameters for the simulator.
+    """
+    section_name: str = "Antenna System 4"
+
+    # Parameters for high elevation beams
+    antenna_parameters_high: ParametersAntennaS1528 = field(
+        # we don't care about frequency and bandwidth here. Just to make validation work.
+        default_factory=lambda: ParametersAntennaS1528(
+            frequency=-1,
+            bandwidth=-1
+        )
+    )
+
+    # Parameters for low elevation beams
+    antenna_parameters_low: ParametersAntennaS1528 = field(
+        # we don't care about frequency and bandwidth here. Just to make validation work.
+        default_factory=lambda: ParametersAntennaS1528(
+            frequency=-1,
+            bandwidth=-1
+        )
+    )
+
+    def load_parameters_from_file(self, config_file: str):
+        """Load the parameters from file an run a sanity check.
+
+        Parameters
+        ----------
+        file_name : str
+            the path to the configuration file
+
+        Raises
+        ------
+        ValueError
+            if a parameter is not valid
+        """
+        super().load_parameters_from_file(config_file)
\ No newline at end of file
diff --git a/sharc/parameters/parameters_antenna.py b/sharc/parameters/parameters_antenna.py
index a8fc8a0e..687483f4 100644
--- a/sharc/parameters/parameters_antenna.py
+++ b/sharc/parameters/parameters_antenna.py
@@ -5,6 +5,7 @@
 from sharc.parameters.antenna.parameters_antenna_s672 import ParametersAntennaS672
 from sharc.parameters.antenna.parameters_antenna_with_freq import ParametersAntennaWithFreq
 from sharc.parameters.imt.parameters_antenna_imt import ParametersAntennaImt
+from sharc.parameters.antenna.parameters_antenna_system4 import ParametersAntennaSystem4
 
 from dataclasses import dataclass, field
 import typing
@@ -29,7 +30,8 @@ class ParametersAntenna(ParametersBase):
         "ITU-R-S.1528-Taylor",
         "ITU-R-S.1528-Section1.2",
         "ITU-R-S.1528-LEO",
-        "MSS Adjacent"]
+        "MSS Adjacent",
+        "Antenna System 4"]
 
     # chosen antenna radiation pattern
     pattern: typing.Literal["OMNI",
@@ -44,7 +46,8 @@ class ParametersAntenna(ParametersBase):
                             "ITU-R-S.1528-Taylor",
                             "ITU-R-S.1528-Section1.2",
                             "ITU-R-S.1528-LEO",
-                            "MSS Adjacent"] = None
+                            "MSS Adjacent",
+                            "Antenna System 4"] = None
 
     # antenna gain [dBi]
     gain: float = None
@@ -94,6 +97,10 @@ class ParametersAntenna(ParametersBase):
         default_factory=ParametersAntennaS672,
     )
 
+    antenna_system_4: ParametersAntennaSystem4 = field(
+        default_factory=ParametersAntennaSystem4,
+    )
+
     def set_external_parameters(self, **kwargs):
         """
         Set external parameters for all sub-parameters of the antenna.
@@ -198,6 +205,8 @@ def validate(self, ctx):
                 self.itu_r_s_672.validate(f"{ctx}.itu_r_s_672")
             case "MSS Adjacent":
                 self.mss_adjacent.validate(f"{ctx}.mss_adjacent")
+            case "Antenna System 4":
+                self.antenna_system_4.validate(f"{ctx}.antenna_system_4")
             case _:
                 raise NotImplementedError(
                     "ParametersAntenna.validate does not implement this antenna validation!", )
diff --git a/sharc/station_factory.py b/sharc/station_factory.py
index 18bf86a3..1005d43b 100644
--- a/sharc/station_factory.py
+++ b/sharc/station_factory.py
@@ -1784,21 +1784,9 @@ def generate_mss_d2d(
             antenna_pattern = AntennaS1528Taylor(params.antenna.itu_r_s_1528)
         elif params.antenna.pattern == "MSS Adjacent":
             antenna_pattern = AntennaMSSAdjacent(params.frequency)
-        elif params.antenna.pattern == "Satellite Beamforming":
-            params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
-            params.antenna.itu_r_s_1528.antenna_gain = 34.6
-            params.antenna.itu_r_s_1528.antenna_3_dB_bw = 3.06
-            params.antenna.itu_r_s_1528.antenna_l_s = -35.0
-            params.antenna.itu_r_s_1528.far_out_side_lobe = -25
-            antenna_pattern_high = AntennaS1528(params.antenna.itu_r_s_1528)
-
-            # Different antenna parameters for low elevation beams
-            params.antenna.itu_r_s_1528.antenna_pattern = "ITU-R-S.1528-Section1.2"
-            params.antenna.itu_r_s_1528.antenna_gain = 40.0
-            params.antenna.itu_r_s_1528.antenna_3_dB_bw = 1.68
-            params.antenna.itu_r_s_1528.antenna_l_s = -20.0
-            params.antenna.itu_r_s_1528.far_out_side_lobe = -15
-            antenna_pattern_low = AntennaS1528(params.antenna.itu_r_s_1528)
+        elif params.antenna.pattern == "Antenna System 4":
+            antenna_pattern_high = AntennaS1528(params.antenna.antenna_system_4.antenna_parameters_high)
+            antenna_pattern_low = AntennaS1528(params.antenna.antenna_system_4.antenna_parameters_low)
 
         else:
             raise ValueError(

From 89def3791e7e522c6a8da9583cff80ed7eef76b4 Mon Sep 17 00:00:00 2001
From: Bruno Faria <brunohcfaria@gmail.com>
Date: Wed, 21 Jan 2026 12:01:05 -0300
Subject: [PATCH 9/9] fix(topology): Fixed beam ground elevation shape.

---
 sharc/parameters/antenna/parameters_antenna_system4.py | 2 +-
 sharc/topology/topology_imt_mss_dc.py                  | 4 ++++
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/sharc/parameters/antenna/parameters_antenna_system4.py b/sharc/parameters/antenna/parameters_antenna_system4.py
index b98109b9..ca96cb82 100644
--- a/sharc/parameters/antenna/parameters_antenna_system4.py
+++ b/sharc/parameters/antenna/parameters_antenna_system4.py
@@ -45,4 +45,4 @@ def load_parameters_from_file(self, config_file: str):
         ValueError
             if a parameter is not valid
         """
-        super().load_parameters_from_file(config_file)
\ No newline at end of file
+        super().load_parameters_from_file(config_file)
diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py
index 88917282..66b25139 100644
--- a/sharc/topology/topology_imt_mss_dc.py
+++ b/sharc/topology/topology_imt_mss_dc.py
@@ -592,6 +592,10 @@ def get_satellite_pointing(
             (total_active_satellites, orbit_params.num_beams)
         )
 
+        beams_ground_elev = beams_ground_elev.reshape(
+            (total_active_satellites, orbit_params.num_beams)
+        )
+
         # Rotate and set the each beam azimuth and elevation angles - only for
         # the visible satellites
         nadir_elev = all_nadir_elev[active_satellite_idxs]