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/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/sharc/parameters/antenna/parameters_antenna_system4.py b/sharc/parameters/antenna/parameters_antenna_system4.py new file mode 100644 index 00000000..ca96cb82 --- /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) 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/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/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 diff --git a/sharc/station_factory.py b/sharc/station_factory.py index 02384d4f..1005d43b 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, @@ -1783,11 +1784,21 @@ 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 == "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( 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: + antenna_pattern = antenna_pattern_high + else: + antenna_pattern = antenna_pattern_low + 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 bd0f8cbb..66b25139 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,147 @@ 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, beams_ground_elev, 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)) + 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) + 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,)) @@ -346,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) @@ -408,7 +418,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 +427,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,13 +460,18 @@ 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 beams_azim = [] beams_elev = [] + beams_ground_elev = [] n = 0 for act_sat in active_satellite_idxs: @@ -464,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([]) @@ -473,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( @@ -565,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) @@ -574,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] @@ -589,7 +611,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( diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml index 1cd0d9fc..86631172 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: @@ -998,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" @@ -1096,6 +1100,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..b959acfb 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( @@ -871,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) 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()