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/imt/parameters_imt_mss_dc.py b/sharc/parameters/imt/parameters_imt_mss_dc.py index c14079f4..4ad0d409 100644 --- a/sharc/parameters/imt/parameters_imt_mss_dc.py +++ b/sharc/parameters/imt/parameters_imt_mss_dc.py @@ -6,7 +6,11 @@ from sharc.parameters.parameters_base import ParametersBase from sharc.parameters.parameters_orbit import ParametersOrbit -from sharc.parameters.imt.parameters_grid import ParametersSatelliteWithServiceGrid, ParametersSelectActiveSatellite +from sharc.parameters.imt.parameters_grid import ( + ParametersSatelliteWithServiceGrid, + ParametersSelectActiveSatellite, + ParametersZone, +) SHARC_ROOT_DIR = (Path(__file__) / ".." / ".." / ".." / "..").resolve() @@ -189,6 +193,30 @@ def validate(self, ctx): ) +@dataclass +class ParametersPowerControlZone(ParametersBase): + """Dataclass for a power control zone in the IMT MSS-DC topology.""" + geometry: ParametersZone = field(default_factory=ParametersZone) + power_backoff_db: float = None + + +@dataclass +class ParametersPowerControl(ParametersBase): + """Dataclass for power control parameters in the IMT MSS-DC topology.""" + zones: list[ParametersPowerControlZone] = field( + default_factory=lambda: [ParametersPowerControlZone()]) + + def validate(self, ctx): + """ + Validate the power control parameters. + """ + super().validate(ctx) + for i in range(len(self.zones)): + self.zones[i].geometry.validate(ctx + f"zones.{i}.geometry") + if not isinstance(self.zones[i].power_backoff_db, float): + raise ValueError("power_backoff_db is not properly defined.") + + @dataclass class ParametersImtMssDc(ParametersBase): """Dataclass for the IMT MSS-DC topology parameters.""" @@ -211,6 +239,9 @@ class ParametersImtMssDc(ParametersBase): # for IMT Space Stations beam_radius: float = 36516.0 + power_control_zones: ParametersPowerControl = field( + default_factory=ParametersPowerControl) + sat_is_active_if: ParametersSelectActiveSatellite = field( default_factory=ParametersSelectActiveSatellite) 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/parameters/parameters_mss_d2d.py b/sharc/parameters/parameters_mss_d2d.py index c566269c..7b0e5961 100644 --- a/sharc/parameters/parameters_mss_d2d.py +++ b/sharc/parameters/parameters_mss_d2d.py @@ -2,7 +2,11 @@ from dataclasses import dataclass, field, asdict from sharc.parameters.parameters_base import ParametersBase from sharc.parameters.parameters_orbit import ParametersOrbit -from sharc.parameters.imt.parameters_imt_mss_dc import ParametersSelectActiveSatellite, ParametersSectorPositioning +from sharc.parameters.imt.parameters_imt_mss_dc import ( + ParametersSelectActiveSatellite, + ParametersSectorPositioning, + ParametersPowerControl +) from sharc.parameters.parameters_p619 import ParametersP619 from sharc.parameters.parameters_antenna import ParametersAntenna from sharc.parameters.antenna.parameters_antenna_s1528 import ParametersAntennaS1528 @@ -40,6 +44,9 @@ class ParametersMssD2d(ParametersBase): # active beams_load_factor: float = 1.0 + power_control_zones: ParametersPowerControl = field( + default_factory=ParametersPowerControl) + # Central beam positioning beam_positioning: ParametersSectorPositioning = field( default_factory=ParametersSectorPositioning) 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/simulation_downlink.py b/sharc/simulation_downlink.py index 58e9231b..f39b791d 100644 --- a/sharc/simulation_downlink.py +++ b/sharc/simulation_downlink.py @@ -152,20 +152,23 @@ def power_control(self): """ # Currently, the maximum transmit power of the base station is equaly # divided among the selected UEs - total_power = self.parameters.imt.bs.conducted_power \ + total_power = self.bs.tx_power \ + self.bs_power_gain tx_power = total_power - 10 * math.log10(self.parameters.imt.ue.k) # calculate transmit powers to have a structure such as # {bs_1: [pwr_1, pwr_2,...], ...}, where bs_1 is the base station id, # pwr_1 is the transmit power from bs_1 to ue_1, pwr_2 is the transmit # power from bs_1 to ue_2, etc + # NOTE: Re-writing bs.tx_power here as dict bs_active = np.where(self.bs.active)[0] self.bs.tx_power = dict( - [(bs, tx_power * np.ones(self.parameters.imt.ue.k)) for bs in bs_active]) + [(bs, tx_power[bs] * np.ones(self.parameters.imt.ue.k)) for bs in bs_active]) # Update the spectral mask if self.adjacent_channel: - self.bs.spectral_mask.set_mask(p_tx=total_power) + # FIXME: assumes all BS have same total power + max_tx_power = self.parameters.imt.bs.conducted_power + self.bs_power_gain + self.bs.spectral_mask.set_mask(p_tx=max_tx_power) def calculate_sinr(self): """ @@ -268,8 +271,6 @@ def calculate_sinr_ext(self): # Inteferer transmit power in dBm over the overlapping band # (MHz) with UEs. if self.overlapping_bandwidth > 0: - # in_band_interf_power = self.param_system.tx_power_density + \ - # 10 * np.log10(self.overlapping_bandwidth * 1e6) + 30 with warnings.catch_warnings(): warnings.filterwarnings( "ignore", @@ -277,7 +278,7 @@ def calculate_sinr_ext(self): message="divide by zero encountered in log10", ) in_band_interf_power = \ - self.param_system.tx_power_density + 10 * np.log10( + self.system.tx_power_density[system_interfering] + 10 * np.log10( self.ue.bandwidth[ue] * 1e6 ) + 10 * np.log10(weights) - \ self.coupling_loss_imt_system[ue, :][system_interfering] @@ -286,12 +287,12 @@ def calculate_sinr_ext(self): if self.adjacent_channel: # emissions outside of tx bandwidth and inside of rx bw # due to oob emissions on tx side - tx_oob = np.resize(-500., 1) + tx_oob = np.ones_like(system_interfering) * -500. # emissions outside of rx bw and inside of tx bw # due to non ideal filtering on rx side # will be the same for all UE's, only considering - rx_oob = np.resize(-500., 1) + rx_oob = np.ones_like(system_interfering) * -500. # TODO: M.2101 states that: # "The ACIR value should be calculated based on per UE allocated number of resource blocks" @@ -320,7 +321,8 @@ def calculate_sinr_ext(self): # NOTE: only the power not overlapping is attenuated by ACS # tx_pow_adj_lin = PSD * non_overlap_imt_bw # rx_oob = tx_pow_adj_lin / acs - rx_oob[::] = self.param_system.tx_power_density + 10 * np.log10(non_overlap_sys_bw * 1e6) - acs_dB + rx_oob[::] = self.system.tx_power_density[system_interfering] + \ + 10 * np.log10(non_overlap_sys_bw * 1e6) - acs_dB elif self.parameters.imt.adjacent_ch_reception == "OFF": pass else: @@ -337,7 +339,7 @@ def calculate_sinr_ext(self): warnings.filterwarnings("ignore", category=RuntimeWarning, message="divide by zero encountered in log10") - tx_oob[0] = self.system.spectral_mask.power_calc( + tx_oob[::] = self.system.spectral_mask.power_calc( center_freqs, ue_bws ) - 30 @@ -365,8 +367,7 @@ def calculate_sinr_ext(self): # tx_oob_in_measurement = (tx_pow_lin / aclr) # => approx. PSD = (tx_pow_lin / aclr) / measurement_bw # approximated received tx_oob = PSD * non_overlap_imt_bw - tx_oob[::] = self.param_system.tx_power_density + \ - 10 * np.log10(1e6) - \ + tx_oob[::] = self.system.tx_power_density[system_interfering] + 60 - \ aclr_dB + 10 * np.log10( non_overlap_imt_bw) elif self.param_system.adjacent_ch_emissions == "OFF": @@ -377,9 +378,9 @@ def calculate_sinr_ext(self): self.param_system.adjacent_ch_emissions}") if self.param_system.adjacent_ch_emissions != "OFF": - tx_oob = tx_oob[:] - self.coupling_loss_imt_system[ue, :][system_interfering] + tx_oob = tx_oob - self.coupling_loss_imt_system[ue, system_interfering] - rx_oob = rx_oob[:] - self.coupling_loss_imt_system_adjacent[ue, :][system_interfering] + rx_oob = rx_oob - self.coupling_loss_imt_system_adjacent[ue, system_interfering] # Out of band power # sum linearly power leaked into band and power received in the @@ -388,9 +389,9 @@ def calculate_sinr_ext(self): 10 ** (0.1 * tx_oob) + 10 ** (0.1 * rx_oob) ) # Total external interference into the UE in dBm - ue_ext_int = 10 * np.log10(np.power(10, - 0.1 * in_band_interf_power) + np.power(10, - 0.1 * oob_power)) + ue_ext_int = 10 * np.log10( + np.power(10, 0.1 * in_band_interf_power) + + np.power(10, 0.1 * oob_power)) # Sum all the interferers for each UE self.ue.ext_interference[ue] = 10 * \ @@ -415,7 +416,8 @@ def calculate_sinr_ext(self): self.ue.geom) # shape: [n_tx, n_ue] # EIRP in dBW/MHz per transmitter - eirp_dBW_MHz = self.param_system.tx_power_density + \ + # FIXME: Should be calculated only for active transmitters + eirp_dBW_MHz = self.system.tx_power_density[:, np.newaxis] + \ 60 + self.system_imt_antenna_gain # PFD formula (dBW/m²/MHz) diff --git a/sharc/simulation_uplink.py b/sharc/simulation_uplink.py index ebb30718..f904a34f 100644 --- a/sharc/simulation_uplink.py +++ b/sharc/simulation_uplink.py @@ -237,7 +237,7 @@ def calculate_sinr_ext(self): warnings.filterwarnings("ignore", category=RuntimeWarning, message="divide by zero encountered in log10") - in_band_interf = self.param_system.tx_power_density + \ + in_band_interf = self.system.tx_power_density[system_interfering] + \ 10 * np.log10(beams_bw[:, np.newaxis] * 1e6) + \ 10 * np.log10(weights)[:, np.newaxis] - \ self.coupling_loss_imt_system[active_beams, :][:, system_interfering] @@ -247,11 +247,11 @@ def calculate_sinr_ext(self): if self.adjacent_channel: # emissions outside of tx bandwidth and inside of rx bw # due to oob emissions on tx side - tx_oob = np.resize(-500., len(active_beams)) + tx_oob = np.ones((len(active_beams), len(system_interfering))) * -500. # emissions outside of rx bw and inside of tx bw # due to non ideal filtering on rx side - rx_oob = np.resize(-500., len(active_beams)) + rx_oob = np.ones((len(active_beams), len(system_interfering))) * -500. # NOTE: M.2101 states that: # "The ACIR value should be calculated based on per UE allocated number of resource blocks" @@ -265,7 +265,8 @@ def calculate_sinr_ext(self): ) self._acs_warned = True acs_dB = self.parameters.imt.bs.adjacent_ch_selectivity - rx_oob[::] = self.param_system.tx_power_density + 10 * np.log10(non_overlap_sys_bw * 1e6) - acs_dB + rx_oob[:] = self.system.tx_power_density[system_interfering] + \ + 10 * np.log10(non_overlap_sys_bw * 1e6) - acs_dB elif self.parameters.imt.adjacent_ch_reception == "OFF": pass elif self.parameters.imt.adjacent_ch_reception is False: @@ -281,11 +282,12 @@ def calculate_sinr_ext(self): warnings.filterwarnings("ignore", category=RuntimeWarning, message="divide by zero encountered in log10") + # FIXME: given that each system station may have different power, this may not be correct. for i, center_freq, bw in zip( range(len(self.bs.center_freq[bs])), self.bs.center_freq[bs], beams_bw): # mask returns dBm # so we convert to [dB] - tx_oob[i] = self.system.spectral_mask.power_calc( + tx_oob[i, :] = self.system.spectral_mask.power_calc( center_freq, bw ) - 30 @@ -310,10 +312,8 @@ def calculate_sinr_ext(self): ) # [dB] - tx_oob[::] = self.param_system.tx_power_density + \ - 10 * np.log10(1e6) - \ - aclr_dB + 10 * np.log10( - non_overlap_imt_bw) + tx_oob[:] = self.system.tx_power_density[system_interfering] + 60 - \ + aclr_dB + 10 * np.log10(non_overlap_imt_bw)[:, np.newaxis] elif self.param_system.adjacent_ch_emissions == "OFF": pass else: @@ -323,10 +323,10 @@ def calculate_sinr_ext(self): if self.param_system.adjacent_ch_emissions != "OFF": # oob for system is inband for IMT - tx_oob = tx_oob[:, np.newaxis] - self.coupling_loss_imt_system[active_beams, :][:, system_interfering] + tx_oob = tx_oob - self.coupling_loss_imt_system[active_beams, :][:, system_interfering] # oob for IMT - rx_oob = rx_oob[:, np.newaxis] - self.coupling_loss_imt_system_adjacent[active_beams, :][:, system_interfering] + rx_oob = rx_oob - self.coupling_loss_imt_system_adjacent[active_beams, :][:, system_interfering] # Out of band power # sum linearly power leaked into band and power received in the @@ -337,7 +337,6 @@ def calculate_sinr_ext(self): # [dBm] ext_interference = 10 * np.log10(in_band_interf_lin + oob_interf_lin) + 30 - # Sum all the interferers from each active system transmitters for each bs self.bs.ext_interference[bs] = 10 * np.log10( np.sum(np.power(10, 0.1 * ext_interference), axis=1)) diff --git a/sharc/station_factory.py b/sharc/station_factory.py index 02384d4f..775894b8 100644 --- a/sharc/station_factory.py +++ b/sharc/station_factory.py @@ -59,6 +59,7 @@ from sharc.mask.spectral_mask_3gpp import SpectralMask3Gpp from sharc.mask.spectral_mask_mss import SpectralMaskMSS from sharc.support.sharc_geom import CoordinateSystem +from sharc.support.geometry import SimulatorGeometry from sharc.support.sharc_utils import wrap2_180 @@ -96,6 +97,7 @@ def generate_imt_base_stations( num_bs = topology.num_base_stations imt_base_stations = StationManager(num_bs) imt_base_stations.station_type = StationType.IMT_BS + power_backoff = 0.0 if param.topology.type == "NTN": imt_base_stations.geom.set_global_coords( topology.space_station_x * np.ones(num_bs), @@ -112,6 +114,7 @@ def generate_imt_base_stations( elev=topology.elevation ) imt_base_stations.is_space_station = True + power_backoff = topology.power_backoff else: if topology.determines_local_geometry: imt_base_stations.geom = topology.get_bs_geometry() @@ -132,7 +135,9 @@ def generate_imt_base_stations( imt_base_stations.active = random_number_gen.rand( num_bs, ) < param.bs.load_probability - imt_base_stations.tx_power = param.bs.conducted_power * np.ones(num_bs) + # Conducted power per antenna element. Total power will depend on + # the number of antenna elements and it's configured in power control. + imt_base_stations.tx_power = param.bs.conducted_power * np.ones(num_bs) - power_backoff imt_base_stations.rx_power = dict( [(bs, -500 * np.ones(param.ue.k)) for bs in range(num_bs)], ) @@ -869,10 +874,7 @@ def generate_single_space_station( ) space_station.active = np.array([True]) - space_station.tx_power = np.array( - [param.tx_power_density + 10 * - math.log10(param.bandwidth * 1e6) + 30], - ) + space_station.tx_power_density = param.tx_power_density * np.ones(space_station.num_stations) space_station.rx_interference = -500 space_station.antenna = np.array([ @@ -944,10 +946,7 @@ def generate_fss_space_station(param: ParametersFssSs): ) fss_space_station.active = np.array([True]) - fss_space_station.tx_power = np.array( - [param.tx_power_density + 10 * - math.log10(param.bandwidth * 1e6) + 30], - ) + fss_space_station.tx_power_density = np.ones(fss_space_station.num_stations) * param.tx_power_density fss_space_station.rx_interference = np.array([-500]) if param.antenna_pattern == "OMNI": @@ -1079,10 +1078,7 @@ def generate_fss_earth_station( ) fss_earth_station.active = np.array([True]) - fss_earth_station.tx_power = np.array( - [param.tx_power_density + 10 * - math.log10(param.bandwidth * 1e6) + 30], - ) + fss_earth_station.tx_power_density = np.ones(fss_earth_station.num_stations) * param.tx_power_density fss_earth_station.rx_interference = np.array([-500]) if param.antenna_pattern.upper() == "OMNI": @@ -1253,10 +1249,8 @@ def generate_single_earth_station( single_earth_station.active = np.array([True]) single_earth_station.bandwidth = np.array([param.bandwidth]) - single_earth_station.tx_power = np.array( - [param.tx_power_density + 10 * - math.log10(param.bandwidth * 1e6) + 30], - ) + single_earth_station.tx_power_density = np.ones( + single_earth_station.num_stations) * param.tx_power_density single_earth_station.noise_temperature = np.array( [param.noise_temperature]) @@ -1317,10 +1311,7 @@ def generate_fs_station(param: ParametersFs): ) fs_station.active = np.array([True]) - fs_station.tx_power = np.array( - [param.tx_power_density + 10 * - math.log10(param.bandwidth * 1e6) + 30], - ) + fs_station.tx_power_density = np.ones(fs_station.num_stations) * param.tx_power_density fs_station.rx_interference = -500 if param.antenna_pattern == "OMNI": @@ -1386,6 +1377,8 @@ def generate_haps( haps.antenna = np.empty(num_haps, dtype=Antenna) + haps.tx_power_density = np.ones(haps.num_stations) * param.tx_power_density + if param.antenna_pattern == "OMNI": for i in range(num_haps): haps.antenna[i] = AntennaOmni(param.antenna_gain) @@ -1458,6 +1451,7 @@ def generate_rns( ) sys.exit(1) + rns.tx_power_density = np.ones(rns.num_stations) * param.tx_power_density rns.bandwidth = np.array([param.bandwidth]) rns.noise_temperature = param.noise_temperature rns.thermal_noise = -500 @@ -1569,6 +1563,8 @@ def generate_space_station( space_station.active = np.array([True]) space_station.rx_interference = np.array([-500]) + space_station.tx_power_density = np.ones(space_station.num_stations) * param.tx_power_density + if param.antenna_pattern == "OMNI": space_station.antenna = np.array([AntennaOmni(param.antenna_gain)]) elif param.antenna_pattern == "ITU-R RS.1813": @@ -1635,9 +1631,9 @@ def generate_mss_ss(param_mss: ParametersMssSs): mss_ss.is_space_station = True mss_ss.active = np.ones(num_bs, dtype=int) - mss_ss.tx_power = np.ones( - num_bs, dtype=int) * param_mss.tx_power_density + 10 * np.log10( - param_mss.bandwidth * 10**6) + mss_ss.tx_power_density = np.ones( + num_bs, dtype=int) * param_mss.tx_power_density + mss_ss.antenna = np.empty(num_bs, dtype=AntennaS1528Leo) for i in range(num_bs): @@ -1718,6 +1714,14 @@ def generate_mss_d2d( mss_d2d = StationManager(n=total_satellites) mss_d2d.station_type = StationType.MSS_D2D # Set the station type to MSS D2D mss_d2d.is_space_station = True # Indicate that the station is in space + mss_d2d.geom = SimulatorGeometry( + mss_d2d.num_stations, True, + ( + coordinate_system.ref_lat, + coordinate_system.ref_long, + coordinate_system.ref_alt, + ) + ) if params.spectral_mask == "IMT-2020": mss_d2d.spectral_mask = SpectralMaskImt(StationType.IMT_BS, @@ -1738,20 +1742,33 @@ def generate_mss_d2d( else: raise ValueError( f"Invalid or not implemented spectral mask - {params.spectral_mask}") - mss_d2d.spectral_mask.set_mask( - params.tx_power_density + - 10 * - np.log10( - params.bandwidth * - 1e6) + 30 - ) - # Configure satellite positions in the StationManager + max_tx_power = params.tx_power_density + 10 * np.log10(params.bandwidth * 1e6) + + # Configure satellite positions in the StationManager x = mss_d2d_values["sat_x"] y = mss_d2d_values["sat_y"] 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"] + power_backoff = mss_d2d_values["sat_power_backoff"] + + mss_d2d.tx_power_density = params.tx_power_density * np.ones( + mss_d2d.num_stations, + ) - power_backoff + + mss_d2d.spectral_mask.set_mask(max_tx_power + 30) + + sat_lat = mss_d2d_values["sat_lat"] + sat_lon = mss_d2d_values["sat_lon"] + + mss_d2d.geom.set_local_coord_sys( + sat_lat, + sat_lon, + # put local coord system in subsatellite + np.zeros_like(sat_lon) + ) mss_d2d.geom.set_global_coords( x, y, z, azim, elev, @@ -1783,11 +1800,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 == "Antenna System 4": + 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/station_manager.py b/sharc/station_manager.py index 07ecd3ca..1e8e982a 100644 --- a/sharc/station_manager.py +++ b/sharc/station_manager.py @@ -22,34 +22,34 @@ class StationManager(object): """ def __init__(self, n): - self.num_stations = n - self.idx_orbit = np.empty(n) - self.indoor = np.zeros(n, dtype=bool) - self.active = np.ones(n, dtype=bool) - self.tx_power = np.empty(n) - self.rx_power = np.empty(n) - self.rx_interference = np.empty(n) # Rx interferece in dBW - self.ext_interference = np.empty(n) # External interferece in dBW - self.antenna = np.empty(n, dtype=Antenna) + self.num_stations = n # Number of stations managed + self.idx_orbit = np.empty(n) # Index of the orbit for space stations + self.indoor = np.zeros(n, dtype=bool) # Whether the station is indoor + self.active = np.ones(n, dtype=bool) # Whether the station is active in the simulation + self.tx_power = np.zeros(n) # Station transmit power - dBm for IMT + self.tx_power_density = np.empty(n) # Transmit power density in dB/Hz + self.rx_power = np.empty(n) # Received power + self.rx_interference = np.empty(n) # Rx interferece + self.ext_interference = np.empty(n) # External interferece + self.antenna = np.empty(n, dtype=Antenna) # antenna objects, one for each station self.bandwidth = np.empty(n) # Bandwidth in MHz - self.noise_figure = np.empty(n) - self.noise_temperature = np.empty(n) - self.thermal_noise = np.empty(n) - self.total_interference = np.empty(n) + self.noise_figure = np.empty(n) # Noise figure in dB + self.noise_temperature = np.empty(n) # Noise temperature in K + self.thermal_noise = np.empty(n) # Thermal noise in dBm/MHz + self.total_interference = np.empty(n) # Total received interference self.pfd_external = np.empty(n) # External PFD in dBW/m²/MHz - # Aggregated External PFD in dBW/m²/MHz - self.pfd_external_aggregated = np.empty(n) - self.snr = np.empty(n) - self.sinr = np.empty(n) - self.sinr_ext = np.empty(n) - self.inr = np.empty(n) # INR in dBm/MHz - self.pfd = np.empty(n) # Powerflux density in dBm/m^2 - self.spectral_mask = np.empty(n, dtype=SpectralMask) - self.center_freq = np.empty(n) - self.station_type = StationType.NONE - self.is_space_station = False - self.geom = SimulatorGeometry(n) - self.max_earth_sta_interf_distance = np.inf + self.pfd_external_aggregated = np.empty(n) # Aggregated External PFD in dBW/m²/MHz + self.snr = np.empty(n) # SNR in dB + self.sinr = np.empty(n) # SINR in dB + self.sinr_ext = np.empty(n) # Used to store IMT SINR with external interference in dB + self.inr = np.empty(n) # INR in dB + self.pfd = np.empty(n) # Power flux density in dBm/m^2/MHz + self.spectral_mask = np.empty(n, dtype=SpectralMask) # Spectral mask objects, one for each station + self.center_freq = np.empty(n) # Center frequency in MHz + self.station_type = StationType.NONE # Station type + self.is_space_station = False # Whether the station is a space station - used for path loss calculations + self.geom = SimulatorGeometry(n) # Geometry object to manage station coordinates + self.max_earth_sta_interf_distance = np.inf # Max distance for earth station interference calculations def get_station_list(self, id=None) -> list: """Return a list of Station objects for the given indices. diff --git a/sharc/topology/topology_imt_mss_dc.py b/sharc/topology/topology_imt_mss_dc.py index bd0f8cbb..32cda473 100644 --- a/sharc/topology/topology_imt_mss_dc.py +++ b/sharc/topology/topology_imt_mss_dc.py @@ -63,6 +63,7 @@ def __init__(self, params: ParametersImtMssDc, self.space_station_x = None self.space_station_y = None self.space_station_z = None + self.power_backoff = None # power backoff per beam self.lat = None self.lon = None @@ -75,117 +76,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,126 +182,158 @@ 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_power_backoff, 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_power_backoff = np.array( + functools.reduce( + lambda x, + y: list(x) + + list(y), + beams_power_backoff)) + + 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,)) assert (space_station_z.shape == (num_base_stations,)) + assert (space_station_power_backoff.shape == (num_base_stations,)) assert (lat.shape == (num_base_stations,)) assert (lon.shape == (num_base_stations,)) assert (altitudes.shape == (num_base_stations,)) @@ -338,6 +355,7 @@ def get_coordinates( "num_satellites": num_base_stations, "num_active_satellites": len(active_satellite_idxs), "active_satellites_idxs": active_satellite_idxs, + "sat_power_backoff": space_station_power_backoff, "sat_x": space_station_x, "sat_y": space_station_y, "sat_z": space_station_z, @@ -346,6 +364,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) @@ -383,6 +402,26 @@ def get_satellite_pointing( grid = orbit_params.beam_positioning.service_grid.lon_lat_grid grid_lon = grid[0] grid_lat = grid[1] + + # create points(lon, lat) to compare to country + grid_points = gpd.points_from_xy( + grid_lon, grid_lat, crs=EARTH_DEFAULT_CRS) + + power_backoff_zones = orbit_params.power_control_zones.zones + + grid_power_backoffs = np.zeros_like(grid_lon) + # the first definition takes precedence + # so it should be executed last + for zone in reversed(power_backoff_zones): + mask = grid_points.within( + zone.geometry._polygon + ) + grid_power_backoffs[mask] = zone.power_backoff_db + + # pb = grid_power_backoffs[grid_power_backoffs!= 0.] + # print("pb", pb) + # exit() + grid_ecef = orbit_params.beam_positioning.service_grid.ecef_grid grid_x = grid_ecef[0] grid_y = grid_ecef[1] @@ -408,7 +447,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 +456,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,30 +489,41 @@ 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 + power_backoff = [] for act_sat in active_satellite_idxs: if act_sat in sat_points_towards: 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]) + power_backoff.append(grid_power_backoffs[sat_points_towards[act_sat]]) else: beams_azim.append([]) beams_elev.append([]) + beams_ground_elev.append([]) + power_backoff.append([]) # FIXME: change either this or the transform_ue_xyz to make this correct # we don't currently care sx = np.zeros(n) sy = np.zeros(n) - return beams_elev, beams_azim, sx, sy + return power_backoff, 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 +615,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 +625,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 +644,7 @@ def get_satellite_pointing( nadir_azim[i] ) - return beams_elev, beams_azim, sx, sy + return np.zeros_like(beams_elev), beams_elev, beams_azim, beams_ground_elev, sx, sy @staticmethod def get_distr( @@ -656,6 +711,9 @@ def calculate_coordinates(self, random_number_gen=np.random.RandomState()): self.y = sat_values["sectors_y"] self.z = sat_values["sectors_z"] + # power backoff per beam + self.power_backoff = sat_values["sat_power_backoff"] + self.indoor = np.zeros( self.num_base_stations, dtype=bool) # ofcourse, all are outdoor diff --git a/tests/parameters/parameters_for_testing.yaml b/tests/parameters/parameters_for_testing.yaml index 1cd0d9fc..51f257d1 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: @@ -282,6 +283,15 @@ imt: # if positive, makes border smaller by x km # if negative, makes border bigger by x km margin_from_border: 11.1241 + power_control_zones: + zones: + - power_backoff_db: 0.0 + geometry: + type: FROM_COUNTRIES + from_countries: + country_names: + - Brazil + margin_from_border: 150.0 # Parameters for the orbits orbits: - n_planes: 20 @@ -998,6 +1008,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" @@ -1020,6 +1033,22 @@ mss_d2d: tx_power_density: -30 # Number of sectors num_sectors: 19 + power_control_zones: + zones: + - geometry: + type: CIRCLE + circle: + center_lat: 0 + center_lon: 0 + radius_km: 10 + power_backoff_db: 13.0 + - geometry: + type: CIRCLE + circle: + center_lat: 0 + center_lon: 0 + radius_km: 15.5 + power_backoff_db: 5.0 beam_positioning: # type may be one of # "ANGLE_FROM_SUBSATELLITE", "ANGLE_AND_DISTANCE_FROM_SUBSATELLITE", @@ -1096,6 +1125,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..879818b2 100644 --- a/tests/parameters/test_parameters.py +++ b/tests/parameters/test_parameters.py @@ -30,9 +30,7 @@ def setUp(self): self.parameters.read_params() def test_parameters_imt(self): - """Unit test for ParametersIMT.""" - """Unit test for ParametersIMT - """ + """Unit test for ParametersIMT""" self.assertEqual(self.parameters.imt.topology.type, "INDOOR") self.assertEqual( self.parameters.imt.minimum_separation_distance_bs_ue, 1.3, @@ -287,6 +285,39 @@ def test_parameters_imt(self): ["LAT_LONG_INSIDE_COUNTRY", "MINIMUM_ELEVATION_FROM_ES", "MAXIMUM_ELEVATION_FROM_ES"], ) + # power_control_zones: + # zones: + # - power_backoff_db: 0.0 + # geometry: + # type: FROM_COUNTRIES + # from_countries: + # country_names: + # - Brazil + # margin_from_border: 150.0 + # - power_backoff_db: 0.0 + # geometry: + # type: FROM_COUNTRIES + # from_countries: + # country_names: + # - Brazil + # margin_from_border: 0 + # Test power control zones parameters + self.assertEqual( + self.parameters.imt.topology.mss_dc.power_control_zones.zones[0].power_backoff_db, + 0.0) + self.assertEqual( + self.parameters.imt.topology.mss_dc.power_control_zones.zones[0].geometry.type, + "FROM_COUNTRIES" + ) + self.assertEqual( + self.parameters.imt.topology.mss_dc.power_control_zones.zones[0].geometry.from_countries.country_names[0], + "Brazil" + ) + self.assertEqual( + self.parameters.imt.topology.mss_dc.power_control_zones.zones[0].geometry.from_countries.margin_from_border, + 150. + ) + self.parameters.imt.topology.mss_dc.validate("test_imt") self.assertEqual( self.parameters.imt.topology.mss_dc.beam_positioning.service_grid.beam_radius, @@ -313,6 +344,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 +742,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( @@ -740,6 +775,50 @@ def test_parametes_mss_d2d(self): self.assertEqual( self.parameters.mss_d2d.sat_is_active_if.maximum_elevation_from_es, 1.113) + self.assertEqual( + len(self.parameters.mss_d2d.power_control_zones.zones), + 2 + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[0].geometry.type, + "CIRCLE" + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[0].geometry.circle.radius_km, + 10. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[0].geometry.circle.center_lat, + 0. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[0].geometry.circle.center_lon, + 0. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[0].power_backoff_db, + 13. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[1].geometry.type, + "CIRCLE" + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[1].geometry.circle.radius_km, + 15.5 + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[1].geometry.circle.center_lat, + 0. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[1].geometry.circle.center_lon, + 0. + ) + self.assertEqual( + self.parameters.mss_d2d.power_control_zones.zones[1].power_backoff_db, + 5.0 + ) self.assertEqual(self.parameters.mss_d2d.param_p619.season, 'SUMMER') self.assertTrue(isinstance(self.parameters.mss_d2d.orbits, list)) expected_orbit_params = [ @@ -871,6 +950,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..3ccf5979 100644 --- a/tests/test_topology_imt_mss_dc.py +++ b/tests/test_topology_imt_mss_dc.py @@ -2,10 +2,11 @@ import numpy as np import numpy.testing as npt from sharc.topology.topology_imt_mss_dc import TopologyImtMssDc -from sharc.parameters.imt.parameters_imt_mss_dc import ParametersImtMssDc +from sharc.parameters.imt.parameters_imt_mss_dc import ParametersImtMssDc, ParametersPowerControlZone 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,54 @@ 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 + power_control_zone = ParametersPowerControlZone() + power_control_zone.power_backoff_db = 0.0 + self.params.power_control_zones.zones = [power_control_zone] + 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()