From 67ac410ce84854b9395913e4170eea1e531450c3 Mon Sep 17 00:00:00 2001 From: artistrea Date: Mon, 3 Nov 2025 23:42:16 -0300 Subject: [PATCH 1/8] 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 Date: Wed, 5 Nov 2025 15:24:02 -0300 Subject: [PATCH 2/8] 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 Date: Mon, 10 Nov 2025 17:29:46 -0300 Subject: [PATCH 3/8] 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 Date: Tue, 2 Dec 2025 13:26:07 -0300 Subject: [PATCH 4/8] 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 Date: Thu, 4 Dec 2025 10:52:40 -0300 Subject: [PATCH 5/8] 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 ") - 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 ") - 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 faa33b5f5302d4fc36c58df4cd43b6a03de798e6 Mon Sep 17 00:00:00 2001 From: Bruno Faria <33635129+brunohcfaria@users.noreply.github.com> Date: Thu, 29 Jan 2026 11:30:17 -0300 Subject: [PATCH 6/8] Fixed typo in season error message. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- sharc/propagation/atmosphere.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sharc/propagation/atmosphere.py b/sharc/propagation/atmosphere.py index d4d382fb..9d520803 100644 --- a/sharc/propagation/atmosphere.py +++ b/sharc/propagation/atmosphere.py @@ -387,7 +387,7 @@ 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}.") + raise ValueError(f"Invalid season name {season}.") if latitude <= 22: # low latitude From 6d95bc52311aa977b37932945fc4e3440d8a3dce Mon Sep 17 00:00:00 2001 From: Bruno Faria <33635129+brunohcfaria@users.noreply.github.com> Date: Thu, 29 Jan 2026 11:31:40 -0300 Subject: [PATCH 7/8] "visilibity" in the test docstring is misspelled "visilibity" in the test docstring is misspelled; it should be "visibility" to match the parameter name and improve clarity. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- tests/test_topology_imt_mss_dc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_topology_imt_mss_dc.py b/tests/test_topology_imt_mss_dc.py index cfa669d9..69a715ca 100644 --- a/tests/test_topology_imt_mss_dc.py +++ b/tests/test_topology_imt_mss_dc.py @@ -158,7 +158,7 @@ def test_visible_satellites(self): npt.assert_array_less(min_elevation_angle, xy_plane_elevations) def test_minimum_service_angle(self): - """Test minimum visilibity angle for service grid service.""" + """Test minimum visibility angle for service grid service.""" orbit = ParametersOrbit( n_planes=1, sats_per_plane=1, From e2b9011108b490ddf0630331f55fe0e91e7a4a0d Mon Sep 17 00:00:00 2001 From: Bruno Faria <33635129+brunohcfaria@users.noreply.github.com> Date: Thu, 29 Jan 2026 11:32:50 -0300 Subject: [PATCH 8/8] Update tests/parameters/test_parameters.py These two consecutive assertions both check major_minor_axis_ratio against 1.0, which is redundant and makes the test harder to maintain; one of them can be removed without changing test behavior. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- tests/parameters/test_parameters.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tests/parameters/test_parameters.py b/tests/parameters/test_parameters.py index b959acfb..e5fd967d 100644 --- a/tests/parameters/test_parameters.py +++ b/tests/parameters/test_parameters.py @@ -884,11 +884,6 @@ def test_parameters_single_space_station(self): 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